Solution for Part (b) of the 3rd Question of Section 3.3 in "Introduction to Probability" by Grinstead and Snell

riffle_problem.py

import numpy as np

def riffle(deck, n_times=1):
	for i in range(n_times):
		deck2 = []
		split = np.random.binomial(len(deck), 1/2, 1)[0]
		half1 = deck[:split]
		half2 = deck[split:]
		
		while len(half1) + len(half2) > 0:
			p = np.random.rand()
			weighted = len(half1) / (len(half1) + len(half2))
			if p < weighted:		
				try:
					deck2.append(half1.pop(0))
				except:
					pass
			else:			
				try:
					deck2.append(half2.pop(0))
				except:
					pass
				
		deck = deck2
		
	return deck

def prob_3(after, trials=10000):
	deck = [i for i in range(52)]
	success = 0
	for j in range(trials):
		deck2 = riffle(deck.copy(), after)
		if deck2.index(0) < 26:
			success += 1
	
	success /= trials
	
	return success

print("Probability after 3 riffles:", prob_3(3, 100000))
print("Probability after 7 riffles:", prob_3(7, 100000))