MAT411 Bayesian Analsysis:

Week 7

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


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

lets look at unfair coins

flips = 10
heads_prob = 0.3 # np.random.random()

Heads_num = stats.binom(flips, heads_prob).rvs()
Tails_num = flips - Heads_num

print(Heads_num,Tails_num)

3 7

num_space=100
prob_space = np.linspace(1/num_space,num_space/(num_space+1),num_space)

unifom_space= np.ones(len(prob_space))/num_space

triangle_space = np.minimum(prob_space,1-prob_space)
triangle_space = triangle_space/triangle_space.sum()

prob_like = pf.nCk(flips,Heads_num)*prob_space**(Heads_num)*(1-prob_space)**(Tails_num)
prob_like=prob_like/prob_like.sum()

prob_prior = triangle_space

prob_norm = np.sum(prob_space*prob_like)
prob_norm = prob_norm/prob_norm.sum()


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

$$ P(D) = \sum_i P(D|\theta_i)*P(\theta_i) $$

plt.plot(prob_space,prob_like)
plt.plot(prob_space,prob_prior)
plt.plot(prob_space,prob_post)

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

flips = 6
heads_prob = 0.3 # np.random.random()

Heads_num = stats.binom(flips, heads_prob).rvs()
Tails_num = flips - Heads_num


num_space=100000
prob_space = np.linspace(1/num_space,num_space/(num_space+1),num_space)

unifom_space= np.ones(len(prob_space))/num_space

triangle_space = np.minimum(prob_space,1-prob_space)
triangle_space = triangle_space/triangle_space.sum()

prob_like = pf.nCk(flips,Heads_num)*prob_space**(Heads_num)*(1-prob_space)**(Tails_num)
prob_like=prob_like/prob_like.sum()

prob_prior = triangle_space

prob_norm = np.sum(prob_space*prob_like)
prob_norm = prob_norm/prob_norm.sum()


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_prior, '--', label ='prior')
plt.plot(prob_space,prob_like, label ='likelihood')
plt.plot(prob_space,prob_post, label = 'posterior')


plt.title('Unfair coin toss with '+str(Heads_num)+' Heads, '+str(Tails_num)+' Tails')

plt.legend()
plt.show()

Lets look at a square wave prior

flips = 60
heads_prob = 0.3 # np.random.random()

Heads_num = stats.binom(flips, heads_prob).rvs()
Tails_num = flips - Heads_num


num_space=1000
prob_space = np.linspace(1/num_space,num_space/(num_space+1),num_space)

unifom_space= np.ones(len(prob_space))/num_space

triangle_space = np.minimum(prob_space,1-prob_space)
triangle_space = triangle_space/triangle_space.sum()

wave = 3

def square_wave(t):
    if int(2*wave*t)%2==0:
        return 1
    else:
        return 0.1

square_prob = [square_wave(t) for t in prob_space]
square_prob= np.array(square_prob)
square_prob = square_prob/square_prob.sum()


prob_like = pf.nCk(flips,Heads_num)*prob_space**(Heads_num)*(1-prob_space)**(Tails_num)
prob_like=prob_like/prob_like.sum()

prob_prior = square_prob

prob_norm = np.sum(prob_space*prob_like)
prob_norm = prob_norm/prob_norm.sum()


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_prior, '--', label ='prior')
plt.plot(prob_space,prob_like, label ='likelihood')
plt.plot(prob_space,prob_post, label = 'posterior')


plt.title('Unfair coin toss with '+str(Heads_num)+' Heads, '+str(Tails_num)+' Tails')

plt.legend()
plt.show()

plt.plot(prob_space, prob_post)

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

square_prob = [square_wave(t) for t in prob_space]

plt.plot(square_prob)

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

Lets look at a Beta function

$$ B(x,y) = \int_0^1 t^{x-1}(1-t)^{y-1} dt $$

OR

$$ B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} $$

where $\Gamma(x) = \int_{u=0}^\infty e^{-u}u^{x-1} du $

what do they look like

a = np.linspace(1,6,6)
b = np.linspace(1,6,6)


x = np.linspace(0,1,100)

_, ax = plt.subplots(len(a),len(b),sharex=True, sharey=True, figsize = (15,15))

for i in range(0,len(a)):
    for j in range(0,len(b)):
        aa = a[i]
        bb = b[j]
        
        y = stats.beta.pdf(x,aa,bb,0,1)
        ax[i,j].plot(x,y)

flips = 20
heads_prob = 0.3 # np.random.random()

Heads_num = stats.binom(flips, heads_prob).rvs()
Tails_num = flips - Heads_num


num_space=1000
prob_space = np.linspace(1/num_space,num_space/(num_space+1),num_space)

unifom_space= np.ones(len(prob_space))/num_space

triangle_space = np.minimum(prob_space,1-prob_space)
triangle_space = triangle_space/triangle_space.sum()

wave = 3

def square_wave(t):
    if int(2*wave*t)%2==0:
        return 1
    else:
        return 0.1

square_prob = [square_wave(t) for t in prob_space]
square_prob= np.array(square_prob)
square_prob = square_prob/square_prob.sum()

a, b = 10,2

prob_beta = stats.beta.pdf(prob_space,a,b,0,1)
prob_beta = prob_beta/prob_beta.sum()


prob_like = pf.nCk(flips,Heads_num)*prob_space**(Heads_num)*(1-prob_space)**(Tails_num)
prob_like=prob_like/prob_like.sum()

prob_prior = prob_beta

