Scneario: We toss a coin and get 8 heads and 2 tails. Our objective is to determine probability of getting head in the next toss. Use bayesian approach to estimate parameter i.e. probability of getting head ( $\theta$ )
Assuming $\theta$ has a beta distribution with $\alpha = \beta = 2.0$, expected value for $\theta$ is given as
$E[\theta] = \frac{n_H + \alpha}{\alpha + \beta + N}$
$E[\theta] = \frac{8 + 2}{2 + 2 + 10} = \frac{10}{14} = 0.71$
%matplotlib inline
import random
import scipy.stats
import scipy.special
import matplotlib.pyplot as plt
def integral(theta):
""" This is likelihood X prior after removing any constant terms."""
return theta**9*(1-theta)**3
n = 100 # Number of points to sample
samples = [random.random()] # Random starting point
for i in range(n):
#Grab last sample point
theta = samples[-1]
#Create new sample by randomly selecting a point from a normal distribution
#of mean = 0 and sd = 0.1. if the new sample is outside of the
#domain then ignore it and use existing sample
newTheta = theta + random.normalvariate(0, 0.1)
if newTheta < 0 or newTheta > 1:
newTheta = theta
#If the probability of new sample as compared to last sample is less than uniform distribution
#then ignore it.
acceptanceRatio = integral(newTheta)/integral(theta)
if acceptanceRatio > random.random(): # accept only if going uphill
samples.append(newTheta)
else:
samples.append(theta)
print "Estimate: ", sum(samples)/n
# Plot how sample theta varies with each iteration
ylab = [i for i in xrange(len(samples))]
pylab.plot(samples, ylab)
pylab.title('Random Walk Visualization')
pylab.xlabel('Theta Value')
pylab.ylabel('Time')
pylab.show()
Estimate: 0.649221876089