% Bayesian estimation example

% Input parameters

N = 100;
H = 55;
p = 0.005:0.01:0.995; % Center of intervals in p 
dp = 0.01;

% Panel 1: Likelihood function
mu = N*p;
sigma = sqrt(N*p.*(1-p));
L = exp(-(H-mu).^2./(2*sigma.^2))./(sigma*sqrt(2*pi));
subplot(221)
plot(p,L,'k-')
xlabel('p')
ylabel('P(55H45T|p)')

% Panels 2-4: posterior probability for various priors

% Uniform prior (equal chance of each p in [0,1])

np = length(p);
prior = ones(1,np);  
subplot(222)
post = prior.*L;
post = post/sum(post*dp);  % Bayes formula (divide by dp to turn into PDF)
plot(p,prior,'b',p,post,'r')
xlabel('p')
ylabel('PDF')
legend('prior','posterior','Location','Northwest')
title('Uniform prior for p')

subplot(223)
plot(p,cumsum(prior*dp),'b',p,cumsum(post*dp),'r',[0 1],...
     [0.025 0.025],'k--',[0 1],[0.975 0.975],'k--')
ylim([0 1])
xlim([0.35 0.65])
xlabel('p')
ylabel('CDF')
legend('prior','posterior','Location','Northwest')
title('Uniform prior for p')

subplot(224)
sigprior = 0.05;
prior = exp(-(p-0.5).^2./(2*sigprior.^2))./(sigprior*sqrt(2*pi));
post = prior.*L;
post = post/sum(post*dp);  % Bayes formula (divide by dp to turn into PDF)
plot(p,cumsum(prior*dp),'b',p,cumsum(post*dp),'r',[0 1],...
     [0.025 0.025],'k--',[0 1],[0.975 0.975],'k--')
ylim([0 1])
xlim([0.35 0.65])
xlabel('p')
ylabel('CDF')
legend('prior','posterior','Location','Northwest')
title('Gaussian prior for p, mu=0.5, std 0.05')


