% script to illustrate updating of beta posterior after 
% coin throws

% Chris Williams, December 2000

p = 0.7;  % p = P(heads) ; biased coin!

ah = 1; % alpha_heads = 1
at = 1; % alpha_tails = 1

iter = 50;

rand('state',0);  % set seed
R = rand(1,iter);
X = (R < p);   % puts out a 1 if heads, 0 otherwise
ax = 0:0.01:1; % for plotting [0,1] interval

for i = 1:iter
  for j = 1:i		% print out string so far -- neater than disp
   fprintf(1,'%d',X(1,j));
  end
  fprintf(1,'\n');
  heads = sum(X(1,1:i));
  tails = i - heads;
  f = ax.^(ah+heads) .* (1-ax).^(at + tails);
  hold off
  plot(ax, f/max(f));   % to get nice vertical [0,1] interval for plot
  hold on
  plot([p p], [0.01 0.99], 'r')
  pause(1)
end