Chapter 2: marginalization over nuisance parameters

We begin by computing the 2D posterior over the Gaussian parameters:

mu=linspace(-10,10,41); mu=mu(2:end)-diff(mu(1:2))/2; %equally-spaced mu values

sig=linspace(0,15,31); sig=sig(2:end)-diff(sig(1:2))/2; %sigmas

D=nrand(0,5,11);

p=nan(length(mu),length(sig)); %initialize 2D-probsurface

The dataset (D) stays constant while we compute the posterior probabilities of mu-sigma-pairs:

for snow=1:length(sig), %for each value of sigma

for mnow=1:length(mu), %for each value of mu

Lnow=1; %reset p for ith mu and jth sig

for dnow=1:length(D), %for each datum

%gaussian likelihood

Lnow=Lnow/(sqrt(2*pi)*sig(snow))*exp(-.5*(D(dnow)-mu(mnow))^2/sig(snow)^2); end

p(mnow,snow)=Lnow; end, end %put p at ith mu and jth sig into 2D matrix

figure; meshc(sig,mu,p/max(max(p))); colormap bone

view(-145,45); r=axis; hold on

Note the innermost for-loop computes the likelihood, and the other two loops step through mu-sigma-pairs. Next marginalize over sigma axis, yielding the desired 1D distribution over mu:

pmu=sum(p,2); pmu=pmu/max(pmu); %marginalize over SIGMAs

for mnow=1:length(mu), %plotting

plot3(sig(end)*[1 1],mu(mnow)*[1 1],[0 pmu(mnow)],'k-');

plot3(sig(end),mu(mnow),pmu(mnow),'ko','MarkerFaceColor','k'); end