In [1]:

    

import numpy as np
def roll_me(N=100):
    rolls = 1 # Always have to roll at least once
    r1 = np.random.randint(1, N+1) # Note that randint uses the interval [low, high)
    r2 = r1
    while r2 >= r1:
        r1 = r2
        r2 = np.random.randint(1, N+1)
        rolls += 1
    return rolls


    

In [2]:

    

roll_me()


    

    
Out[2]:





3

In [3]:

    

[roll_me() for i in range(10)]


    

    
Out[3]:





[3, 2, 2, 3, 3, 2, 3, 2, 3, 2]

In [4]:

    

np.mean([roll_me() for i in range(1000)])


    

    
Out[4]:





2.7090000000000001

In [5]:

    

x = np.array([[n+1 for i in range(1000)] for n in range(1,200)])
jitter_x = x + np.random.rand(*x.shape)
trials = np.array([[roll_me(n+1) for i in range(1000)] for n in range(1,200)])
jitter_trials = trials + np.random.rand(*trials.shape)


    

In [6]:

    

import matplotlib.pyplot as plt
%matplotlib inline

plt.xkcd()

means = np.mean(trials, axis=1)
model = [np.exp(1) for n in range(1,200)]

plt.plot(means, label='simulation')
plt.plot(model, '--', label='e')
plt.legend(loc='lower right')
plt.xlabel("sides on die")
plt.ylabel("ave number rolls")
plt.yticks([0,1,2,3,4])
plt.ylim([0,4])
plt.text(20, 2, '1000 simulations per\nnumber of sides')
plt.show()


    

In [7]:

    

plt.xkcd()

plt.plot(jitter_x, jitter_trials, 'k.', alpha=1/100, label=None)
plt.plot(means, label="simulation means")
plt.plot(model, '--', label="e")
plt.legend()
plt.xlabel("sides on die")
plt.ylabel("number of rolls")
plt.ylim([0,14])
plt.xlim([0,200])
plt.show()