cd PLExp1 ; S = PLEdata_by_sbj , save S ;

for k=1:length(S) subplot(3,5,k);s=S(k);plot(mean(s.dprime,2),'.-');axis([0 33 0 3]);set(gca,'xtick',[0:4:32]);grid on;ylabel('d''');title(sprintf('Sbj %d, Gr %d,  med RT = %d',s.sbj,s.group,round(median(s.data(:,8)))));end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   D-PRIME CURVES  (32 blocks)

%%%% By subjects
N=length(S);
for k=1:N subplot(3,5,k);s=S(k);md=mean(s.dprime,2);plot(1,md(1,:),'r.',2:9,md(2:9,:),'b.-',10:17,md(10:17,:),'r.-',18:25,md(18:25,:),'b.-',26:31,md(26:31,:),'r.-',32,md(32,:),'b.');axis([0 33 0 3]);set(gca,'xtick',[0:8:32]);grid on;title(sprintf('S %d, Gr %d, RT %d',s.sbj,s.group,round(mean(s.tmean_RT(:)))));end
%%%% Overall
d3=zeros(32,3,N);for k=1:N d3(:,:,k)=S(k).dprime;end;
md=mean(squeeze(mean(d3,3)),2);plot(1,md(1,:),'ro',2:9,md(2:9,:),'bs-',10:17,md(10:17,:),'ro-',18:25,md(18:25,:),'bs-',26:31,md(26:31,:),'ro-',32,md(32,:),'bs');axis([0 33 0 3]);set(gca,'xtick',[0:4:32]);grid on;ylabel('d''');xlabel('Block number');title([int2str(length(S)),' subjects']);
%%%% By difficulty levels
md=squeeze(mean(d3,3));plot(1,md(1,:),'ro',2:9,md(2:9,:),'bs-',10:17,md(10:17,:),'ro-',18:25,md(18:25,:),'bs-',26:31,md(26:31,:),'ro-',32,md(32,:),'bs');axis([0 33 0 3]);set(gca,'xtick',[0:4:32]);grid on;ylabel('d''');xlabel('Block number');title([int2str(length(S)),' subjects, 3 contrast levels']);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   Z-PROBABILITY CURVES  (32 blocks)

