Chapter 6: numerical assessment of the evidence
Because the (parameter) likelihood term is identical for each of the three models, we can compute model likelihoods by first computing the sampled likelihood over the rate parameter, and then integrating with respect to the prior over the rate parameter defined within each of the three models separately. This yields:
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 3-model comparison %%%
clear
%%% prelim
Nmus=1085; mumin=-12; mumax=12; mus=linspace(mumin,mumax,Nmus)'; Om=ones(Nmus,1);
Nsigs=204; sigmax=20; sigs=linspace(0,sigmax,Nsigs+1); sigs=sigs(2:end); sign2=sigs.^-2;
diffs=[diff(mus(1:2)) diff(sigs(1:2))];
%%% parameter prior for three models
im1=find(and(mus>0,mus<1)); is1=findnearestN(sigs,5,1);
im2=findnearestN(mus,1,1); is2=find(and(sigs>1,sigs<7));
im3=findnearestN(mus,0,1); is3=find(and(sigs>3,sigs<10));
%keep noise constant across simulated datasets
N=24; err0=randn(N,1); err0=(err0-mean(err0))/std(err0);
%%% 2D prior and likelihood over mu / sigma
primat=-Om*log(sigs); %figure(1); clf; meshc(sigs,mus,primat)
lrange2=log(7)-log(1); lrange3=log(10)-log(3);
like=@(xbar,x2bar) -N*(log(Om*sqrt(2*pi))+log(sigs)+.5*(((xbar-mus).^2)*sign2+Om*(x2bar-xbar^2)*sign2)+.5*log(2*pi)); %Gaussian likelihood
mulist=linspace(-.1,1.1,57); siglist=linspace(1,12,35);
ML=zeros(length(mulist),length(siglist),3);
EV=zeros(length(mulist),length(siglist),3);
for im=1:length(mulist), munow=mulist(im);
for is=1:length(siglist), signow=siglist(is);
dat=err0*signow+munow;
likenow=like(mean(dat),mean(dat.^2));
ML(im,is,1)=logsum(likenow(im1,is1),1,diffs(1));
ML(im,is,2)=logsum(likenow(im2,is2)+primat(im2,is2)/lrange2,2,diffs(2));
ML(im,is,3)=logsum(likenow(im3,is3)+primat(im3,is3)/lrange3,2,diffs(2));
EV(im,is,:)=evid(str8n(ML(im,is,:))','natural'); end, end
figure(2); clf
subplot(2,1,1); plot(mus,exp(logsum(likenow,2,diffs(2))),'k.')
subplot(2,1,2); plot(sigs,exp(logsum(likenow,1,diffs(1))),'k.')
figure(3); clf; meshc(siglist,mulist,ML(:,:,1))
figure(4); clf; meshc(siglist,mulist,ML(:,:,2))
figure(5); clf; meshc(siglist,mulist,ML(:,:,3))
figure(6); clf; meshc(siglist,mulist,EV(:,:,1))
figure(7); clf; hold on; imu=round(length(mulist));
plot(siglist,squeeze(EV(imu,:,1)),'ko-','MarkerFaceColor','k','MarkerSize',10)
plot(siglist,squeeze(EV(imu,:,2)),'ko-','MarkerFaceColor',.7*[1 1 1],'MarkerSize',10)
plot(siglist,squeeze(EV(imu,:,3)),'ko-','MarkerFaceColor',[.3 .4 .8],'MarkerSize',10)
ylabel('evidence (dB)','FontName','Arial','FontSize',16)
xlabel(['data sd (\sigma)'],'FontName','Arial','FontSize',16)
axis([0 13 -150 100])
It is not necessary, however, that the parameter likelihood be the same in all models. For example, we might have reason to propose a fourth model that predicts an exponential error function, whose mean is in the range: . This yields:
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 4-model comparison %%%
%%% add tau param and tau prior
Ntaus=501; taumin=0; taumax=5; taus=linspace(taumin,taumax,Ntaus+1); taus=taus(2:end); ltaus=log(taus);
diffs=[diff(mus(1:2)) diff(sigs(1:2)) diff(taus(1:2))];
pritaus=-ltaus;
datasig=2; datamean=mulist(imu);
dG=1; while or(-min(dG)>.25*datamean,sum(dG<0)~=2), dG=randn([N 1])*datasig+datamean; end
dE=0; while any(dE<=0), dE=erand(datamean,[N 1]); dE=dE-mean(dE)+datamean; end
Nrat=21; ratlist=linspace(0,1,Nrat)';
%%% parameter likelihood for exponential model
elike=@(dset) -dset*(1./taus)-ones(N,1)*ltaus;
ML=zeros(Nrat,4);
EV=zeros(Nrat,4);
for irat=1:Nrat, ratnow=ratlist(irat);
dnow=dG*(1-ratnow)+dE*ratnow;
likenow=like(mean(dnow),mean(dnow.^2));
ML(irat,1)=logsum(likenow(im1,is1),1,diffs(1));
ML(irat,2)=logsum(likenow(im2,is2)+primat(im2,is2)/lrange2,2,diffs(2));
ML(irat,3)=logsum(likenow(im3,is3)+primat(im3,is3)/lrange3,2,diffs(2));
likenow=elike(dnow); likenow((dnow*ones(1,Ntaus))<=0)=-inf; likenow=sum(likenow,1);
ML(irat,4)=logsum(likenow+pritaus,1,diffs(3));
EV(irat,:)=evid(ML(irat,:),ones(1,4),'natural'); end
figure(8); clf; hold on
plot(ratlist,EV(:,1),'ko-','MarkerFaceColor','k','MarkerSize',10)
plot(ratlist,EV(:,2),'ko-','MarkerFaceColor',.7*[1 1 1],'MarkerSize',10)
plot(ratlist,EV(:,3),'ko-','MarkerFaceColor',[.3 .4 .8],'MarkerSize',10)
plot(ratlist,EV(:,4),'ko-','MarkerFaceColor',[.7 .4 .3],'MarkerSize',10)
ylabel('evidence (dB)','FontName','Arial','FontSize',16)
xlabel(['combination ratio'],'FontName','Arial','FontSize',16)
axis([0 1 1.05*[min(min(EV)) max(max(EV))]])
Notice that this model comparison yields infinite evidence against the exponential model when any observed datum is negative. Also, you should recognize that the evidence values for the two favorable models would not change in any appreciable way if the remaining two models were removed from the analysis.