function sheet05_01
% Comparing Bayesian estimation versus maximum-likelihood estimation.

% generate some Bernoulli distributed data
N = 300;
p = 0.3;
X = rand(1, N) < p;

% estimate p using the ml estimator
ml = zeros(1, N);
for I = 1:N
ml(I) = ml_estimate(X(1:I));
end

% estimate p using a "non-informative"
[B11A, B11B] = beta_estimate(1, 1, X);
[B7030A, B7030B] = beta_estimate(70, 30, X);
[B300700A, B300700B] = beta_estimate(300, 700, X);

B11A

plot(1:N, ml, 'g');
hold on;
plot_beta_estimates(B11A, B11B, 'r');
plot_beta_estimates(B7030A, B7030B, 'k');
plot_beta_estimates(B300700A, B300700B, 'b');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% only make changes blow

% 1. Compute the maximum-likelihood estimate for Bernoulli random
% variables from the samples contained in X
function p = ml_estimate(X)
% ...

% Compute the continuous updates of the beta distribution prior
% parameters for X.
function [A, B] = beta_estimate(a, b, X)
N = length(X);
A = zeros(1, N);
B = zeros(1, N);
for I = 1:N
[a, b] = beta_update(a, b, X(I));
A(I) = a;
B(I) = b;
end

% 2. Update the a and b parameters of the beta distribution prior
% for a single X.
function [a, b] = beta_update(a, b, X)
% ...

% 3. Compute the mean and variance of the beta distribution for
% the given vector of parameters
function [M, V] = beta_stats(A, B)
% ...

% Plot mean and two times standard deviation.
function plot_beta_estimates(A, B, mark)
N = length(A);
[M, V] = beta_stats(A, B);
plot(1:N, M, [mark, '-'], ...
1:N, M + 2*sqrt(V), [mark, ':'], ...
1:N, M - 2*sqrt(V), [mark, ':']);