%%%% By subjects and "congruence", averaged over the three contrast levels
for k=1:N subplot(3,5,k);s=S(k);zP=squeeze(mean(s.zP_congr,2));plot(zP,'.-');axis([0 33 -1 4]);set(gca,'xtick',[0:8:32]);grid on;title(sprintf('S %d, Gr %d, RT %d',s.sbj,s.group,round(mean(s.tmean_RT(:)))));end
%%%% Same, but averaged across pairs of blocks to reduce variability
zP2=zeros([1+15+1,2,N]);idx=[3:2:31];for k=1:N zP=squeeze(mean(S(k).zP_congr,2));foo=(zP(idx-1,:)+zP(idx,:))/2;zP2(:,:,k)=[zP(1,:);foo;zP(31,:)];end;idx=[1,idx,32];
for k=1:N subplot(3,5,k);s=S(k);zP=zP2(:,:,k);plot(idx(1),zP(1,1),'r.',idx(1),zP(1,2),'g.',idx(2:5),zP(2:5,1),'r.-',idx(2:5),zP(2:5,2),'g.-',idx(6:9),zP(6:9,1),'r.-',idx(6:9),zP(6:9,2),'g.-',idx(10:13),zP(10:13,1),'r.-',idx(10:13),zP(10:13,2),'g.-',idx(14:16),zP(14:16,1),'r.-',idx(14:16),zP(14:16,2),'g.-',idx(17),zP(17,1),'r.',idx(17),zP(17,2),'g.');axis([0 33 -1 4]);set(gca,'xtick',[0:8:32]);grid on;title(sprintf('S %d, Gr %d, RT %d',s.sbj,s.group,round(mean(s.tmean_RT(:)))));end
%%%% All subjects together
zP=zeros([32,3,2]);for k=1:N zP=zP+S(k).zP_congr;end;zP=zP./N;
mz=squeeze(mean(zP,2));plot(1,mz(1,1),'ro',1,mz(1,2),'gs',2:9,mz(2:9,1),'ro-',2:9,mz(2:9,2),'gs-',10:17,mz(10:17,1),'ro-',10:17,mz(10:17,2),'gs-',18:25,mz(18:25,1),'ro-',18:25,mz(18:25,2),'gs-',26:31,mz(26:31,1),'ro-',26:31,mz(26:31,2),'gs-',32,mz(32,1),'ro',32,mz(32,2),'gs');axis([0 33 0 3]);legend('Incongr','Congr',4);set(gca,'xtick',[0:4:32]);grid on;ylabel('2 \Phi^{-1}(P_{correct})');xlabel('Block number');title([int2str(N),' subjects, 3 contrast levels']);
%%%% By congruence, all subjects together
subplot(2,1,1);mz=squeeze(zP(:,:,1));plot(1:32,mz(:,1),'bs-',1:32,mz(:,2),'ro-',1:32,mz(:,3),'g^-');axis([0 33 -1 3]);set(gca,'xtick',[0:4:32]);grid on;ylabel('2 \Phi^{-1}(P_{correct})');title('Incongruent stimuli');
subplot(2,1,2);mz=squeeze(zP(:,:,2));plot(1:32,mz(:,1),'bs-',1:32,mz(:,2),'ro-',1:32,mz(:,3),'g^-');axis([0 33 -1 3]);legend('C_p=.10','C_p=.15','C_p=.23',4);set(gca,'xtick',[0:4:32]);grid on;ylabel('2 \Phi^{-1}(P_{correct})');title('Congruent stimuli');xlabel('Block number');
%%%% By contrast levels, all subjects together
for k=1:3 subplot(3,1,k);mz=squeeze(zP(:,k,:));plot(1,mz(1,1),'ro',1,mz(1,2),'gs',2:9,mz(2:9,1),'ro-',2:9,mz(2:9,2),'gs-',10:17,mz(10:17,1),'ro-',10:17,mz(10:17,2),'gs-',18:25,mz(18:25,1),'ro-',18:25,mz(18:25,2),'gs-',26:31,mz(26:31,1),'ro-',26:31,mz(26:31,2),'gs-',32,mz(32,1),'ro',32,mz(32,2),'gs');axis([0 33 -1 3]);legend('Incongr','Congr',4);set(gca,'xtick',[0:4:32]);grid on;ylabel('2 \Phi^{-1}(P_{correct})');if (k==3) xlabel('Block number');end;title([int2str(N),' subjects, contrast ',int2str(k)]);end
%%%% By subjects and contrast levels
for f=0:3 figure;for n=1:4 s=S(n+4*f);for k=1:3 subplot(4,3,k+3*n-3);mz=squeeze(s.zP_congr(:,k,:));mz2=[mz(1,:);(mz(2:2:30,:)+mz(3:2:31,:))/2;mz(32,:)];plot(idx(1),mz2(1,1),'r.',idx(1),mz2(1,2),'g.',idx(2:5),mz2(2:5,1),'r.-',idx(2:5),mz2(2:5,2),'g.-',idx(6:9),mz2(6:9,1),'r.-',idx(6:9),mz2(6:9,2),'g.-',idx(10:13),mz2(10:13,1),'r.-',idx(10:13),mz2(10:13,2),'g.-',idx(14:16),mz2(14:16,1),'r.-',idx(14:16),mz2(14:16,2),'g.-',idx(17),mz2(17,1),'r.',idx(17),mz2(17,2),'g.');axis([0 33 -1 5]);set(gca,'xtick',[0:4:32]);grid on;title(sprintf('S %d, Gr %d, Contrast %d',s.sbj,s.group,k));end;end;end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   ROC PLOTS  by Context and Contrast

