BayesRuleBinomial

clear all;
% set random number seeds.
s=999; rand('seed',s); randn('seed',s);

% Make vector of possible bias values
inc = 0.001; % set resolution of distributions.
bmin = inc;
bmax = 1-inc;
b = bmin:inc:bmax;	% range of bias values p(head).
a = 1-b;          		% range of p(tails) values.

% make prior either uniform or binomial.
% set switch for uniform prior.
uprior=0;
if uprior     % Make uniform prior. 
    ht=1./(bmax-bmin);
    prior = ones(size(b)).*ht;
else  % Make prior=binomial
    C = nchoosek(4,2); % not really needed here.
    p = C.* b.^2 .* a.^2;
    prior = p/max(p);
end


% plot prior
figure(1);set(gca,'FontSize',20);
plot(b,prior,'k-','Linewidth',2); set(gca,'Linewidth',2);
xlabel('Coin bias, \theta'); ylabel('Prior, p(\theta)');
set(gca,'YLim',[0 1.1],'FontName','Ariel'); grid on;

ni=8; if uprior==0 ni=3; end % make graphs for figure 4.3.
NN = 2^ni; % max number of coin flips in loop below is NN.
flips = rand(NN,1); % get NN random numbers between zero and one.
x0 = flips<0.6; % get number of heads for a coin with bias 0.6.
% Cheat a bit as want to show distribution for specific outcomes.
x0(1)=1; x0(2)=0; x0(3)=1; x0(4)=1; fprintf('\n');

% find likelihood and posterior for different numbers N of coin flips.
for i=1:ni % step up number of coin flips in powers of 2.
    N=2^i; % get number of coin flips.
    x = x0(1:N); % get data for first N coin flips.
    k = sum(x); % get number of heads.
    % C = nchoosek(N,k); % binmial coefficient, not needed here.
    C = 1; % set binomial coefficient to one.
    nh = k; % nh = number of heads.
    nt = N-k; % nt = number of tails.
    % likelihood function = probability of nh heads and nt tails.
    lik = C.* b.^nh .* a.^nt; 
    lik = lik/max(lik); % scale likelihood to have max value of one
    % plot likelihood function.
    figure(2);set(gca,'FontSize',20);
    plot(b,lik,'k-','Linewidth',2); set(gca,'Linewidth',2);
    xlabel('Coin bias, \theta'); ylabel('Likelihood, p({\bfx}|\theta)');
    set(gca,'YLim',[0 1.1],'FontName','Ariel'); grid on;
    p = lik.*prior;           % find posterior = likelihood * prior
    maxp = max(p);             % find max value of p
    ind = find(maxp==p);     % find index of max value of p
    p = p/maxp;                   % make max value of p be one.
    best = b(ind);                 % find estimated bias, 
    % same as (number of heads)/(number of coin flips)
    % plot posterior
    figure(3); set(gca,'FontSize',20); plot(b,p,'k','Linewidth',2); 
    set(gca,'YLim',[0 1.1],'FontName','Ariel');
    xlabel('Coin bias, \theta'); 
    ylabel('Posterior, p({\theta|\bfx})'); set(gca,'Linewidth',2);
    ss=sprintf('Heads = %d/%d',k,N);
    text(0.1,0.9,ss,'FontSize',20); grid on;
    % print value of bias associated with max value of p
    fprintf('Number of heads = %d. Estimated bias = %.3f\n',k ,best);
end