1 function sheet05_01
   2 % Comparing Bayesian estimation versus maximum-likelihood estimation.
   3 
   4 % generate some Bernoulli distributed data
   5 N = 300;
   6 p = 0.3;
   7 X = rand(1, N) < p;
   8 
   9 % estimate p using the ml estimator
  10 ml = zeros(1, N);
  11 for I = 1:N
  12   ml(I) = ml_estimate(X(1:I));
  13 end
  14 
  15 % estimate p using a "non-informative"
  16 [B11A, B11B] = beta_estimate(1, 1, X);
  17 [B7030A, B7030B] = beta_estimate(70, 30, X);
  18 [B300700A, B300700B] = beta_estimate(300, 700, X);
  19 
  20 B11A
  21 
  22 plot(1:N, ml, 'g');
  23 hold on;
  24 plot_beta_estimates(B11A, B11B, 'r');
  25 plot_beta_estimates(B7030A, B7030B, 'k');
  26 plot_beta_estimates(B300700A, B300700B, 'b');
  27 hold off;
  28 
  29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  30 % only make changes blow
  31 
  32 % 1. Compute the maximum-likelihood estimate for Bernoulli random
  33 % variables from the samples contained in X
  34 function p = ml_estimate(X)
  35 % ...
  36 
  37 % Compute the continuous updates of the beta distribution prior
  38 % parameters for X.
  39 function [A, B] = beta_estimate(a, b, X)
  40 N = length(X);
  41 A = zeros(1, N);
  42 B = zeros(1, N);
  43 for I = 1:N
  44   [a, b] = beta_update(a, b, X(I));
  45   A(I) = a;
  46   B(I) = b;
  47 end
  48 
  49 % 2. Update the a and b parameters of the beta distribution prior
  50 % for a single X.
  51 function [a, b] = beta_update(a, b, X)
  52 % ...
  53 
  54 % 3. Compute the mean and variance of the beta distribution for
  55 % the given vector of parameters
  56 function [M, V] = beta_stats(A, B)
  57 % ...
  58 
  59 % Plot mean and two times standard deviation.
  60 function plot_beta_estimates(A, B, mark)
  61 N = length(A);
  62 [M, V] = beta_stats(A, B);
  63 plot(1:N, M, [mark, '-'], ...
  64      1:N, M + 2*sqrt(V), [mark, ':'], ...
  65      1:N, M - 2*sqrt(V), [mark, ':']);