Chapter 6: multisource slope comparisons

%%% SLM with global dispersion and no intercepts

%%% [zero-one slope comparison]

clear %#ok<*UNRCH>

reset=1;

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

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

slist=linspace(0,50,Ns+1)'; slist=slist(2:end); slistn2=slist.^(-2); lnslist=log(slist);%noise level

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

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

Nrot=22;

rotlist=linspace(-pi/6,pi/4+pi/6,Nrot); %rotations in radians

%keep noise samples identical across simulations for clearer comparison

brel=urand(-.1,.1,[S,1]);

dorig=[thetas(:) thetas(:)]; eorig=nrand(0,sig,[S N]);

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

margL=zeros(Nrot,2); sest=zeros(Nrot,2);

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

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

%%% param likelihood by condition %%%

for s=1:S,

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

%%% DATA SIM using rotation ratio %%%

drot=rot(dorig,rotnow-pi/4+atan(brel(s))); tnow=drot(:,1)'; dnow=drot(:,2)'+eorig(s,:);

likeFN0=-(N/2)*slistn2*mean(dnow.^2,2)-(N+1)*lnslist; %output is Ns x 1

likeFN1=-(N/2)*slistn2*mean((dnow-tnow).^2,2)-(N+1)*lnslist; %output is Ns x 1

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

%%% verify global sigma %%%

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

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

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

%%% Model likelihoods %%%

margL(irot,1)=margL(irot,1)+logsum(likeFN0,1,diffs(1));

margL(irot,2)=margL(irot,2)+logsum(likeFN1,1,diffs(1)); end, end

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

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(30); clf; hold on;

plot(tan(rotlist),margL(:,1),'ko:','MarkerFaceColor','k');

plot(tan(rotlist),margL(:,2),'ko:','MarkerFaceColor',COL(1,:))

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

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