#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')


# # Monte Carlo and de Mere Die Games

# What is the probability of rolling at least one 6 over 4 die? 
# 
# $(1 - \frac{5}{6}^4) \approx 0.517$

# In[11]:


a = np.random.randint(1, 7, 4000)


# In[12]:


a


# In[13]:


allprobs = []
trials = 1000
experiments = 300
for i in range(experiments):
    a = np.random.randint(1, 7, 4 * trials)
    allprobs.append(len([x for x in a.reshape(-1, 4) if 6 in x]) / trials)


# In[14]:


sum(allprobs) / len(allprobs)


# In[17]:


plt.hist(allprobs, 20)


# For students: What is probability of getting at least one pair of 6 in 24 rolls of pairs of die?
# 
# $(1 - \frac{35}{36}^{24}) \approx 0.491$

# In[164]:


[e for e in np.random.randint(1, 7, 48000).reshape((-1, 24, 2))]


# In[174]:


allprobs = []
trials = 1000
experiments = 300
for i in range(experiments):
    count = 0
    for t in [e for e in np.random.randint(1, 7, 48 * trials).reshape((-1, 24, 2))]:
        found = False
        for p in t:
            if p[0] == 6 and p[1] == 6:
                found = True
        if found:
            count += 1    
    allprobs.append(count / trials)


# In[175]:


sum(allprobs) / len(allprobs)


# In[176]:


plt.hist(allprobs, 20)


# # Random walks and Gambler's Ruin
# 
# Start with $20, flip a coin each turn, heads win 1, tails lose 1, for 200 timesteps. What is expected value?

# In[150]:


trials = 500
time = 500
init = 20
overall = []
for t in range(trials):
    initial = [init]
    for i in range(time):
        if np.random.random() > 0.517:
            initial.append(initial[-1] + 1)
        else:
            initial.append(initial[-1] - 1)
        if initial[-1] == 0:
            break
    if len(initial) < time + 1:
        initial += [0] * (time + 1 - len(initial))
    overall.append(initial)


# In[151]:


b = np.linspace(0, time, time + 1)
for o in overall:
    plt.plot(b, o)


# In[152]:


ends = []
for o in overall:
    ends.append(o[-1])


# In[153]:


plt.hist(ends, 20)


# In[154]:


sum(ends) / len(ends)


# In[155]:


aves = []
n = np.array(overall)
for i in range(time + 1): 
    t = n[:, i]
    aves.append(sum(t) / trials)


# In[156]:


plt.plot(aves)


# In[ ]: