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.