C = analyze_PLE_context(S) ; N=length(C);
for k=1:N zHit=norminv(C(k).P_right(:,:,2));zFA=norminv(C(k).P_right(:,:,1));subplot(3,4,k);plot(zFA(:,1),zHit(:,1),'ro-',zFA(:,2),zHit(:,2),'bs-');legend('C1','C2',4);axis([-2 +2 -2 +2]);hh=line([-2 +2]',[-2 +2]');set(hh,'Color','k');hh=line([-2 0]',[+2 0]');set(hh,'Color','k');title(sprintf('Sbj %d, Gr %d',S8(k).sbj,S8(k).group));end
zHit=zeros(3,2);zFA=zeros(3,2);for k=1:N zHit=zHit+norminv(C(k).P_right(:,:,2));zFA=zFA+norminv(C(k).P_right(:,:,1));end;zHit=zHit./N;zFA=zFA./N;
plot(zFA(:,1),zHit(:,1),'ro-',zFA(:,2),zHit(:,2),'bs-');legend('Context 1','Context 2',4);axis([-2 +2 -2 +2]);axis('square');hh=line([-2 +2]',[-2 +2]');set(hh,'Color','k');hh=line([-2 0]',[+2 0]');set(hh,'Color','k');title([int2str(length(C)),' subjects']);xlabel('z(FA)');ylabel('z(Hit)');grid on;
d=dprime(normcdf(zHit),normcdf(zFA))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  PROBABILITY and RT PLOTS by Context and Contrast

pR=zeros([3,2,2]);mRT=zeros([3,2,2]);for k=1:N pR=pR+C(k).P_right;mRT=mRT+C(k).tmean_RT;end;pR=pR./N;mRT=mRT./N;
for k=1:3 subplot(2,2,k);p=squeeze(pR(k,:,:));plot(1:2,p(1,:),'ro-',1:2,p(2,:),'bs-');legend('Context 1','Context 2',4);axis([0 3 0 1]);xlabel('Stimulus');ylabel('Prob(Resp=2)');title(['Difficulty level ' int2str(k)]);set(gca,'xtick',[1 2],'ytick',[0:0.1:1]);hh=refline(0,0.5);set(hh,'color','k','linestyle','-');end
for k=1:3 subplot(2,2,k);p=squeeze(mRT(k,:,:));plot(1:2,p(1,:),'ro-',1:2,p(2,:),'bs-');legend('Context 1','Context 2',4);axis([0 3 400 800]);xlabel('Stimulus');ylabel('mean RT');title(['Difficulty level ' int2str(k)]);set(gca,'xtick',[1 2]);end


for k=1:4 for d=1:3 subplot(4,3,3*k+d-3);foo=squeeze(p(k,d,:,:));plot(1:2,foo(2,:),'ro-',1:2,foo(1,:),'bs:');axis([0 3 0 1]);set(gca,'xtick',[1 2]);hh=refline(0,0.5);set(hh,'color','k','linestyle','-');end;end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   RESPONSE TIMES   (32 blocks)

N=length(S);for k=1:N subplot(3,5,k);s=S(k);mRT=mean(squeeze(mean(s.tmean_RT,3)),2);plot(1,mRT(1,:),'r.',2:9,mRT(2:9,:),'b.-',10:17,mRT(10:17,:),'r.-',18:25,mRT(18:25,:),'b.-',26:31,mRT(26:31,:),'r.-',32,mRT(32,:),'b.');axis([0 33 400 1100]);set(gca,'xtick',[0:4:32]);grid on;title(sprintf('Sbj %d, Gr %d',s.sbj,s.group));end
r = zeros(size(S(1).congr_RT));for k=1:N r=r+S(k).congr_RT;end;r=r./N;
mr=mean(squeeze(mean(r,3)),2);plot(1,mr(1),'ro',2:9,mr(2:9),'bs-',10:17,mr(10:17),'ro-',18:25,mr(18:25),'bs-',26:31,mr(26:31),'ro-',32,mr(32),'bs');axis([0 33 500 900]);set(gca,'xtick',[0:4:32]);grid on;ylabel('Trimmed mean RT');
%%%%  By difficulty levels
mr3=squeeze(mean(r,3));plot(mr3-repmat(mr,1,3),'.-');axis([0 33 -25 +25]);set(gca,'xtick',[0:4:32]);grid on;ylabel('Residual RT');legend('C_p=0.10','C_p=0.15','C_p=0.23');
%%%%  By congruence
mr2=squeeze(mean(r,2));plot(mr2,'.-');axis([0 33 500 900]);set(gca,'xtick',[0:4:32]);grid on;ylabel('Trimmed mean RT');xlabel('Block number');title([int2str(N) ' subjects']);foo=round(mean(mr2));legend(['Incongruent: mean RT=',int2str(foo(1))],['Congruent  : mean RT=',int2str(foo(2))]);
%%%%  By difficulty and congruence
foo = squeeze(mean(r,1));Cp=[0.10 0.15 0.23]';plot(Cp,foo(:,1),'ro-',Cp,foo(:,2),'gs-');axis([0.05 0.28 590 640]);xlabel('Contrast level');ylabel('Trimmed mean RT');set(gca,'xtick',Cp);legend('Incongruent','Congruent');
for k=1:N  s=S(k);foo=squeeze(mean(s.congr_RT,1));subplot(3,5,k);plot(Cp+0.001,foo(:,1),'r.-',Cp+0.001,foo(:,2),'g.-');set(gca,'xtick',Cp+0.001,'xticklabel','10|15|23');title(sprintf('Sbj %d, Gr %d',s.sbj,s.group));end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   D-PRIME and RT COMBINED

