Chapter 3: posterior distributions for the gaussian mean and dispersion

% Gaussian-priors-for-mean-and-sigma example

mureal=-2.1; sdreal=1;

load GaussDat.mat

mumin=-6.0; mumax=6.0;

sigmin=.1; sigmax=2.5;

clear D

muvec=mumin:.1:mumax; sigvec=sigmin:.1:sigmax;

pmu=ones(length(muvec),1)/length(muvec);

psig=1./sigvec; psig=psig/max(psig);

pi1=pmu(:)*psig(:)';

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

Ns=[0 10 20 50 100 500];

%Ns=[0 10 15 20 25 30];

[M V]=meshgrid(sigvec,muvec);

steps1=1; steps2=4;

figure(1); clf; hold on

mesh(M(1:steps2:end,3:steps1:end),V(1:steps2:end,3:steps1:end),pi1(1:steps2:end,3:steps1:end)), colormap(bone), brighten(-.8)

plot3(M(1:end,3),V(1:end,3),pi1(1:end,3),'k-','LineWidth',.8)

for n=1:steps2:length(muvec), plot3(M(n,[3 3+steps1]),V(n,[3 3+steps1]),pi1(n,[3 3+steps1]),'k-','LineWidth',1.2), end

stem3(M(1:steps2:end,3:steps1:end),V(1:steps2:end,3:steps1:end),pi1(1:steps2:end,3:steps1:end),'k.')

r=axis; axis([r(1:4) 0 .00275]), view(30,25)

pi1=log(pi1)-log(max(max(pi1)));

muvec=mumin:.01:mumax; sigvec=sigmin:.01:sigmax; sigind=find(sigvec==sdreal);

[M V]=meshgrid(sigvec,muvec);

pmu=ones(length(muvec),1)/length(muvec);

psig=1./sigvec; psig=psig/max(psig);

pi1=log(pmu(:)*psig(:)'); pi1=pi1-max(max(pi1));

fnow=zeros(length(muvec),length(sigvec));

for n=2:length(Ns),

Nnew=Ns(n-1)+1:Ns(n); Dnew=Gdraw(Nnew); %size(D)

fnew=zeros(length(muvec),length(sigvec),length(Dnew));

for nn=1:length(Dnew),

fnew(:,:,nn)=floglike(muvec,sigvec,Dnew(nn)); end

fnow=fnow+sum(fnew,3);

post1=fnow+pi1; post1=exp(post1-max(max(post1)));

Ncut=150; skips=4;

if n==2, xs=[]; cutoff=1; preal=sum(post1,2)/max(sum(post1,2));

while length(xs)<Ncut, cutoff=cutoff/2;

xs=find(preal>cutoff); end

ys=[]; cutoff=1; preal=sum(post1,1)/max(sum(post1,1));

while length(ys)<Ncut, cutoff=cutoff/2;

ys=find(sum(preal,1)>cutoff); end, end

figure(n); clf; hold on

mesh(M(xs(1:skips:end),ys(1:skips:end)),V(xs(1:skips:end),ys(1:skips:end)),post1(xs(1:skips:end),ys(1:skips:end))), colormap(bone)

stem3(M(xs(1:skips:end),ys(1:skips:end)),V(xs(1:skips:end),ys(1:skips:end)),post1(xs(1:skips:end),ys(1:skips:end)),'k.'); r=axis;

stem3(ones(1,length(xs(1:skips:end)))*sigvec(ys(1)),muvec(xs(1:skips:end)),sum(post1(xs(1:skips:end),:),2)/max(sum(post1(xs(1:skips:end),:),2)),'ko','filled','k')

stem3(sigvec(ys(1:skips:end)),ones(1,length(ys(1:skips:end)))*muvec(xs(end)),sum(post1(:,ys(1:skips:end)),1)/max(sum(post1(:,ys(1:skips:end)),1)),'ko','filled','k')

axis([r(1:4) 0 1]), view(30,25); end