Chapter 6: multisource slope comparisons

%%% SLM with source-specific intercept and global dispersion

%%% [zero-one slope comparison]

clear %#ok<*UNRCH>

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

N=6; S=5; sig=10;

Na=201; alist=linspace(-3*sig,3*sig,Na); %data slope

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

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

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

rotlist=[-.3:.1:1.3]; Nrot=length(rotlist); %rotations in radians

%keep noise samples identical across simulations for clearer comparison

dorig=[thetas(:) thetas(:)]; eorig=nrand(0,sig,[S N]); alphas=urand(-2*sig,2*sig,[S 1]);

for s=1:S, eorig(s,:)=eorig(s,:)-mean(eorig(s,:))+alphas(s); end

Ov=ones(Na,1); Oh=ones(1,Ns); ON=ones(1,N);

modlike=zeros(Nrot,2);

sest0=zeros(Nrot,S); sest1=zeros(Nrot,S);

aest0=zeros(Nrot,S); aest1=zeros(Nrot,S);

Lall0=zeros(S,Ns,Na); Lall1=zeros(S,Ns,Na);

slistmat=slist*ones(1,Na); slistn2=slist.^(-2);

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

margL0=zeros(Ns,1); margL1=zeros(Ns,1);

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

%%% DATA SIM using rotation ratio %%%

drot=rot(dorig,rotnow-pi/4); tnow=drot(:,1)';

likeFN0=@(dnow) -(N/2)*slistn2*mean((Ov*(dnow)-alist'*ON).^2,2)'-(N+1)*log(slistmat); %output is Ns x Na

likeFN1=@(dnow) -(N/2)*slistn2*mean((Ov*(dnow-tnow)-alist'*ON).^2,2)'-(N+1)*log(slistmat); %output is Ns x Na

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

%%% param likelihood by condition %%%

for s=1:S, dnow=drot(:,2)'+eorig(s,:);

tmpL0=likeFN0(dnow);

tmpL1=likeFN1(dnow);

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

%%% verify global sigma %%%

margL0=margL0+logsum(tmpL0,2,diffs(2));

margL1=margL1+logsum(tmpL1,2,diffs(2));

sest0(irot,s)=logpeakval([logsum(tmpL0,2,diffs(2)) slist]);

aest0(irot,s)=logpeakval([logsum(tmpL0,1,diffs(1)) alist']);

sest1(irot,s)=logpeakval([logsum(tmpL1,2,diffs(2)) slist]);

aest1(irot,s)=logpeakval([logsum(tmpL1,1,diffs(1)) alist']); end

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

%%% Model likelihoods %%%

modlike(irot,1)=logsum(margL0,1,diffs(1));

modlike(irot,2)=logsum(margL1,1,diffs(1)); end

Elist=10*(modlike(:,2)-modlike(:,1));

figure(1); clf; hold on

plot(tan(rotlist),Elist,'ko-','MarkerFaceColor','k','MarkerSize',11);

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

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

figure(2); clf; hold on; plot(tan(rotlist),modlike(:,1),'ko:','MarkerFaceColor','k'); plot(tan(rotlist),modlike(:,2),'ko:','MarkerFaceColor',COL(1,:))

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

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