%%%%%%%%%%%  Fit the three learning curves simultaneously with a common time constant
%%%%   April 03, 2003


 load HUMAN_ZCONGR ; edpr=mean(HUMAN_ZCONGR)'.*2;
 blk1=[1,10:17,26:31]';blk2=[2:9,18:25,32]';

%% Three time variables: Clock time T (in blocks), time since last switch Ts, and time in each context Tc
 T=[1:32]';Ts=[Inf,1:8,1:8,1:8,1:6,1]';Tc=[1,1:8,2:9,9:16,10:15,17]';


%%%% Null hypothesis: the context has no effect whatsoever
 par0=lsqcurvefit('PLE_exponent',[2 0.5 10],[T Ts],edpr), par0=lsqcurvefit('PLE_exponent',par0,[T Ts],edpr)
par0 = [1.7165  0.5161  8.3497]
  F_asympt=1.7165, F_init=D*(1-g)=0.8306, g=0.5161, tau=8.3497

 mdpr0=PLE_exponent(par0);resid0=edpr-mdpr0;mrse0=std(resid0,1),cc=corrcoef(mdpr0,edpr);cc(2)^2
mrse0 = 0.1413 ; R2 = 0.7425

%old|  par0=lsqcurvefit('PLE_2exp',[2 1 .1],T,edpr), par0=lsqcurvefit('PLE_2exp',par0,T,edpr)
%old| 
%old| par0:  F_asympt=1.7184, F_init=0.8323, alpha=0.1187, tau=1/alpha=8.4246
%old| 
%old|  mdpr0=PLE_2exp(par0,T);resid0=edpr-mdpr0;mrse0=std(resid0,1),cc=corrcoef(mdpr0,edpr);cc(2)^2
%old| mrse0 = 0.1413 ; R2 = 0.7425


%%%% Switch-cost hypothesis: general curve + switch costs
 par1=lsqcurvefit('PLE_exponent',[2 0.5 10 .2 2],[T Ts],edpr),par1=lsqcurvefit('PLE_exponent',par1,[T Ts],edpr)
par1 = [1.8415  0.4749  10.3764  0.1799  1.1393]
  F_asympt=1.8415, F_init=D*(1-g)=0.9670, g=0.4749, tau=10.3764
  s=0.1799, s*D=0.3313, tau_s = 1.1393

 mdpr1=PLE_exponent(par1);resid1=edpr-mdpr1;mrse1=std(resid1,1),cc=corrcoef(mdpr1,edpr);cc(2)^2
mrse1 = 0.0825 ; R2 = 0.9123

%old|  par1=lsqcurvefit('PLE_exponent',[2 1 .1 -.3 .5],[T Ts],edpr),par1=lsqcurvefit('PLE_exponent',par1,[T Ts],edpr)
%old| 
%old| par1:  F_asympt=1.8448, F_init=0.9696, alpha=0.0952, tau=1/alpha=10.5042
%old|        swcost=-0.3320, beta=0.8761, tau_s=1/beta=1.1414
%old| 
%old|  mdpr1=PLE_exponent(par1);resid1=edpr-mdpr1;mrse1=std(resid1,1),cc=corrcoef(mdpr1,edpr);cc(2)^2
%old| mrse1 = 0.0825 ; R2 =  0.9123


%%%% Parallel learning hypothesis: general curve + specific curve
%old|  par2=lsqcurvefit('PLE_2exp',[2 1 .1 .5 .5],[T Tc],edpr),par2=lsqcurvefit('PLE_2exp',par2,[T Tc],edpr)
%old| 
%old| par2:  F_asympt=1.2741, F_init=0.7844, alpha=0.0781, tau=1/alpha=12.8041
%old|        G_asympt=0.4752, beta=0.4251, tau_c=1/beta=2.3524
%old| 
%old|  mdpr2=PLE_2exp(par2);resid2=edpr-mdpr2;mrse2=std(resid2,1),cc=corrcoef(mdpr2,edpr);cc(2)^2
%old| mrse2 = 0.1247 ; R2 = 0.7993



%%%%%  Test whether the increase of H1 relative to H0 is significant
%%
%%  Adopt the formula for nested multiple regression models:
 cc=corrcoef(edpr,mdpr1);R2f=cc(2)^2, cc=corrcoef(edpr,mdpr0);R2r=cc(2)^2

R2f=0.9123 ; R2r = 0.7425

 F=((32-5-1)*(R2f-R2r))/((5-3)*(1-R2f)), 1-fcdf(F,2,26)

F(2,26) = 25.1460,  significance = 8.3643e-007 < 10e-6


%% Ditto for the independence hypothesis
%old|  cc=corrcoef(edpr,mdpr2);R2f_=cc(2)^2,F_=((32-5-1)*(R2f_-R2r))/((5-3)*(1-R2f_)),1-fcdf(F_,2,26),clear R2f_ F_
%old| 
%old| R2f_ = 0.7993 ; F_ = 3.6796 ; signif = 0.0392



