MAT411 Bayesian Data Analysis

Spring 2021: A New Hope

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 Random numbers

np.random.random()

0.9398514154371357

np.random.random(3)

array([0.08871139, 0.11951492, 0.60484217])

np.random.random([2,3])

array([[0.83377845, 0.113502  , 0.48395982],
       [0.96972018, 0.13614564, 0.01108457]])

np.random.random([2,3,2])

array([[[0.74973199, 0.19211608],
        [0.6692399 , 0.70883508],
        [0.34076713, 0.58215675]],

       [[0.53063272, 0.54821062],
        [0.91736751, 0.28088648],
        [0.15885051, 0.31746666]]])

plt.figure(figsize=(13,5))
plt.plot(np.random.random(100),'.')

[<matplotlib.lines.Line2D at 0x7fed71e8e7c0>]

plt.figure(figsize=(13,5))
bob=plt.hist(np.random.random(1000000),bins=100)

np.random.normal(0,1,3)

array([-0.51844007, -0.13443273, -0.50652096])

plt.plot(np.random.normal(0,1,1000),'.')

[<matplotlib.lines.Line2D at 0x7fed206e5eb0>]

plt.figure(figsize=(13,5))
bob=plt.hist(np.random.normal(0,1,1000000),bins=200)

sns.displot(np.random.normal(0,1,100000),kind="kde")

<seaborn.axisgrid.FacetGrid at 0x7fed60960760>

PDF: Probability Density Function

$$ f(x) = \frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} $$

def f(x,mu,sig):
    return 1/(sig*(2*np.pi)**0.5)*np.exp(-0.5*((x-mu)/sig)**2)

f(1,0,1)

0.24197072451914337

plt.figure(figsize=(7,6))

x=np.linspace(-3,3,9999)
x1 = np.linspace(-2,0,9999)

plt.plot(x,f(x,0,1),',')
plt.fill_between(x1,f(x1,0,1))

plt.xlim(-4,4)
plt.ylim(0,0.45)
plt.show()

import Pythonic_files as pf

def g(x):
    return f(x,0,1)
pf.Midpoint_int(-2,0,g,1000)

0.47724988604880725

Enter Scipy

x = np.linspace(-4,4,9999)

y = stats.norm(0,1).pdf(x)
plt.figure(figsize=(7,6))
plt.plot(x,y,',')

plt.show()

stats.norm(0,1).pdf(1)

0.24197072451914337

What about discrete distributions

stats.binom(n=4,p=0.5).pmf(4)

0.0625

trials =3

x= np.arange(0,trials+1)

y = stats.binom(n=trials,p=0.5).pmf(x)

plt.vlines(x,0,y,lw=10)

<matplotlib.collections.LineCollection at 0x7fed3865a910>

trials =20

x= np.arange(0,trials+1)

y = stats.binom(n=trials,p=0.3).pmf(x)

plt.vlines(x,0,y,lw=10)

<matplotlib.collections.LineCollection at 0x7fed41d6b6d0>

sns.displot(stats.norm(0,1).rvs(20000))

<seaborn.axisgrid.FacetGrid at 0x7fed81b58e50>

Flipping a coin

np.random.random()>0.5

True

if np.random.random()>0.5:
    toss = "H"
        
else:
    toss = "T"
print(toss)

T

flips=10
for f in range(0,flips):
    if np.random.random()>0.5:
        print('H')
    else:
        print("T")
    

T
H
T
T
H
H
H
T
T
H

flips=15
outcome=[]
for f in range(0,flips):
    if np.random.random()>0.5:
        outcome.append('H')
    else:
        outcome.append('T')
outcome

['H', 'T', 'T', 'T', 'T', 'H', 'T', 'H', 'T', 'H', 'H', 'T', 'H', 'T', 'H']

probability_H = outcome.count('H')/len(outcome)
print(probability_H)

0.4666666666666667

Can we plot this probability as the experiment runs

flips = 10
moving_prob = []
outcome = []
for f in range(0,flips):
    if np.random.random()>0.5:
        outcome.append('H')
    else:
        outcome.append('T')
        
    probability = outcome.count('H')/len(outcome)
    moving_prob.append(probability)

moving_prob

[0.0,
 0.0,
 0.3333333333333333,
 0.5,
 0.6,
 0.6666666666666666,
 0.7142857142857143,
 0.625,
 0.5555555555555556,
 0.6]

flips = 3000
moving_prob = []
outcome = []
for f in range(0,flips):
    if np.random.random()>0.5:
        outcome.append('H')
    else:
        outcome.append('T')
        
    probability = outcome.count('H')/len(outcome)
    moving_prob.append(probability)

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

plt.plot(moving_prob,',-')
plt.hlines(0.5,0,flips,color='r')
plt.ylim(-0.1,1.1)
plt.show()

Lets have a look at the resources needed for this experiment

flips = 30000
moving_prob = []
outcome = []
for f in range(0,flips):
    if np.random.random()>0.5:
        outcome.append('H')
    else:
        outcome.append('T')
        
    probability = outcome.count('H')/len(outcome)
    moving_prob.append(probability)

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

plt.plot(moving_prob,',-')
plt.hlines(0.5,0,flips,color='r')
plt.ylim(-0.1,1.1)
plt.show()

%%time

flips = 3000
moving_prob = []
outcome = []
for f in range(0,flips):
    if np.random.random()>0.5:
        outcome.append('H')
    else:
        outcome.append('T')
        
    probability = outcome.count('H')/len(outcome)
    moving_prob.append(probability)

CPU times: user 53.1 ms, sys: 3.29 ms, total: 56.4 ms
Wall time: 59.5 ms

np.array([1,2,3,4,5])/np.array([2,2,2,2,2])

array([0.5, 1. , 1.5, 2. , 2.5])

%%time
flips=300
results = np.random.random(flips)*2//1
results_sum = results.cumsum()
results_prob = results_sum/np.ones(flips).cumsum()

CPU times: user 197 µs, sys: 4.73 ms, total: 4.92 ms
Wall time: 5 ms

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

plt.plot(results_prob,'-')
plt.hlines(0.5,0,flips,color='r')
plt.ylim(-0.1,1.1)
plt.show()