Chapter 6: frequentist algorithm versus 2-model comparison

clear; close all %#ok<*UNRCH>

plotFLAG=1;

load rawslopestatloop.mat

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

N=6; Nrep=251; Nrot=25; Nmod=2;

sig=15; Ns=401; slist=sig; %noise level

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

rotlist=atan(linspace(0,1,Nrot));

%keep noise samples identical across simulations for clearer comparison

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

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

margL=nan(Nrot,Nmod,Nrep); Elist=nan(Nrot,Nmod,Nrep); Plist=nan(Nrot,Nrep);

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)

%or it is (deviation from unity line)

tnow=dnow(:,1)'; dnow=dnow(:,2)'+eorig(irep,:); deltStr{1}='dnow.^2';

deltStr{2}='(dnow-tnow).^2'; for nn=1:size(Elist,2),

eval(['D2=' deltStr{nn} ';'])

eval(['likeFN' num2str(nn-1) '=-(N/2)*(mean(D2,2)./slist.^2)-(N+1)*log(slist);']); end

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

%%% Stat'l significance %%%

try

tbl=table(tnow',dnow'-mean(dnow)); %,'Variable Names',{'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)=ttt.Coefficients.pValue(2); %p-value for null zero-slope model

catch load stat1; end

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

%%% Model likelihoods %%%

margL(irot,2,irep)=likeFN1; %marginal likelihood of the

%unity-slope model, over unknown SD

margL(irot,1,irep)=likeFN0;

%marginal likelihood of the zero-slope model, over %unknown SD

Elist(irot,1:2,irep)=evid(margL(irot,:,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 -900 900]);

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([-150 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);

axis([-.05 1.05 -.02 1.02]); end