%% Naive (brute-force) estimation of a Bernoulli variable

n = 10; %number of trials to simulate
q = .3; %true parameter value
o = rand(n,1)<q; %observation sequence
r = .01; %resolution of discrete approximation for distributions
d = zeros(floor(1+1/r),n+1); %prior/posterior distribution
p = (0:r:1)'; %parameter values associated with rows of d
d(:,1) = 1/size(d,1); %initialize prior (set to be uniform)

%run bayesian model
for i=1:n
    if o(i)==1
        d(:,i+1) = d(:,i).*p; %prior times likelihood
    else
        d(:,i+1) = d(:,i).*(1-p); %prior times likelihood
    end
    d(:,i+1) = d(:,i+1)/sum(d(:,i+1)); %normalize
end

%plot evolution of posterior
s = .7; %scaling of max of each distribution
t = 1; %spacing of trials to display
figure(1),clf,hold on
for i=1:t:n+1
    plot(p,i-1+d(:,i)*t*s/max(d(:,i))) %prior for trial i (posterior from i-1)
    plot(0:1,[i-1,i-1],'k') %baseline for plot
    if t==1 && i<n+1
        plot(o(i),i,'r.','MarkerSize',20) %current observation
    end
end
axis([-.001 1 0 n+t-.01])

%% Bernoulli estimation with Beta conjugate prior

%first a look at the beta distribution
figure(2),clf
plot(p,betapdf(p,2,5))

%setup for simulation
n = 10; %number of trials to simulate
q = .3; %true parameter value
o = rand(n,1)<q; %observation sequence
a = zeros(n,1); b = zeros(n,1); %parameters of prior, across trials
a(1) = 1; b(1) = 1; %parameters of initial prior

%run bayesian model
for i=1:n
    if o(i)==1 %increment 1s count (alpha)
        a(i+1) = a(i) + 1;
        b(i+1) = b(i);
    else %increment 0s count (beta)
        a(i+1) = a(i);
        b(i+1) = b(i) + 1;
    end
end

%even faster (gives same result as above)
a(2:n+1) = a(1) + cumsum(o==1);
b(2:n+1) = b(1) + cumsum(o==0);

%plot evolution of posterior
s = .7; %scaling of max of each distribution
t = 1; %spacing of trials to display
figure(3),clf,hold on
for i=1:t:n+1
    plot(p,i-1+betapdf(p,a(i),b(i))*t*s/max(betapdf(p,a(i),b(i)))) %prior for trial i (posterior from i-1)
    plot(0:1,[i-1,i-1],'k') %baseline for plot
    if t==1 && i<n+1
        plot(o(i),i,'r.','MarkerSize',20) %current observation
    end
end
axis([0 1 0 n+t-.01])