%%%%%  Plot the data and H1 and H2

 plot(blk1,edpr(blk1),'ro',blk2,edpr(blk2),'bs',2:9,mdpr1(2:9),'kx-',2:9,mdpr2(2:9),'k.-');
 axis([0 33 0 2.5]);set(gca,'xtick',[0:2:32],'xticklabel','0||4||8||12||16||20||24||28||32');
 xlabel('Block number');ylabel('d''');legend('Data, context A','Data, context B','Switch-cost hyp.','Independence hyp.',4);
 hold on;plot(1,mdpr1(1),'kx',10:17,mdpr1(10:17),'kx-',18:25,mdpr1(18:25),'kx-',26:31,mdpr1(26:31),'kx-',32,mdpr1(32),'kx');
 plot(1,mdpr2(1),'k.',10:17,mdpr2(10:17),'k.-',18:25,mdpr2(18:25),'k.-',26:31,mdpr2(26:31),'k.-',32,mdpr2(32),'k.');hold off;
% Set page size=[1 3.25 6.5 4.5]; Marker size for data=8; Line width=2 and 1.


%%%%%  Raw data, 04/03/2003
 [T Ts Tc , edpr mdpr1 mdpr2 mdpr0]

 T  Ts  Tc    empir.   sw.cost  indep.   null
----------------------------------------------
 1 Inf   1    0.8866    0.9669  0.7844  0.8307

 2   1   1    0.7737    0.7160  0.8212  0.9307
 3   2   2    0.9890    0.9825  1.0198  1.0194
 4   3   3    1.1810    1.1292  1.1588  1.0980
 5   4   4    1.2232    1.2228  1.2582  1.1678
 6   5   5    1.4224    1.2914  1.3311  1.2298
 7   6   6    1.2681    1.3468  1.3861  1.2847
 8   7   7    1.3896    1.3943  1.4287  1.3334
 9   8   8    1.3713    1.4362  1.4629  1.3767

10   1   2    1.0998    1.1428  1.1962  1.4150
11   2   3    1.3254    1.3701  1.3220  1.4490
12   3   4    1.4194    1.4812  1.4091  1.4792
13   4   5    1.4515    1.5425  1.4707  1.5060
14   5   6    1.6663    1.5817  1.5151  1.5298
15   6   7    1.7963    1.6104  1.5481  1.5508
16   7   8    1.5526    1.6337  1.5733  1.5695
17   8   9    1.7668    1.6536  1.5931  1.5861

18   1   9    1.3599    1.3403  1.6036  1.6008
19   2  10    1.6280    1.5494  1.6188  1.6139
20   3  11    1.5417    1.6441  1.6314  1.6255
21   4  12    1.6677    1.6904  1.6421  1.6357
22   5  13    1.7359    1.7160  1.6514  1.6448
23   6  14    1.6360    1.7324  1.6595  1.6529
24   7  15    1.6175    1.7444  1.6668  1.6601
25   8  16    1.7880    1.7542  1.6733  1.6665

26   1  10    1.4582    1.4316  1.6694  1.6721
27   2  11    1.6135    1.6324  1.6782  1.6771
28   3  12    1.7350    1.7194  1.6854  1.6816
29   4  13    1.7868    1.7588  1.6914  1.6855
30   5  14    1.9867    1.7781  1.6965  1.6890
31   6  15    1.7061    1.7888  1.7010  1.6921

32   1  17    1.4064    1.4661  1.7052  1.6948
----------------------------------------------


%%%%%%%%%%  PLOT ALL SIX LEARNING CURVES  %%%%%%%%%%%%%%%%%%%
%%
 z=HUMAN_ZCONGR;subplot(2,1,1);hold on;
 idx=[2:9]';plot(idx,z(2,idx),'b^-',idx,z(4,idx),'ro-',idx,z(6,idx),'gs-');
 idx=[10:17]';plot(idx,z(2,idx),'b^-',idx,z(4,idx),'ro-',idx,z(6,idx),'gs-');
 idx=[18:25]';plot(idx,z(2,idx),'b^-',idx,z(4,idx),'ro-',idx,z(6,idx),'gs-');
 idx=[26:31]';plot(idx,z(2,idx),'b^-',idx,z(4,idx),'ro-',idx,z(6,idx),'gs-');
 idx=[1 32]';plot(idx,z(2,idx),'b^',idx,z(4,idx),'ro',idx,z(6,idx),'gs');
 title('Incongruent stimuli');axis([0 33 -.5 2]);set(gca,'box','on','xtick',[0:4:32],'ytick',[-.5:.5:2]);ylabel('z-probability correct');

 subplot(2,1,2);plot(1:32,z(1,:),'b^',1:32,z(3,:),'ro',1:32,z(5,:),'gs');
 title('Congruent stimuli');axis([0 33 -.5 2]);set(gca,'box','on','xtick',[0:4:32],'ytick',[-.5:.5:2]);ylabel('z-probability correct');
 xlabel('Block number');legend('Contrast 0.245','Contrast 0.160','Contrast 0.106',3);
% Set page size=[1.25 1 6 9]; Line width=1; Open symbols.

