% Implementation of Coin Toss example in 
% Smith, Uncertainty Quantification, Section 4.8

 q = 0.5;        % probability of heads

 N = input('number of trials:  ');

 Y = (rand(N,1) <= q);    % Y = 1 if heads, 0 if tails

 for k=1:N, if Y(k)==1, fprintf('H'), else, fprintf('T'); end, end
 fprintf('\n\n')
 
 N0 = sum(Y==0);   % number of tails 
 N1 = sum(Y==1);   % number of heads

 clear pi0 pi
% define a prior density (not a good one)
 mu = 0.3; 
 sig = 0.1;
 pi0 = @(q) (1/(sig*sqrt(2*pi))).*exp(-(((q-mu).^2)/(2*sig^2)));
 pi0int = quad(@(q) pi0(q),0,1);

% define the conditional density
 pi = @(q) (q.^(N1)).*((1-q).^(N0));
 jointint = quad(@(q) pi(q).*pi0(q),0,1)

 Q = linspace(0,1,500);
 figure(1), clf
 plot(Q,pi0(Q)/pi0int,'r-','linewidth',2)
 hold on
 plot(Q,pi(Q).*pi0(Q)/jointint,'b-','linewidth',2)
 title(sprintf('%d Heads, %d Tails', N1, N0),'interpreter','latex','fontsize',20)
 xlabel('$q$','fontsize',18,'interpreter','latex')
 leg = legend('prior','posterior','location','northeast','orientation','vertical');
 set(leg,'fontsize',18,'interpreter','latex')
  