N=length(S);allD=zeros(32,N);allR=zeros(32,N);for k=1:N s=S(k);d=mean(s.dprime,2);r=mean(s.tmean_RT,2);allD(:,k)=d;allR(:,k)=r;subplot(4,4,k);plot(r,d,'.-');axis([400 1200 0 3]);cc=corrcoef([d,r]);cc=cc(2)^2;title(sprintf('Sbj %d, Gr %d,  R^2=.%02d',s.sbj,s.group,round(cc*100)));end
d=mean(allD,2);r=mean(allR,2);figure;plot(r(1),d(1),'ro',r(2:9),d(2:9),'bs-',r(10:17),d(10:17),'ro-',r(18:25),d(18:25),'bs-',r(26:31),d(26:31),'ro-',r(32),d(32),'bs');axis([500 900 0 3]);grid on;cc=corrcoef([d,r]);cc=cc(2)^2;title(sprintf('%d Ss,  R^2=.%02d',N,round(cc*100)));ylabel('d''');xlabel('Mean RT  (10% trimmed)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%   CURVE FITTING

xdata = [1:32 ; Inf,1:8,1:8,1:8,1:6,1]';
%%  Accuracy
N=length(S);d=zeros(32,N);for k=1:N d(:,k)=mean(S(k).dprime,2);end;md=mean(d,2);
x0=[2,1,0.1,-0.3,0.7];x=lsqcurvefit('PLE_exponent',x0,xdata,md),RMSE=RMSE=sqrt(mean((md-PLE_exponent(x)).^2))
F=PLE_exponent(x);plot(1,F(1),'k.',2:9,F(2:9),'k.-',10:17,F(10:17),'k.-',18:25,F(18:25),'k.-',26:31,F(26:31),'k.-',32,F(32),'k.');axis([0 33 0 3]);set(gca,'xtick',[0:4:32]);grid on;
hold on;plot(1,md(1),'ro',2:9,md(2:9),'bs',10:17,md(10:17),'ro',18:25,md(18:25),'bs',26:31,md(26:31),'ro',32,md(32),'bs');hold off;xlabel('Block number');ylabel('d''');title(sprintf('Exponential fit, RMSE=%.3f,  %d Ss',RMSE,N));
%%  Latency
r = zeros(size(S(1).tmean_RT));for k=1:N r=r+S(k).tmean_RT;end;r=r./N;mr=mean(squeeze(mean(r,3)),2);
y0=[500,800,0.2,+60,0.1];y=lsqcurvefit('PLE_exponent',y0,xdata,mr),RMSE=sqrt(mean((mr-PLE_exponent(y)).^2))
F=PLE_exponent(y);plot(1,F(1),'k.',2:9,F(2:9),'k.-',10:17,F(10:17),'k.-',18:25,F(18:25),'k.-',26:31,F(26:31),'k.-',32,F(32),'k.');axis([0 33 400 900]);set(gca,'xtick',[0:4:32]);grid on;
hold on;plot(1,mr(1),'ro',2:9,mr(2:9),'bs',10:17,mr(10:17),'ro',18:25,mr(18:25),'bs',26:31,mr(26:31),'ro',32,mr(32),'bs');hold off;xlabel('Block number');ylabel('Trimmed mean RT');title(sprintf('Exponential fit, RMSE=%.3f,  %d Ss',sqrt(mean((mr-F).^2)),N));
