Chapter 6: multisource dispersion

%%% Biased Gaussian

%%% [models compare unbiased models that differ in whether dispersion is source specific or globsal]

clear %#ok<*UNRCH>

N=10; S=5;

COL={.4*[1 1 1],.8*[1 1 1]};

sigS=linspace(.2,10,S); sigG=mean(sigS)*ones(S,1);

Ns=501; slist=linspace(0,20,Ns+1)'; slist=slist(2:end); sn2=slist.^(-2);

slnrange=log(max(slist))-log(min(slist));

err0=nrand([S,N]);

for s=1:S, err0(s,:)=(err0(s,:)-mean(err0(s,:))); end

diffs=diff(slist(1:2));

ratlist=[0:.05:1]; Nrat=length(ratlist);

ML1=zeros(Nrat,1); ML2=zeros(Nrat,1);

for irat=1:Nrat, rnow=ratlist(irat);

Lall1=zeros(1); Lall2=zeros(Ns,1);

for s=1:S,

%%%%%%%%R%%%%%%%

%%% DATA SIM %%%

dnow=err0(s,:)*(sigG(s)*(1-rnow)+sigS(s)*(rnow));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% param likelihood by condition %%%

likeFN1=-(N/2)*sn2*mean((dnow).^2,2)'-(N+1)*log(slist); %output is Ns x 1

likeFN2=-(N/2)*sn2*mean((dnow).^2,2)'-(N )*log(slist); %output is Ns x 1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% marginalize source-specific sig %%%

Lall1=Lall1+logsum(likeFN1,1,diffs)-log(slnrange);

Lall2=Lall2+likeFN2; end

ML1(irat)=Lall1; ML2(irat)=logsum(Lall2-log(slist)-log(slnrange),1,diffs); end

%%% compute evidence

EV=10*(ML2-ML1);

%%% plotting

figure(4); clf; hold on

plot(ratlist,ML1,'ko','MarkerFaceColor',COL{1},'MarkerSize',10)

plot(ratlist,ML2,'ko','MarkerFaceColor',COL{2},'MarkerSize',10)

plot([ratlist(1)-.1 ratlist(end)+.1],[0 0],'k:')

xlabel(['ratio of xSource to global \sigma'],'FontSize',15);

ylabel('model likelihood','FontName','Arial','FontSize',15);

r=axis; axis([ratlist(1)-.1 ratlist(end)+.1 r(3) max([r(4) 10])])

figure(5); clf; hold on

plot(ratlist([1 end]),[0 0],'k--','LineWidth',2);

plot(ratlist(EV<0),EV(EV<0),'ko','MarkerFaceColor',COL{1},'MarkerSize',12,'LineWidth',1);

plot(ratlist(EV>0),EV(EV>0),'ko','MarkerFaceColor',COL{2},'MarkerSize',12,'LineWidth',1);

r=axis; axis([ratlist(1)-.1 ratlist(end)+.1 r(3:4)]);

xlabel(['ratio of xSource to global \sigma'],'FontSize',15);

ylabel('evidence (db)','FontName','Arial','FontSize',15);