prob_norm = np.sum(prob_space*prob_like)
prob_norm = prob_norm/prob_norm.sum()


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_prior, '--', label ='prior')
plt.plot(prob_space,prob_like, label ='likelihood')
plt.plot(prob_space,prob_post, label = 'posterior')


plt.title('Unfair coin toss with '+str(Heads_num)+' Heads, '+str(Tails_num)+' Tails')

plt.legend()
plt.show()

lets take a look at different prior of beta

flips = 10
heads_prob = 0.3 # np.random.random()

Heads_num = stats.binom(flips, heads_prob).rvs()
Tails_num = flips - Heads_num


num_space=100
prob_space = np.linspace(1/num_space,num_space/(num_space+1),num_space)


a, b = 10,2

prob_beta = stats.beta.pdf(prob_space,a,b,0,1)
prob_beta = prob_beta/prob_beta.sum()


prob_like = pf.nCk(flips,Heads_num)*prob_space**(Heads_num)*(1-prob_space)**(Tails_num)
prob_like=prob_like/prob_like.sum()

prob_prior = prob_beta

prob_norm = np.sum(prob_space*prob_like)
prob_norm = prob_norm/prob_norm.sum()


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_like, c='r',label ='likelihood')
plt.plot(prob_space,prob_prior, '--',c='b', label ='prior')

plt.plot(prob_space,prob_post,c='b', label = 'posterior')


plt.title('Unfair coin toss with '+str(Heads_num)+' Heads, '+str(Tails_num)+' Tails')


a, b = 2,2

prob_beta = stats.beta.pdf(prob_space,a,b,0,1)
prob_beta = prob_beta/prob_beta.sum()



prob_prior = prob_beta


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.plot(prob_space,prob_prior, '-.',c='g', label ='prior')

plt.plot(prob_space,prob_post,c='g', label = 'posterior')




a, b = 0.5,0.5

prob_beta = stats.beta.pdf(prob_space,a,b,0,1)
prob_beta = prob_beta/prob_beta.sum()



prob_prior = prob_beta


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.plot(prob_space,prob_prior, '-.',c='y', label ='prior')

plt.plot(prob_space,prob_post,c='y', label = 'posterior')





plt.title('Unfair coin toss with '+str(Heads_num)+' Heads, '+str(Tails_num)+' Tails')


plt.legend()
plt.show()


plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_post)
plt.vlines(prob_space[prob_post.argmax()],0, prob_post.max(), color='r')

<matplotlib.collections.LineCollection at 0x7fa2818c3b50>

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_post)
plt.vlines(prob_space[prob_post.argmax()],0, prob_post.max(), color='r')

<matplotlib.collections.LineCollection at 0x7fa24ad5cfa0>

Lets find mode, median, mean

prob_post.argmax()

27

prob_space[prob_post.argmax()]

0.2772997299729973

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_post)
plt.vlines(prob_space[prob_post.argmax()],0, prob_post.max(), color='r')

plt.title('mode')
plt.show()

prob_median = prob_post.cumsum()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_post)
plt.vlines(prob_space[(prob_median*2//1%2).argmax()], 0, prob_post.max(),color='r')
plt.title('median')
plt.show()

prob_mean = prob_space*prob_post
prob_mean.sum()

0.22728548267998466

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_post)
plt.vlines(prob_mean.sum(), 0, prob_post.max(),color='r')
plt.title('mean')
plt.show()

flips = 20
heads_prob = 0.3 # np.random.random()

Heads_num = stats.binom(flips, heads_prob).rvs()
Tails_num = flips - Heads_num


num_space=100
prob_space = np.linspace(1/num_space,num_space/(num_space+1),num_space)

unifom_space= np.ones(len(prob_space))/num_space

triangle_space = np.minimum(prob_space,1-prob_space)
triangle_space = triangle_space/triangle_space.sum()

wave = 2

def square_wave(t):
    if int(2*wave*t)%2==0:
        return 1
    else:
        return 0.1

square_prob = [square_wave(t) for t in prob_space]
square_prob= np.array(square_prob)
square_prob = square_prob/square_prob.sum()

a, b = 10,2

prob_beta = stats.beta.pdf(prob_space,a,b,0,1)
prob_beta = prob_beta/prob_beta.sum()


prob_like = pf.nCk(flips,Heads_num)*prob_space**(Heads_num)*(1-prob_space)**(Tails_num)
prob_like=prob_like/prob_like.sum()

prob_prior = square_prob

prob_norm = np.sum(prob_space*prob_like)
prob_norm = prob_norm/prob_norm.sum()


prob_post = prob_prior*prob_like/prob_norm
prob_post = prob_post/prob_post.sum()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_prior, '--', label ='prior')
plt.plot(prob_space,prob_like, label ='likelihood')
plt.plot(prob_space,prob_post, label = 'posterior')

plt.vlines(prob_space[prob_post.argmax()],0, prob_post.max(), color='r', label = 'mode')
prob_median = prob_post.cumsum()
plt.vlines(prob_space[(prob_median*2//1%2).argmax()], 0, prob_post.max(),color='r', label ='median')
prob_mean = prob_space*prob_post
plt.vlines(prob_mean.sum(), 0, prob_post.max(),color='r', label= 'mean')


plt.title('Unfair coin toss with '+str(Heads_num)+' Heads, '+str(Tails_num)+' Tails')

plt.legend()
plt.show()

plt.figure(figsize=(13,6))
plt.plot(prob_space,prob_post)

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

looking ahead

x,y =pf.elipse_noise(10000,27,(20,10),(100,200))
plt.figure(figsize=(13,6))

plt.scatter(x,y, alpha=0.051)

<matplotlib.collections.PathCollection at 0x7fa2316aedf0>

len(x)

10000