In [1]:

    

import numpy as np
def cake_me(N=30):
    puffs = 0
    while N > 0:
        N -= np.random.randint(1, N+1) # Note that randint uses the interval [low, high)
        puffs += 1
    return puffs


    

In [2]:

    

cake_me()


    

    
Out[2]:





3

In [3]:

    

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


    

    
Out[3]:





[5, 3, 3, 3, 3, 4, 2, 5, 5, 2]

In [4]:

    

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


    

    
Out[4]:





4.0199999999999996

In [5]:

    

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


    

In [6]:

    

import matplotlib.pyplot as plt
%matplotlib inline
#plt.style.use('fivethirtyeight')
plt.xkcd()

means = np.mean(trials, axis=1)
model = [np.log(2*(n+1)) for n in range(100)]

plt.plot(means, label='simulation')
plt.plot(model, '--', label='log(2(N-1))')
plt.xlabel("number of candles")
plt.ylabel("ave number of puffs")
# plt.yticks([1,2,3,4,5])
plt.ylim([0,6])
plt.xlim([0,100])
plt.text(40, 3, '1000 simulations per\nnumber of candles')
plt.legend(loc='lower right')
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="log(2(N-1))")
plt.legend()
plt.xlabel("number of candles")
plt.ylabel("number of puffs")
plt.ylim([0,14])
plt.xlim([0,100])
plt.show()