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