MAT 411 Bayesian Data Analysis

The Prior strikes back

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy import stats


%matplotlib inline
%config InlineBackend.figure_format='retina'

Lets look at this problems except, whats the probability of flipping 4 heads out of 10

def Num_heads(flips):
    results = np.random.random(flips)
    
    results=results*2//1
    results = results.sum()
    return results

Num_heads(10)

6.0

Samples = []
experiment = 2000000
for x in range(0,experiment):
    if Num_heads(10)==4:
        Samples.append(1)
    else:
        Samples.append(0)
Samples = np.array(Samples)


results = Samples.cumsum()
results = results/np.ones(experiment).cumsum()

plt.plot(results,',')
plt.show()

results[-1]

0.205282

We know this is Binomial Experiment and thus not the exact answer

Binomial Probability = $P(x) = \left(_x^n\right)p^xq^{n-x} = \frac{n!}{(n-x)!x!}p^xq^{n-x}$

def nCk(n,k):
    return np.math.factorial(n)/(np.math.factorial(n-k)*np.math.factorial(k))

nCk(10,3)

120.0

Answer = nCk(10,4)*0.5**4*0.5**(10-4)
Answer

0.205078125

print("The sampling result is :",results[-1])
print('The exact answer is :', Answer)

The sampling result is : 0.20473
The exact answer is : 0.205078125

A Difference between Frequentist and Bayesian

flips = 30
results = []
for x in range(0,10000):
    results.append(Num_heads(flips)/flips)

plt.hist(results,bins=30)
plt.show()

sns.displot(results)

<seaborn.axisgrid.FacetGrid at 0x7fbdf07e2dc0>

flips = 30
results = []
experiments=300
for x in range(0,experiments):
    results.append(Num_heads(flips)/flips)

print('This is the mean',np.mean(results))
sns.displot(results)
plt.vlines(np.mean(results),0,.3*experiments,color='g')
plt.show()

This is the mean 0.4976666666666666

Confidence interval = $ \mu \pm z^*\frac{\sigma}{\sqrt{n}}$

std = np.std(results)
Zstar = 1.96 # 95% confidence interval

flips = 30
results = []
experiments=20000
for x in range(0,experiments):
    results.append(Num_heads(flips)/flips)

print('This is the mean',np.mean(results))
print('standard deviation:', np.std(results))

mu = np.mean(results)
std = np.std(results)
Zstar = 1.96

sns.displot(results)
plt.vlines(mu,0,.3*experiments,color='g')

plt.vlines(mu-Zstar*std/np.sqrt(experiments), 0 ,.3*experiments, color='r')
plt.vlines(mu+Zstar*std/np.sqrt(experiments), 0 ,.3*experiments, color='r')

plt.show()

This is the mean 0.4995866666666667
standard deviation: 0.0922890762766645

Monty Hall Problem

Game show:

The contenstant chooses a door from 3, then the host picks another door and opens it( it will contain a goat)

So either the door that is not picked or the contestants door contains the big prize and the other a donkey.

Question What should the contestant do to improve his chances, either stick with the door or choose the last door?

Choices = ['G',"D","C"]

np.random.choice(Choices)

'D'

Choices.remove('G')

Choices

['D', 'C']

"G" in Choices

False

So lets run the game, using this type of code and find out what is the best stratergy

def game_stay():
    Choices = ['G',"D","C"]
    Play = np.random.choice(Choices)
    if Play == "C":
        result = 1
        
    else:
        result = 0
        
    return result

Winning_prob= []
experiments=100000
for x in range(0,experiments):
    Choices = ['G',"D","C"]
    Play = np.random.choice(Choices)
    if Play == "C":
        result = 1
        
    else:
        result = 0
    Winning_prob.append(result)

    
End_prob = np.cumsum(np.array(Winning_prob))
results = End_prob/np.ones(experiments).cumsum()

plt.figure(figsize=(13,6))
plt.plot(results)

plt.hlines(0.5,0,experiments, color='r')

plt.ylim(0,1)

plt.show()

results[-1]

0.33341

Diabetic Robber

We know a person robbed a store, because the left some blood, we now know they are diabetic.

In America the prevalence of diabetes is 14.6% in men and 9.1% in females. Given we know the blood from the robber, whats the probability of them being male.

def Robbery():
    no_lead=True
    count=0
    while no_lead:
        sex = np.random.choice(["M","F"])
        Diabetic_Chance = np.random.random()
        if sex =="M" and Diabetic_Chance<=0.146:
            break
        elif sex=="F" and Diabetic_Chance<=0.091:
            break
        count+=1
    return sex  

Robbery()

'F'

Diabetic_Robbers =[]
moving_prob=[]

for x in range(0,2000):
    Diabetic_Robbers.append(Robbery())
    
    probability =Diabetic_Robbers.count('M')/len(Diabetic_Robbers)
    moving_prob.append(probability)
plt.plot(moving_prob)
plt.hlines(0.5,0,len(moving_prob),color='r')
plt.ylim(0,1)

(0.0, 1.0)

for x in range(0,2000):
    Diabetic_Robbers.append(Robbery())
    
    probability =Diabetic_Robbers.count('M')/len(Diabetic_Robbers)
    moving_prob.append(probability)
    

0.61

def Robbery(num):
    R_sample=[]
    while len(R_sample)<num:
        sex = np.random.choice(["M","F"])
        Diabetic_Chance = np.random.random()
        if sex =="M" and Diabetic_Chance<=0.146:
            R_sample.append(1)
        elif sex=="F" and Diabetic_Chance<=0.091:
            R_sample.append(0)
        
    return R_sample  

Robbery(20)

[1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1]

end_prob=[]
for c in range(0,100):
    Sample_output = np.array(Robbery(1000))
    results=Sample_output.cumsum()
    results=results/np.ones(len(results)).cumsum()
    plt.plot(results, '-',c='b',alpha=0.3)

end_prob=[]
for c in range(0,100):
    Sample_output = np.array(Robbery(1000))
    results=Sample_output.cumsum()
    results=results/np.ones(len(results)).cumsum()
    end_prob.append(results[-1])

sns.displot(end_prob,kind='kde')

<seaborn.axisgrid.FacetGrid at 0x7fbe32340d90>

Bayes Theorem

$$ P(A|B) = ? $$

What is the blue?

So blue divided by all of original green (B) ....

$$ P(A|B) = \frac{P(A\cap B)}{P(B)} $$

So our analytical solution to our Robbery problem is

$$ P(\text{Male Given Diabetic}) = \frac{P(\text{Male and Diabetic})}{P(\text{Diabetic}) }$$$$ P(\text{Male Given Diabetic}) = \frac{0.5\times 0.146}{0.5\times 0.146 + 0.5\times0.091} $$

sns.displot(end_prob,kind='kde')

<seaborn.axisgrid.FacetGrid at 0x7fbe41d48880>

0.5*.146/(0.5*.146+0.5*0.091)

0.6160337552742616

np.mean(end_prob)

0.61786