Chapter 6: frequentist algorithm versus 3-model comparison

clear %#ok<*UNRCH>

plotFLAG=1; verifysigFLAG=0;

COL=[.55*[1 1 1]; .3 .4 .85; .8 .5 .4];

N=6; Nrep=251; Nrot=25; Nmod=3; sig=15;

Ns=401; slist=linspace(0,101,Ns+1)'; slist=slist(2:end); %noise level

blist=linspace(0,1,201); %slopes used to approximate the underlying continuous slope dimension

rotlist=atan(linspace(0,1,Nrot)); %slopes to use

diffs=[diff(slist(1:2)) diff(blist(1:2))];

thetas=linspace(-100,100,N);

%keep noise samples identical across simulations for clearer comparison

eorig=nrand(0,sig,[Nrep N]);

dorig=[thetas(:) zeros(N,1)];

%load sim params

%load(['rawslopestatloop3.mat']);

if verifysigFLAG, sest=zeros(Nrot,Nmod); end

margL=nan(Nrot,Nmod,Nrep); Elist=nan(Nrot,Nmod,Nrep); Emax=nan(Nrot,Nmod); Plist=nan(Nrot,Nrep); p1blistmat2=(1+ones(Ns,1)*blist.^2); Ipri=[.4636 .3217];

for irot=1:Nrot, rotnow=rotlist(irot);

for irep=1:Nrep,

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

%%% DATA SIM %%%

dnow=rot(dorig,rotnow);

%delta-squared is either just (deviation from 0)^2, or it is (deviation from unity-line)^2

tnow=dnow(:,1)'; dnow=dnow(:,2)'+eorig(irep,:);

deltStr{1}='dnow.^2'; deltStr{2}='(dnow-betas(n)*tnow).^2';

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

%%% param likelihood by condition %%%

likeFN0=nan(Ns,sum(blist<.5)); likeFN1=nan(Ns,sum(blist>.5));

for n=1:sum(blist<=.5),

tmp=-(N/2)*(mean((dnow-blist(n)*tnow).^2,2)./slist.^2)-(N+1)*log(slist); %L term

likeFN0(:,n)=tmp-log(Ipri(1)*p1blistmat2(:,n)); end %add prior term (normalized by Ipri)

for nn=0:sum(blist>.5)-1,

tmp=-(N/2)*(mean((dnow-blist(nn+n+1)*tnow).^2,2)./slist.^2)-(N+1)*log(slist); %L term

likeFN1(:,nn+1)=tmp-log(Ipri(2)*p1blistmat2(:,n)); end %add prior term (normalized by Ipri)

Lall0=-(N/2)*(mean((dnow-0.0*tnow).^2,2)./slist.^2)-(N+1)*log(slist);

Lall1=-(N/2)*(mean((dnow-1.0*tnow).^2,2)./slist.^2)-(N+1)*log(slist);

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

%%% param likelihood by condition %%%

%Lall0=logsum(likeFN0,2,diffs(2));

%Lall1=logsum(likeFN1,2,diffs(2));

Lall5=-(N/2)*(mean((dnow-0.5*tnow).^2,2)./slist.^2)-(N+1)*log(slist);

if verifysigFLAG,

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

%%% verify global sigma %%%

sest(irot,2)=logpeakval([Lall1 slist]);

sest(irot,1)=logpeakval([Lall0 slist]); end

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

%%% Stat'l significance %%%

try tbl=table(tnow',dnow'-mean(dnow),'VariableNames',{'x','y'}); %fit the null model (fitlm requires table datatype)

ttt=fitlm(tbl,'linear'); %returns ttt struct, which contains 'Coefficients' table for the fit of the null model

Plist(irot,irep,1)=ttt.Coefficients.pValue(2); %p-vallue for null (zero-slope) model

tbl=table(tnow',(dnow'-mean(dnow))-.5*tnow','VariableNames',{'x','y'}); ttt=fitlm(tbl,'linear'); %unity slope model

Plist(irot,irep,2)=ttt.Coefficients.pValue(2); %p-vallue for null (unity-slope) model

catch load stat3; end

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

%%% Model likelihoods %%%

margL(irot,3,irep)=logsum(Lall5,1,diffs(1)); %marginal likelihood of the 0.5-slope model

margL(irot,2,irep)=logsum(Lall1,1,diffs(1)); %marginal likelihood of the unity-slope model, over unknown SD and slope

margL(irot,1,irep)=logsum(Lall0,1,diffs(1)); %marginal likelihood of the zero-slope model, over unknown SD and slope

Elist(irot,1:Nmod,irep)=evid(margL(irot,1:Nmod,irep)); end, end

if plotFLAG,

Hlist=Plist<.05;

BLU=[.3 .4 .65]; Glist=[.1 .5 .75];

figure(1); clf; hold on

plot([.2 .8],[0 0],'k--','LineWidth',1.75); axis([-.05 1.05 -150 150]);

for nn=1:Nmod, emed=median(Elist(:,nn,:),3); imed3=find(emed>3);

plot(tan(rotlist(imed3))',emed(imed3),'o','MarkerFaceColor',BLU,'MarkerEdgeColor', BLU,'MarkerSize',13);

plot(tan(rotlist)'*ones(1,Nrep),squeeze(Elist(:,nn,:)),'k.');

plot(tan(rotlist)',emed,'ko:','MarkerFaceColor',Glist(nn)*[1 1 1],'MarkerSize',9); end

xlabel('data slope','FontName','Arial','FontSize',15);

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

figure(2); clf; hold on; plot([0 1],[0 0],'k--')

for nn=1:Nmod,

plot(tan(rotlist)'*ones(1,Nrep),squeeze(margL(:,nn,:)),'.','Color',Glist(nn)*[1 1 1])

plot(tan(rotlist),median(squeeze(margL(:,nn,:)),2),'ko:','MarkerFaceColor',Glist(nn)*[1 1 1],'MarkerSize',13); end

xlabel('data slope','FontName','Arial','FontSize',15); xlim([-.05 1.05])

ylabel('mod likelihood','FontName','Arial','FontSize',15); ylim([-30 0])

figure(3); clf; hold on

plot(tan(rotlist)'*ones(1,Nrep),Plist,'.','Color',Glist(1)*[1 1 1])

plot(tan(rotlist),mean(Plist,2),'ko:','MarkerFaceColor',Glist(1)*[1 1 1],'MarkerSize',13);

plot([0 1],.05*[1 1],'k--','LineWidth',1.75)

xlabel('data slope','FontName','Arial','FontSize',15);

ylabel('p-value','FontName','Arial','FontSize',15);

xlim([-.05 1.05]); ylim([-.02 1.02]), end