%% Now the fun part --- the inset
 m=mean(z'),h=axes('position',[0.69 0.12 0.2 0.14]);plot(1:2,m([2 1]),'b^-',1:2,m([4 3]),'ro-',1:2,m([6 5]),'gs-');
 axis([.5 2.5 -.25 1.5]);set(gca,'xtick',[1 2],'xticklabel','Incon|Congr','XAxisLocation','top','ytick',[0 .5 1],'yticklabel','0||1');

m = [0.8210  1.1264 , 0.9559  0.5653 , 1.0308 -0.0696]



%%%%%%%%%%  Subject-by-subject analyses  %%%%%%%%%%%%%%%%%%%

 edpr13=zeros(32,13);for k=1:13 edpr13(:,k)=mean(S(k).dprime,2);end
 [par1_13,rmse1_13] = PLE_curvefit('PLE_exponent',edpr13,par1);[par1_13 rmse1_13]
 foo=[par1_13 rmse1_13];foo=[[1:13]' foo(:,[1 2]) foo(:,1).*(1-foo(:,2)) foo(:,[3 4]) foo(:,1).*foo(:,4) foo(:,[5 6])]

 S   F_asym     g     F_init     tau        s      s*D      tau_s     RMSE
---------------------------------------------------------------------------
 1   2.4237   0.8336  0.4033    6.3231    0.1617  0.3919    2.4906   .2627
 2   1.7975   0.5130  0.8754   16.8901    0.0879  0.1580    2.2744   .2616
 3   2.4419   0.6517  0.8505    1.3976    0.1048  0.2559    1.5292   .1725
 4   1.9444   0.1893  1.5764   10.1537    0.3221  0.6263    1.0089   .3072
 5   1.4431   0.7477  0.3641    0.3620    0.3165  0.4567    1.8962   .2629
 6   2.3437   0.7609  0.5603    9.8307    0.0769  0.1803    1.7448   .1951
 7   1.4537   0.2679  1.0643   20.3079    0.1778  0.2585    1.7313   .2251
 8   2.4009   0.5712  1.0296   21.6605    0.1391  0.3340    0.5444   .2885
 9   1.5245   0.7589  0.3676   28.6507    0.1420  0.2166    0.4414   .2038
10   1.1333   0.5121  0.5530   13.6100    0.2665  0.3020    0.8358   .2050
11   2.5163   0.3856  1.5460    8.9995    0.1375  0.3460    1.3163   .3132
12   1.6009   0.2515  1.1983   25.1856    0.3610  0.5779    1.6636   .2706
13   2.0791   0.6803  0.6647   14.4044    0.2035  0.4231    1.6111   .1896

 describe(foo)
          Mean  Std.dev       Min     Q25  Median     Q75     Max
-----------------------------------------------------------------
F_asym   1.931    0.470      1.13    1.51    1.94    2.41    2.52
g        0.548    0.217      0.19    0.36    0.57    0.75    0.83
F_init   0.850    0.417      0.36    0.52    0.85    1.10    1.58
tau     13.675    8.691      0.36    8.33   13.61   20.65   28.65
s        0.192    0.095      0.08    0.13    0.16    0.28    0.36
s*D      0.348    0.144      0.16    0.25    0.33    0.43    0.63
tau_s    1.468    0.621      0.44    0.97    1.61    1.78    2.49
RMSE     0.243    0.047      0.17    0.20    0.26    0.28    0.31


%old|  S  F_asym    F_init    alpha     swcost     beta      RMSE
%old| ------------------------------------------------------------
%old|  1  2.4253    0.4040    0.1584   -0.3930    0.3948    0.2627
%old|  2  2.0210    0.9167    0.0383   -0.1666    0.4667    0.2611
%old|  3  2.4420    0.8506    0.7158   -0.2560    0.6524    0.1725
%old|  4  1.9497    1.5802    0.0938   -0.6274    0.9920    0.3072
%old|  5  1.4433    0.3646    2.7543   -0.4568    0.5252    0.2629
%old|  6  2.3449    0.5604    0.1017   -0.1807    0.5606    0.1951
%old|  7  3.5150    1.0887    0.0045   -0.2616    0.6032    0.2225
%old|  8  3.7937    1.1046    0.0159   -0.3423    1.8391    0.2841
%old|  9  3.1508    0.4097    0.0107   -0.2308    2.3406    0.1976
%old| 10  1.1445    0.5580    0.0699   -0.3034    1.1910    0.2050
%old| 11  2.5180    1.5476    0.1105   -0.3467    0.7566    0.3132
%old| 12  3.6004    1.2101    0.0044   -0.5843    0.6204    0.2681
%old| 13  2.1109    0.6768    0.0656   -0.4255    0.6251    0.1895

%old|  describe([par1_13 rmse1_13])
%old| 
%old|           Mean  Std.dev       Min     Q25  Median     Q75     Max
%old| -----------------------------------------------------------------
%old| F_asym   2.497    0.817      1.14    2.00    2.43    3.24    3.79
%old| F_init   0.867    0.418      0.36    0.52    0.85    1.13    1.58
%old| alpha    0.319    0.755      0.00    0.01    0.07    0.12    2.75
%old| swcost  -0.352    0.143     -0.63   -0.43   -0.34   -0.25   -0.17
%old| beta     0.890    0.582      0.39    0.55    0.63    1.04    2.34
%old| RMSE     0.242    0.047      0.17    0.20    0.26    0.27    0.31


%old|  [par2_13,rmse2_13] = PLE_curvefit('PLE_2exp',edpr13,par2);describe([par2_13 rmse2_13])
%old| 
%old|           Mean  Std.dev       Min     Q25  Median     Q75     Max
%old| -----------------------------------------------------------------
%old| F_asym   1.759    1.147      0.00    1.02    1.86    2.29    3.74
%old| F_init   0.814    0.435      0.16    0.43    0.82    1.21    1.41
%old| alpha    0.098    0.162      0.00    0.01    0.04    0.12    0.61
%old| G_asym   0.649    0.540      0.00    0.18    0.47    1.06    1.69
%old| beta     0.462    0.329      0.01    0.15    0.45    0.71    1.06
%old| RMSE     0.253    0.059      0.18    0.21    0.24    0.29    0.38


% The switch-cost model is (marginally) better for 11 out of 13 subjects. Exceptions: subjects 1 and 2.
 find(rmse1_13 > rmse2_13)'
ans = [1 2]    
 plot(rmse1_13,rmse2_13,'o');axis([.1 .4 .1 .4]);axis('square');refline(1,0);describe(rmse2_13-rmse1_13)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   0.010    0.026     -0.05    0.00    0.01    0.02    0.07



%%%%%%%%%%  By difficulty level, April 3, 2003  %%%%%%%%%%%%%%%%%%%

%%%%%%  Each level fitted separately, 15=3*5 parameters total
%
 d=HUMAN_ZCONGR';edpr3=[d(:,1)+d(:,2),d(:,3)+d(:,4),d(:,5)+d(:,6)];clear d
 [par1_3,rmse1_3] = PLE_curvefit('PLE_exponent',edpr3,par1);
 foo=par1_3;foo=[[.23;.15;.10] foo(:,[1 2]) foo(:,1).*(1-foo(:,2)) foo(:,[3 4]) foo(:,1).*foo(:,4) foo(:,5) rmse1_3]

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.3941  0.4343  1.3543    9.9001    0.2032  0.4865  1.1881    0.1280
.15   1.8690  0.4964  0.9412    8.8682    0.1595  0.2980  1.4397    0.1073
.10   1.2684  0.5339  0.5912   13.4278    0.1690  0.2143  0.5999    0.1072


%old| Cp    F_asym  F_init  alpha     swcost   beta    1/alpha  1/beta     RMSE
%old| --------------------------------------------------------------------------
%old| .23   2.3977  1.3572  0.0997   -0.4866  0.8424   10.0305  1.1871    0.1280
%old| .15   1.8677  0.9400  0.1133   -0.2975  0.6945    8.8275  1.4400    0.1073
%old| .10   1.2946  0.6025  0.0676   -0.2167  1.6398   14.8035  0.6098    0.1071

%old|  [par2_3,rmse2_3] = PLE_curvefit('PLE_2exp',edpr3,par2);[par2_3 rmse2_3]
%old| 
%old| Cp    F_asym    F_init    alpha     G_asym     beta      RMSE
%old| --------------------------------------------------------------
%old| .23   1.5529    1.0794    0.0793    0.6970    0.4335    0.1877
%old| .15   1.3542    0.7829    0.1195    0.4090    0.3459    0.1372
%old| .10   0.9584    0.4809    0.0444    0.3199    0.4339    0.1230


 mdpr1_3=[PLE_exponent(par1_3(1,:)),PLE_exponent(par1_3(2,:)),PLE_exponent(par1_3(3,:))];
 cc=corrcoef([edpr3 mdpr1_3]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
R2 = [0.8715    0.8681    0.7639]   % variance accounted for


%%%%%%  Three levels fitted together with common tau, s, and tau_s 
%% 9 parameters total:  params=[D1 D2 D3 , g1 g2 g3 , tau s tau_s]
%
 foo = [2.5 2.0 1.5 , 0.5 0.5 0.5 , 10 .2 1] ;
 [par1_9,rmse1_9] = PLE_curvefit('PLE_exponent',edpr3(:),foo) ;
 [par1_9,rmse1_9] = PLE_curvefit('PLE_exponent',edpr3(:),par1_9)
par1_9 = [2.3932 1.9023 1.2230 , 0.4382 0.4854 0.5276 , 10.1879  0.1835 1.1772]
rmse1_9 = 0.1165

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.3932  0.4382  1.3446   10.1879    0.1835  0.4392  1.1772    0.1291
.15   1.9023  0.4854  0.9788   10.1879    0.1835  0.3491  1.1772    0.1096
.10   1.2230  0.5777  0.6453   10.1879    0.1835  0.2244  1.1772    0.1098

 mdpr1_9=reshape(PLE_exponent(par1_9),32,3);
 foo=(edpr3-mdpr1_9).^2;sqrt(mean(foo))
RMSE = [0.1291  0.1096  0.1098]
 cc=corrcoef([edpr3 mdpr1_9]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
32-point R2s = [0.8696  0.8630  0.7524]
 cc=corrcoef([edpr3(:) mdpr1_9(:)]);R2r=cc(2).^2
96-point R2 = 0.9459


%%%%%  Test whether the 9-parameter model is superior to the 15-parameter one
 cc=corrcoef(edpr3(:),mdpr1_3(:));R2f=cc(2)^2

R2f=0.9477 ; R2r = 0.9459

 F=((96-15-1)*(R2f-R2r))/((15-9)*(1-R2f)), 1-fcdf(F,6,80)

F(6,80) = 0.4536,  significance = 0.8404


%%%%%%  Constrain g across the three difficulty levels 
%% 7 parameters total:  params=[D1 D2 D3 , g , tau s tau_s]
%
 foo = [2.5 2.0 1.5 , 0.5 , 10 .2 1] ;
 [par1_7,rmse1_7] = PLE_curvefit('PLE_exponent',edpr3(:),foo) ;
 [par1_7,rmse1_7] = PLE_curvefit('PLE_exponent',edpr3(:),par1_7)
par1_7 = [2.4145 1.8892 1.1974 , 0.4660 , 10.1551  0.1837  1.1784]
rmse1_7 = 0.1181

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.4145  0.4660  1.2894   10.1551    0.1837  0.4435  1.1784    0.1308
.15   1.8892  0.4660  1.0089   10.1551    0.1837  0.3470  1.1784    0.1102
.10   1.1974  0.4660  0.6395   10.1551    0.1837  0.2200  1.1784    0.1124

 par1_7([1 2 3]) ./ [.23 .15 .10]
ans = [10.4978  12.5947  11.9740]     % simple proportionality?
 plot([.106 .160 .245],par1_7([3 2 1]),'bo-',[0 .106],[0 par1_7(3)],'b:');axis([0 0.3 0 3]);
 set(gca,'xtick',[.106 .160 .245]);xlabel('Stimulus contrast');ylabel('Asymptotic d''');title('Empirical group data,  7-parameter switch cost fits')

 mdpr1_7=reshape(PLE_exponent(par1_7),32,3);
 foo=(edpr3-mdpr1_7).^2;sqrt(mean(foo))
RMSE = [0.1308  0.1102  0.1124]
 cc=corrcoef([edpr3 mdpr1_7]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
32-point R2s = [0.8672  0.8611  0.7485]
 cc=corrcoef([edpr3(:) mdpr1_7(:)]);R2r=cc(2).^2
96-point R2 = 0.9444

%%%%%  Test whether the 7-parameter model is superior to the 9-parameter one

R2f=0.9459 ; R2r = 0.9444

 F=((96-9-1)*(R2f-R2r))/((9-7)*(1-R2f)), 1-fcdf(F,2,86)

F(2,86) = 1.2062,  significance = 0.3043


**********  SUMMARY TABLE : EMPIRICAL GROUP DATA : 3 APRIL 2003  ***********

 Cp   F_asym    g     F_init    tau       s      s*D    tau_s     RMSE   R2
----------------------------------------------------------------------------
Null hypothesis, 3 parameters
---   1.7165  0.5161  0.8306    8.3497   ---     ---     ---     .1413 .7425

Switch-cost hypothesis, 5 parameters
---   1.8415  0.4749  0.9670   10.3764  0.1799  0.3313  1.1393   .0825 .9123

Independence hypothesis, 5 parameters            C      tau_c
---   1.2741  0.3843  0.7844   12.8041   ---    0.4752  2.3524   .1247 .7993


---------------- Switch-cost fits at separate contrast levels --------------
Average curve, 5 parameters
---   1.8415  0.4749  0.9670   10.3764  0.1799  0.3313  1.1393   .0825 .9123

15-parameter model,R2=0.9477
.23   2.3941  0.4343  1.3543    9.9001  0.2032  0.4865  1.1881   .1280 .8715
.15   1.8690  0.4964  0.9412    8.8682  0.1595  0.2980  1.4397   .1073 .8681
.10   1.2684  0.5339  0.5912   13.4278  0.1690  0.2143  0.5999   .1072 .7639

9-parameter model, R2=0.9459
.23   2.3932  0.4382  1.3446   10.1879  0.1835  0.4392  1.1772   .1291 .8696
.15   1.9023  0.4854  0.9788   10.1879  0.1835  0.3491  1.1772   .1096 .8630
.10   1.2230  0.5276  0.5777   10.1879  0.1835  0.2244  1.1772   .1098 .7524

7-parameter model, R2=0.9444
.23   2.4145  0.4660  1.2894   10.1551  0.1837  0.4435  1.1784   .1308 .8672
.15   1.8892  0.4660  1.0089   10.1551  0.1837  0.3470  1.1784   .1102 .8611
.10   1.1974  0.4660  0.6395   10.1551  0.1837  0.2200  1.1784   .1124 .7485
----------------------------------------------------------------------------
 Cp   F_asym    g     F_init    tau       s      s*D    tau_s     RMSE   R2




%%%%%  Plot the data by difficulty level, with switch-cost curves
 T1=[1, 2:.5:9, 10:.5:17, 18:.5:25, 26:.5:31, 32]';Ts1=[Inf, 1:.5:8, 1:.5:8, 1:.5:8, 1:.5:6, 1]';
% Ts1=[Inf Inf Inf, .5:.5:8.5, .5:.5:8.5 ,.5:.5:8.5 ,.5:.5:6.5, .5:.5:1.5]';T1=[0.5:.5:1.5, 1.5:.5:9.5, 9.5:.5:17.5, 17.5:.5:25.5, 25.5:.5:31.5, 31.5:.5:32.5]';
 m=reshape(PLE_exponent(par1_7,[T1,Ts1]),58,3);

 plot(T,edpr3(:,1),'b^',T,edpr3(:,2),'bo',T,edpr3(:,3),'bs',1,m(1,:),'kx-',2:.5:9,m(2:16,:),'k-',10:.5:17,m(17:31,:),'k-',18:.5:25,m(32:46,:),'k-',26:.5:31,m(47:57,:),'k-',T1([2 17 32 47 58]),m([2 17 32 47 58],:),'kx');
%%%%...2:9,mdpr1_3(2:9,:),'kx-',10:17,mdpr1_3(10:17,:),'kx-',18:25,mdpr1_3(18:25,:),'kx-',26:31,mdpr1_3(26:31,:),'kx-',[1 32]',mdpr1_3([1 32],:),'kx')
 legend('Contrast 0.245','Contrast 0.160','Contrast 0.106','Switch cost fits',4);
 axis([0 33 0 3]);set(gca,'xtick',[0:4:32]);xlabel('Block number  (30 trials/block, 4 blocks/day)');ylabel('d''');
% Set page size=[1 3.25 6.5 4.5]; Line width=1; Open symbols.




%%%%%%%%%%  Response times,  7 Apr 2003 %%%%%%%%%%%%%%%%%%%
%
 clear all ; load S ; home
 T=[1:32]';Ts=[Inf,1:8,1:8,1:8,1:6,1]';
 rt3_13=zeros(32,3,13);for k=1:13 rt3_13(:,:,k)=squeeze(mean(S(k).tmean_RT,3));end;
 rt3=squeeze(mean(rt3_13,3));plot(rt3,'o-');axis([0 33 450 900]);
% Order: Cp=[.10 .15 .23]

 r=mean(rt3,2);round(r')
r=[798   810   779   742   723   676   656   628   626   629   605   598 ...
   594   579   583   586   564   585   579   559   575   563   567   561 ...
   551   576   549   557   555   548   544   575]


%% Baseline fit of the aggregate RT curve
 [RTpar0,rmse]=PLE_curvefit('PLE_exponent',r,[500 -0.5 10])
RTpar0 = [552.3163  -0.5063  6.4955]
rmse = 12.4162

RTpar0:  RT_asympt=552 msec, g=0.5063, RT_init=832 msec, tau=6.4955

 mRT0=PLE_exponent(RTpar0);cc=corrcoef([r mRT0]);cc(2).^2
R2 = 0.9729


%% Switch-cost fit of the aggregate RT curve
 [RTpar1,rmse]=PLE_curvefit('PLE_exponent',r,[500 -0.5 10 -.1 2])
RTpar1 = [499.4110  -0.6014  5.8382  -0.1286   22.8004]
rmse = 7.5849 msec
 RTpar1(1).*[(1-RTpar1(2)) RTpar1(4)*(1-exp(-7/RTpar1(5)))]
ans =  799.7608  -16.9783

RTpar1:  RT_asympt=500, g=0.6014, RT_init=800, tau=5.84
         s=0.1286, deltaRT=17, tau_s=22.8

 mRT1=PLE_exponent(RTpar1);cc=corrcoef([r mRT1]);cc(2).^2
R2 = 0.9899


%%%%%  Test whether the increase of H1 relative to H0 is significant
 cc=corrcoef(r,mRT1);rtR2f=cc(2)^2, cc=corrcoef(r,mRT0);rtR2r=cc(2)^2

rtR2f = 0.9899 ; rtR2r = 0.9729

 rtF=((32-5-1)*(rtR2f-rtR2r))/((5-3)*(1-rtR2f)), 1-fcdf(rtF,2,26)

rtF(2,26) = 21.9322,  significance = 2.6262e-006 < 10e-5


%%%%%  Switch-cost fits at three contrast levels
%
 [RTpar1_3,rmse]=PLE_curvefit('PLE_exponent',rt3(:),[500 500 500 -0.6 6 -.1 5])
RTpar3 = [500.1728  498.8331  492.6201 , -0.6065  5.8101  -0.1334  24.1767]
rmse = 8.6995

 RTpar1_3(1:3).*(1-RTpar1_3(4))
RT0 = [803.5433  801.3910  791.4097]

 RTpar1_3(1:3).*RTpar1_3(6)*(1-exp(-7/RTpar1_3(7)))
delta RT = -[16.7713  16.7264  16.5180]

 mRT1_3=reshape(PLE_exponent(RTpar1_3),32,3);
 foo=(rt3-mRT1_3).^2;sqrt(mean(foo))
RMSE = [9.3720  7.9029  8.7609]
 cc=corrcoef([rt3 mRT1_3]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
32-point R2s = [0.9852    0.9890    0.9861]
 cc=corrcoef([rt3(:) mRT1_3(:)]);R2r=cc(2).^2
96-point R2 = 0.9868




%%% OLD OLD OLD OLD OLD OLD OLD OLD OLD OLD OLD OLD OLD
%%%%%%%%%%  Response times,  Dec 2002 %%%%%%%%%%%%%%%%%%%

 rt=zeros(32,13);for k=1:13 foo=squeeze(mean(S(k).tmean_RT,3));rt(:,k)=mean(foo,2);end;r=mean(rt,2);
 round(r')
r=[798   810   779   742   723   676   656   628   626   629   605   598 ...
   594   579   583   586   564   585   579   559   575   563   567   561 ...
   551   576   549   557   555   548   544   575]

%% Null hypothesis
 [RTpar0,RTrmse0]=PLE_curvefit('PLE_2exp',r,[550 850 .2],T);[RTpar0,RTrmse0]

RTpar0:  RT_asympt=552, RT_init=832, alpha=0.1539, tau=1/alpha=6.4977

 mRT0=PLE_2exp(RTpar0,T);resRT0=r-mRT0;std(resRT0,1)
RTmrse0 = 12.4161

 plot(blk1,resRT0(blk1),'ro',blk2,resRT0(blk2),'bs');axis([0 33 -40 40]);set(gca,'xtick',[0:4:32]);refline(0,0);


%% Switch cost hypothesis
 [RTpar1,RTrmse1]=PLE_curvefit('PLE_exponent',r,[550 850 .2 +30 .1]);[RTpar1,RTrmse1]

RTpar1:  RT_asympt=516, RT_init=808, alpha=0.1631, tau=1/alpha=6.1312
         swcost=+46.62, beta=0.0676, tau_s=1/beta=14.79

 mRT1=PLE_exponent(RTpar1);resRT1=r-mRT1;std(resRT1,1)
RTmrse1 = 7.9591

 (1-exp(-7/14.79)).*[1 46.62]     % the real switch cost is about 17.6=46.6*0.37
    0.3771   17.5782

 plot(blk1,resRT1(blk1),'ro',blk2,resRT1(blk2),'bs');axis([0 33 -40 40]);set(gca,'xtick',[0:4:32]);refline(0,0);
 plot(blk1,r(blk1),'ro',blk2,r(blk2),'bs',2:9,mRT1(2:9),'kx-',10:17,mRT1(10:17),'kx-',18:25,mRT1(18:25),'kx-',26:31,mRT1(26:31),'kx-',[1,32],mRT1([1,32]),'kx');
 axis([0 33 450 900]);set(gca,'xtick',[0:4:32],'ytick',[500:50:900]);xlabel('Block number');ylabel('Trimmed mean RT  [msec]');
% Set page size=[1 3.25 6.5 4.5]; Marker size for data=8; Line width=2; Set property 'MarkerFaceColor'=[1 0 0] and [0 0 1].



%% The independence hypothesis does not yield any improvement over the null hypothesis.
 [RTpar2,RTrmse2]=PLE_curvefit('PLE_2exp',r,[550 850 .2 +30 .1]);[RTpar2,RTrmse2]

RTpar2:  RT_asympt=552, RT_init=832, alpha=0.1537, tau=1/alpha=6.5062
         swcost=+28.67, beta=0.0006, tau_s=1/beta=1622
RTmrse2 = 12.4160

%% If we impose a lower bound lb=0.05 on beta (by editing PLE_CURVEFIT.M, line 47)
%% the switch cost parameter becomes zero!!
RTpar2:  RT_asympt=552, RT_init=832, alpha=0.1538, tau=1/alpha=6.5020
         swcost=00.00, beta=0.050*, tau_s=1/beta=20*
RTmrse2 = 12.4161


%%%%%  Test whether the increase of H1 relative to H0 is significant
 cc=corrcoef(r,mRT1);rtR2f=cc(2)^2, cc=corrcoef(r,mRT0);rtR2r=cc(2)^2

rtR2f = 0.9888 ; rtR2r = 0.9729

 rtF=((32-5-1)*(rtR2f-rtR2r))/((5-3)*(1-rtR2f)), 1-fcdf(rtF,2,26)

rtF(2,26) = 18.6396,  significance = 9.5117e-006 < 10e-5



*************************************************************************
*************************************************************************
*****                                                               *****
*****    M O D E L    D - P R I M E     F I T S       04/19/2003    *****
*****                                                               *****
*************************************************************************
*************************************************************************

 cd PLModel1 ; load final_fits

 fdpr=mean(Z_full)'.*2;rdpr=mean(Z_reduc)'.*2;edpr=mean(HUMAN_ZCONGR)'.*2;
 blk1=[1,10:17,26:31]';blk2=[2:9,18:25,32]';

%% Two time variables: Clock time T (in blocks) and time since last switch Ts
 T=[1:32]';Ts=[Inf,1:8,1:8,1:8,1:6,1]';


%%%% Switch-cost hypothesis: general curve + switch costs
 fpar1=lsqcurvefit('PLE_exponent',[2 0.5 10 .2 2],[T Ts],fdpr),fpar1=lsqcurvefit('PLE_exponent',fpar1,[T Ts],fdpr)
fpar1 = [1.8064  0.4499  7.9331  0.1531  2.1577]
  F_asympt=1.8064, F_init=D*(1-g)=0.9938, g=0.4499, tau=7.9331
  s=0.1531, s*D=0.2765, tau_s=2.1577

 mfdpr1=PLE_exponent(fpar1);resid=edpr-mfdpr1;mfrse1=std(resid,1),cc=corrcoef(mfdpr1,edpr);cc(2)^2
mfrse1 = 0.0930 ; R2 = 0.8950

% The same for the reduced model (that is, when the running-averaging is switched off)
 rpar1=lsqcurvefit('PLE_exponent',[2 0.5 10 .2 2],[T Ts],rdpr),rpar1=lsqcurvefit('PLE_exponent',rpar1,[T Ts],rdpr)
rpar1 = [1.7573  0.4020  6.7641  0.1883  2.7163]
  F_asympt=1.7573, F_init=D*(1-g)=1.0509, g=0.4020, tau=6.7641
  s=0.1883, s*D=0.3309, tau_s=2.7163

 mrdpr1=PLE_exponent(rpar1);resid=edpr-mrdpr1;mrrse1=std(resid,1),cc=corrcoef(mrdpr1,edpr);cc(2)^2
mrrse1 = 0.1077 ; R2 = 0.8729


%%%%%  Plot the data and F1
 plot(blk1,edpr(blk1),'ro',blk2,edpr(blk2),'bs',2:9,mfdpr1(2:9),'kx-',blk1,fdpr(blk1),'r.',blk2,fdpr(blk2),'b.');
 axis([0 33 0 2.5]);set(gca,'xtick',[0:2:32],'xticklabel','0||4||8||12||16||20||24||28||32');
 xlabel('Block number');ylabel('d''');legend('Data, context A','Data, context B','Switch-cost hyp.','Full model, A', 'Full model, B',4);
 hold on;plot(1,mfdpr1(1),'kx',10:17,mfdpr1(10:17),'kx-',18:25,mfdpr1(18:25),'kx-',26:31,mfdpr1(26:31),'kx-',32,mfdpr1(32),'kx');hold off;
% Set page size=[1 3.25 6.5 4.5]; Marker size for data=8; Line width=2 and 1.



%%%%%%%%%   Three difficulty levels  %%%%%%%%%%%%%%%%%
%% 7 parameters total:  params=[D1 D2 D3 , g , tau s tau_s]
%
 d=HUMAN_ZCONGR';edpr3=[d(:,1)+d(:,2),d(:,3)+d(:,4),d(:,5)+d(:,6)];
 d=Z_full';fdpr3=[d(:,1)+d(:,2),d(:,3)+d(:,4),d(:,5)+d(:,6)];
 d=Z_reduc';rdpr3=[d(:,1)+d(:,2),d(:,3)+d(:,4),d(:,5)+d(:,6)];clear d

 foo = [2.5 2.0 1.5 , 0.5 , 10 .2 1] ;
 [fpar3,frmse3]=PLE_curvefit('PLE_exponent',fdpr3(:),foo),[fpar3,frmse3]=PLE_curvefit('PLE_exponent',fdpr3(:),fpar3)
fpar3 = [2.4806  1.7999  1.1300  0.4519  8.0361  0.1484  2.0552]
frmse3 = 0.0333

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.4806  0.4519  1.3596    8.0361    0.1484  0.3681  2.0552    0.1810
.15   1.7999  0.4519  0.9865    8.0361    0.1484  0.2671  2.0552    0.1171
.10   1.1300  0.4519  0.6194    8.0361    0.1484  0.1677  2.0552    0.1232

 mfdpr3=reshape(PLE_exponent(fpar3),32,3);foo=(edpr3-mfdpr3).^2;sqrt(mean(foo))
RMSE = [0.1810  0.1171  0.1232]
 cc=corrcoef([edpr3 mfdpr3]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
32-point R2s = [0.8415  0.8652  0.7273]
 cc=corrcoef([edpr3(:) mfdpr3(:)]);cc(2).^2
96-point R2 = 0.9263


% The same for the reduced model (that is, when the running-averaging is switched off)
 foo = [2.5 2.0 1.5 , 0.5 , 10 .2 1] ;
 [rpar3,rrmse3]=PLE_curvefit('PLE_exponent',rdpr3(:),foo),[rpar3,rrmse3]=PLE_curvefit('PLE_exponent',rdpr3(:),rpar3)
rpar3 = [2.3865  1.7614  1.1215  0.4036  6.8696  0.1819  2.7435]
rrmse3 = 0.0471

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.3865  0.4036  1.4234    6.8696    0.1819  0.4342  2.7435    0.1527
.15   1.7614  0.4036  1.0505    6.8696    0.1819  0.3204  2.7435    0.1398
.10   1.1215  0.4036  0.6689    6.8696    0.1819  0.2040  2.7435    0.1331

 mrdpr3=reshape(PLE_exponent(rpar3),32,3);foo=(edpr3-mrdpr3).^2;sqrt(mean(foo))
RMSE = [0.1810  0.1171  0.1232]
 cc=corrcoef([edpr3 mrdpr3]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
32-point R2s = [0.8377  0.8434  0.6925]
 cc=corrcoef([edpr3(:) mrdpr3(:)]);cc(2).^2
96-point R2 = 0.9204

