Chapter 3: posterior distributions for the gaussian mean

mureal=21; sd=7.5;

Gdraw=nrand(mureal,sd,[500 1]); Gdraw(1:5)=[20 15 20 22 18];

mumin=5; mumax=25;

clear D

mu0=mumin-5:.001:mumax+5;

mu=mumin:.001:mumax;

ld=(length(mu0)-length(mu))/2;

theta=(mu-mumin)/(mumax-mumin);

pbi=theta.^10.*(1-theta).^20+theta.^20.*(1-theta).^10; pbi=[zeros(ld,1); pbi(:)/sum(pbi); zeros(ld,1)];

pmi=theta.^5.*(1-theta).^5; pmi=[zeros(ld,1); pmi(:)/sum(pmi); zeros(ld,1)];

pu=ones(size(mu))/length(mu); pu=[zeros(ld,1); pu(:)/sum(pu); zeros(ld,1)];

pu0=ones(size(mu0))/length(mu0); pu0=pu0(:)/sum(pu0);

floglike=@(mus,D) (-(D(:)*ones(1,length(mus))-ones(length(D),1)*mus).^2/sd^2);

Ns=[0:5 10:5:20 50 100 500];

figure; subplot(4,3,1); hold on

pis=[pu0 pu pmi pbi]/max(pmi);

plot(mu0,pis(:,1),'k-','LineWidth',1.05)

plot(mu0,pis(:,2),'k:','LineWidth',1.1)

plot(mu0,2*pis(:,3),'k-.','LineWidth',2.1)

plot(mu0,2*pis(:,4),'k--','LineWidth',1.1)

%axis([min(mu0) max(mu0) 0 1])

for n=2:length(Ns),

N=Ns(n-1)+1:Ns(n); D=Gdraw(N);

fnow=sum(floglike(mu0,D),1);

fnow=exp(fnow-max(fnow));

fnow=fnow(:)*ones(1,size(pis,2));

post=fnow.*pis;

pnorm=max(post);

post=post./(ones(size(post,1),1)*pnorm);

%if n<5, post=post./max(post(:,4)); else post=post./max(max(post)); end

subplot(4,3,n), hold on

plot(mu0,post(:,1),'k-','LineWidth',1.05)

plot(mu0,post(:,2),'k:','LineWidth',1.1)

plot(mu0,post(:,3),'k-.','LineWidth',2.1)

plot(mu0,post(:,4),'k--','LineWidth',1.1)

%meandef(post,theta)

%sqrt(vardef(post,theta))

axis([min(mu0) max(mu0) 0 1]); end

xlabel('mu'); ylabel('p(mu|data)');