
******  08 March 2005  -- ANOVA analyses on the PLExp2 data
******

 cd('C:\work\PLExp2');load estats18
 edpr18=cat(3,estats18(:).dprime);edpr18(37:48,:,:)=[];edpr=mean(edpr18,3);size(edpr)

%%% Average over time and contrast, to get the grand within-subject means
 d18=squeeze(mean(edpr18,2));d18=mean(d18)';plot(d18,'o')
 gr=[estats18(:).group]';xtab1(gr)
 anova1(d18,gr)
F(1,16) = 0.9151 , p = 0.3530 -- not significant

%%% Within-subject confidence intervals for the empirical d' data
% [36x3x18]: (block,Cpeak,sbj)
 S=reshape(1:18,[1,1,18]);S=repmat(S,[36,3,1]);
 gr=reshape([estats18(:).group],[1 1 18]);gr=repmat(gr,[36,3,1]);
 Cp=reshape([.23 .15 .10],[1 3 1]);Cp=repmat(Cp,[36,1,18]);
 bl=reshape(1:36,[36,1,1]);bl=repmat(bl,[1,3,18]);

 [SS,df,MS,lbl] = anova(edpr18(:),[S(:),bl(:),Cp(:)],'SBC');

Source         SumSq   eta2[%]     df        MeanSq
---------------------------------------------------
   S         284.691    21.40      17       16.7465
   B         132.570     9.97      35        3.7877
   C         476.173    35.80       2      238.0863
  SB         221.496    16.65     595        0.3723
  SC          59.437     4.47      34        1.7481
  BC          19.645     1.48      70        0.2806
 SBC         136.185    10.24    1190        0.1144   % <-- error term
Totl        1330.196   100.00    1943        0.6846
---------------------------------------------------

% Within-subject confidence interval based on MSerr = MSsbc above.
% See Eq.2 in Loftus & Masson (1994). Psychon Bull & Review, 1(4), 476-490.
sqrt(0.1144/18) = 0.0797    % RMSE over Ss
tinv([.025 .975],1190) = [-1.9620 1.9620]
CI95 = mean +/- .1564  , where .1564=1.962*.0797

% Alternative confidence interval based on the 5-parameter switch-cost fits
rmse = .1213 ; df = 108-5 ; tinv(.975,103)=1.983 ; CI95 = +/-.241

% Yet another alternative: Average all across subjects and look at
% the time-by-contrast interaction
 Cp=repmat([.23 .15 .10],36,1);bl=repmat([1:36]',1,3);
 [SS,df,MS,lbl] = anova(edpr(:),[bl(:),Cp(:)],'BC');

Source         SumSq   eta2[%]     df        MeanSq      F
---------------------------------------------------------------
   B           7.365    21.10      35        0.2104     13.4967
   C          26.454    75.78       2       13.2270    848.3693
  BC           1.091     3.13      70        0.0156
Totl          34.910   100.00     107        0.3263
---------------------------------------------------

% rmse = sqrt(MSerr) = sqrt(.0156) = .1249
% CI95 = sqrt(MS(3)) * tinv(.975,70) = .2490


%%  Confidence intervals for the z-Probability profiles, using the last method
 zPcongr18=cat(3,estats18(:).zP_congr);zPcongr18(37:48,:,:)=[]; % [36x3x18]
 zPincon18=cat(3,estats18(:).zP_incon);zPincon18(37:48,:,:)=[]; % [36x3x18]
 zP=[mean(zPcongr18,3),mean(zPincon18,3)];plot(zP,'.-')         % [36x6]

 Cp=repmat([.23 .15 .10 -.23 -.15 -.10],36,1);bl=repmat([1:36]',1,6);
 [SS,df,MS,lbl] = anova(zP(:),[bl(:),Cp(:)],'BC');

Source         SumSq   eta2[%]     df        MeanSq
---------------------------------------------------
   B           3.683     3.69      35        0.1052
   C          92.632    92.80       5       18.5264
  BC           3.503     3.51     175        0.0200
Totl          99.817   100.00     215        0.4643
---------------------------------------------------

% rmse = sqrt(MSerr) = sqrt(.0200) = .1415
% CI95 = sqrt(MS(3)) * tinv(.975,175) = .2792


%% Calculate the exact p value of the reversal for congruent stimuli.
 zPc=mean(zPcongr18,3);plot(zPc,'.-')         % [36x3]
 Cp=repmat([.23 .15 .10],36,1);bl=repmat([1:36]',1,3);
 [SS,df,MS,lbl] = anova(zPc(:),[bl(:),Cp(:)],'BC');

Source         SumSq   eta2[%]     df        MeanSq
---------------------------------------------------
   B           3.601    80.00      35        0.1029
   C           0.600    13.32       2        0.2999
  BC           0.300     6.67      70        0.0043
Totl           4.501   100.00     107        0.0421
---------------------------------------------------

 F=MS(2)/MS(3),p = 1-fcdf(F,df(2),df(3))
F = 69.9113 , p < 1e-16



***************************************************************************
***  March 24, 2005:  Copy of Congruent/Incongruent results
***  See C:\work\PLExp2\modelfits\Hebb_fits.txt for details.
***  See also C:\work\PLExp2\programs\PLE2_stats.m.
***  Original analysis conducted on April 8, 2004
***

 cd('C:\work\PLExp2\modelfits');load HUMAN_ZCONGR; MODEL_ZCONGR;
 sumsq2 = PLM_sumsqfit2(MODEL_ZCONGR,HUMAN_ZCONGR)
sumsq2 = 7.5112 = 0.8346 * 9 sessions

 mZ=MODEL_ZCONGR';describe(mZ);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % model, 500+500 runs
------------------------------------------------------------
   1.223    0.189      0.91    1.11    1.22    1.38    1.60    % congr, Cp=.23
   1.069    0.250      0.47    0.90    1.11    1.24    1.51    % incon, Cp=.23
   1.319    0.169      1.06    1.22    1.32    1.44    1.67    % congr, Cp=.15
   0.187    0.176     -0.24    0.07    0.20    0.32    0.50    % incon, Cp=.15
   1.364    0.153      1.13    1.28    1.35    1.49    1.68    % congr, Cp=.10
  -0.477    0.131     -0.78   -0.55   -0.47   -0.39   -0.27    % incon, Cp=.10
------------------------------------------------------------
   0.781    0.178      0.42    0.67    0.79    0.91    1.11

 mdpr=mZ(:,[1 3 5])+mZ(:,[2 4 6]);describe(mdpr)
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % model, 500+500 runs
------------------------------------------------------------
   2.293    0.274      1.72    2.10    2.31    2.52    2.69    % Cp=.23
   1.507    0.191      1.13    1.36    1.53    1.66    1.78    % Cp=.15
   0.887    0.122      0.64    0.79    0.90    0.99    1.10    % Cp=.10
------------------------------------------------------------
   1.562    0.196      1.16    1.42    1.58    1.72    1.86

 mZ=MODEL_ZCONGR';describe(mZ);eZ=HUMAN_ZCONGR';describe(eZ);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % empirical, 9+9 observers
------------------------------------------------------------
   1.208    0.190      0.79    1.08    1.20    1.37    1.50    % congr, Cp=.23
   0.963    0.250      0.25    0.87    1.00    1.12    1.33    % incon, Cp=.23
   1.300    0.200      0.66    1.17    1.28    1.45    1.63    % congr, Cp=.15
   0.284    0.143     -0.20    0.23    0.28    0.38    0.54    % incon, Cp=.15
   1.390    0.188      0.85    1.26    1.44    1.52    1.70    % congr, Cp=.10
  -0.432    0.106     -0.68   -0.52   -0.45   -0.35   -0.26    % incon, Cp=.10
------------------------------------------------------------
   0.785    0.179      0.28    0.68    0.79    0.91    1.07

 edpr=eZ(:,[1 3 5])+eZ(:,[2 4 6]);describe(edpr)
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % empirical, 9+9 observers
------------------------------------------------------------
   2.171    0.367      1.33    1.95    2.25    2.43    2.70    % Cp=.23
   1.584    0.277      0.79    1.45    1.65    1.78    1.96    % Cp=.15
   0.959    0.174      0.57    0.85    0.98    1.11    1.22    % Cp=.10
------------------------------------------------------------
   1.571    0.273      0.90    1.41    1.63    1.77    1.96


***** 24 March 2005 -- within-block confidence intervals for HUMAN_ZCONGR
%
 load('C:\work\PLExp2\modelfits\HUMAN_ZCONGR.mat');size(HUMAN_ZCONGR)
ans = [6 36]
 cond=repmat([1:6]',1,36);bl=repmat([1:36],6,1);
 [SS,df,MS,lbl] = anova(HUMAN_ZCONGR(:),[bl(:),cond(:)],'BC');

Source         SumSq   eta2[%]     df        MeanSq
---------------------------------------------------
   B           3.683     3.69      35        0.1052
   C          92.632    92.80       5       18.5264
  BC           3.503     3.51     175        0.0200
Totl          99.817   100.00     215        0.4643
---------------------------------------------------

 CI95 = sqrt(MS(3)/36) * tinv(.975,175)
CI95 = 0.0465


***** 24 March 2005 -- overall response asymmetry
%
 load('C:\work\PLExp2\D.mat');D(1)
ans =
         sbj: 1
       group: 1
        stim: [14489x7 double]
        data: [14489x10 double]
    miss_idx: [3x1 double]
    demo_idx: [38x1 double]
    skip_idx: [41x1 double]

%- Retrieve the relevant data for the first 9 sessions only.
%  See import_PLE2_data.m for details.
 idx=1:9*1200;resp=zeros(10800,18);Gabor_degr=zeros(10800,18);noise_degr=zeros(10800,18);
 for k=1:18 d=D(k).data;d(D(k).skip_idx,:)=[];resp(:,k)=d(idx,7);Gabor_degr(:,k)=d(idx,3);noise_degr(:,k)=d(idx,5);end
 xtab2(Gabor_degr(:),noise_degr(:))

         |   -15    15 | Total Percent
---------+-------------+--------------
     -10 | 48600 48600 | 97200    50.0
      10 | 48600 48600 | 97200    50.0
---------+-------------+--------------
 Total   | 97200 97200 |194400
 Percent |  50.0  50.0 |         100.0


 x=xtab2(resp(:),noise_degr(:))

         |   -15    15 | Total Percent
---------+-------------+--------------
       1 | 60966 32791 | 93757    48.2
       2 | 36234 64409 |100643    51.8
---------+-------------+--------------
 Total   | 97200 97200 |194400
 Percent |  50.0  50.0 |         100.0

 x=x./repmat(sum(x),2,1) , [x(1,1)+x(2,2) , x(1,2)+x(2,1)]./2

noise_degr   -15    +15               probability
---------------------------------------------------
resp="L"   .6272  .3374      resp=congr     .6449
resp="R"   .3728  .6626      resp=incongr   .3551
---------------------------------------------------


***** 24 March 2005 -- individual differences in the average d' levels
%
 load('C:\work\PLExp2\estats18.mat');
 edpr18=cat(3,estats18(:).dprime);edpr18(37:48,:,:)=[];edpr=squeeze(mean(edpr18,1))'
 describe(edpr)';Serr = [std(edpr) std(mean(edpr,2))]./sqrt(18)

Sbj   Cp=.23    Cp=.15    Cp=.10       Mean     NO FEEDBACK
--------------------------------------------
01    2.8279    1.8168    0.9087      1.8511
02    1.4841    0.9694    0.5326      0.9953
03    2.8728    2.1294    1.1224      2.0415
04    1.5369    1.2519    0.8558      1.2149
05    2.0548    1.4747    0.7717      1.4337
06    1.8214    1.4006    1.1163      1.4461
07    1.9477    1.4057    0.7844      1.3793
08    1.6696    1.3314    0.8009      1.2673
09    2.2168    1.6421    0.9981      1.6190
10    1.1522    0.8895    0.5601      0.8673
11    2.2306    1.8938    1.3247      1.8164
12    2.3859    1.8010    1.3308      1.8392
13    1.6927    1.2954    0.8804      1.2895
14    2.8381    2.0238    1.0832      1.9817
15    2.1513    1.3157    0.6076      1.3582
16    2.4173    1.7284    1.1081      1.7512
17    2.3357    1.7206    1.0199      1.6921
18    3.4361    2.4170    1.4484      2.4338
--------------------------------------------
mean  2.1707    1.5837    0.9586      1.5710
std   0.5784    0.3979    0.2635      0.3938
Serr  0.1363    0.0938    0.0621      0.0928   <-- Std error of the mean, b/n Ss

min   1.1522    0.8895    0.5326      0.8673
Q1    1.6927    1.3157    0.7844      1.2895
med   2.1841    1.5584    0.9534      1.5326
Q3    2.4173    1.8168    1.1163      1.8392
max   3.4361    2.4170    1.4484      2.4338
--------------------------------------------
Sbj   Cp=.23    Cp=.15    Cp=.10       Mean     NO FEEDBACK


%%% Compare with the feedback data from PLExp1 (Petrov, Dosher, Lu; Psych Review)
%
 clear all ; load('C:\work\PLModel1\estats.mat');estats
estats = 
       group: [13x1 double]
    zP_congr: [32x3x13 double]
    zP_incon: [32x3x13 double]
      dprime: [32x3x13 double]
       accZP: [13x6 double]
    contextc: [13x20 double]
    contexti: [13x20 double]
      otherc: [13x20 double]
      otheri: [13x20 double]

 dpr=squeeze(mean(estats.dprime,1))'; dpr=[dpr mean(dpr,2)]
 describe(dpr)';Serr = std(dpr)./sqrt(13)

Sbj   Cp=.23    Cp=.15    Cp=.10       Mean     TRIAL-BY-TRIAL FEEDBACK
--------------------------------------------
01    2.3915    1.8652    1.2699      1.8422
02    1.7434    1.3382    0.8581      1.3133
03    2.9422    2.2826    1.5869      2.2706
04    2.2745    1.7994    0.9752      1.6830
05    1.7272    1.3006    0.7415      1.2564
06    2.2912    1.7915    1.1170      1.7332
07    1.6822    1.1514    0.6820      1.1719
08    2.0116    1.7014    1.1068      1.6066
09    1.0005    0.7982    0.5361      0.7783
10    1.0813    0.8809    0.5521      0.8381
11    2.6676    2.3250    1.4408      2.1445
12    1.7401    1.1505    0.6861      1.1923
13    1.7624    1.3908    0.9425      1.3652
--------------------------------------------
mean  1.9474    1.5212    0.9612      1.4766
std   0.5641    0.4859    0.3321      0.4557
Serr  0.1565    0.1348    0.0921      0.1264   <-- Std error of the mean, b/n Ss

min   1.0005    0.7982    0.5361      0.7783
Q1    1.7159    1.1512    0.6851      1.1872
med   1.7624    1.3908    0.9425      1.3652
Q3    2.3162    1.8158    1.1552      1.7605
max   2.9422    2.3250    1.5869      2.2706
--------------------------------------------
Sbj   Cp=.23    Cp=.15    Cp=.10       Mean     TRIAL-BY-TRIAL FEEDBACK

% T test, independent samples
 t = (1.5710-1.4766)/sqrt(.0928^2+.1264^2)
t = 0.6020   % not significant

 tinv([.90 .95 .975],18+13-2)
ans =  1.3114  1.6991  2.0452




***************************************************************************
**********  Analysis of PLExp2 data  ---  Nov 10, 2003  *******************
**********

 !copy C:\work\PLExp2\data\sbj0??.dat C:\work\PLExp2\raw_data.dat
 D = import_PLE2_data('C:\work\PLExp2\raw_data.dat') ; st=PLE2_stats(D);

 idx=[1 3 4 6:18];d=mean(cat(3,st(idx).dprime),3); plot(d,'o');axis([0 49 0 3]);set(gca,'xtick',[0:4:48]);grid on
 md=mean(d,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:33,md(26:33),'ro-',34:36,md(34:36),'bs-',37:40,md(37:40),'bd-',41:44,md(41:44),'r^-',45:48,md(45:48),'k.-');axis([0 49 0 2.5]);set(gca,'xtick',[0:4:48]);grid on;ylabel('Mean d''');

 idx=[1 3 4 6:18];RT=mean(mean(cat(3,st(idx).RT),3),2);plot(RT,'o');axis([0 49 0 1200]);set(gca,'xtick',[0:4:48]);grid on;title('Mean of median RTs across Ss and conditions');


%% Individual plots
 figure(1);for k=1:18 subplot(6,3,k);plot(mean(st(k).dprime,2),'.-');axis([0 49 0 3]);set(gca,'xtick',[0:4:48],'xticklabel',[]);grid on;title(sprintf('Sbj %d  Group %d',st(k).sbj,st(k).group));if (rem(k,3)==1) ylabel('d''');else set(gca,'yticklabel',[]);end;end
 figure(2);for k=1:18 subplot(6,3,k);plot(mean(st(k).RT,2),'.-');axis([0 49 0 1200]);set(gca,'xtick',[0:4:48],'xticklabel',[]);grid on;title(sprintf('Sbj %d  Group %d',st(k).sbj,st(k).group));if (rem(k,3)==1) ylabel('median RT');else set(gca,'yticklabel',[]);end;end


%% By difficulty level
 plot(1,d(1,:),'ro',2:9,d(2:9,:),'bs-',10:17,d(10:17,:),'ro-',18:25,d(18:25,:),'bs-',26:33,d(26:33,:),'ro-',34:36,d(34:36,:),'bs-',37:40,d(37:40,:),'bd-',41:44,d(41:44,:),'r^-',45:48,d(45:48,:),'k.-');axis([0 49 0 3]);set(gca,'xtick',[0:4:48]);grid on;ylabel('Mean d''');


%%%%   d' fits       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
 edpr = mean(cat(3,st(:).dprime),3); edpr(37:48,:)=[];size(edpr)

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

%%% Null hypothesis: the context has no effect whatsoever
 edpr1=mean(edpr,2);
 par0=lsqcurvefit('PLE_exponent',[2 0.5 10],[T Ts],edpr1), par0=lsqcurvefit('PLE_exponent',par0,[T Ts],edpr1)
par0 = [1.7921  0.4566  9.4381]
  F_asymp=1.7921, F_init=D*(1-g)=0.9739, g=0.4566, tau=9.4381
  % PLExp1 fits: F_asympt=1.7165, F_init=D*(1-g)=0.8306, g=0.5161, tau=8.3497

 mdpr0=PLE_exponent(par0,[T Ts]);resid0=edpr1-mdpr0;mrse0=std(resid0,1),cc=corrcoef(mdpr0,edpr1);cc(2)^2
mrse0 = 0.1405 ; R2 = 0.7105           % PLExp1 fits: 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],edpr1),par1=lsqcurvefit('PLE_exponent',par1,[T Ts],edpr1)
par1 = [1.8895  0.4227  10.6962  0.1546  1.4360]
  F_asympt=1.8895, F_init=D*(1-g)=1.0908, g=0.4227, tau=10.6962
  s=0.1546, s*D=0.2922, tau_s=1.4360
  % PLExp1 fits: 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,[T Ts]);resid1=edpr1-mdpr1;mrse1=std(resid1,1),cc=corrcoef(mdpr1,edpr1);cc(2)^2
mrse1 = 0.1016 ; R2 = 0.8488
  % PLExp1 fits: mrse1 = 0.0825 ; R2 = 0.9123


%%%  1000 bootstrap fits, resampling the 18 subjects:
%    (See PLExp2/var_estim.txt for details.)

    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.909    0.113      1.53    1.83    1.91    1.99    2.27  % D
   0.429    0.050      0.23    0.40    0.43    0.46    0.59  % g
  11.310    3.450      1.45    8.96   10.65   12.99   29.98  % tau
   0.153    0.025      0.08    0.14    0.15    0.17    0.25  % s
   1.397    0.423      0.00    1.24    1.44    1.63    3.08  % tau_s



%%%%%%  Each level fitted separately, 15=3*5 parameters total, switch-cost hyp
 [par1_3,rmse1_3] = PLE_curvefit('PLE_exponent',edpr,par1,[T Ts]);
 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.6309  0.4298  1.5001   11.3965    0.1403  0.3691  1.5013    0.1417
.15   1.8997  0.4258  1.0908   10.8365    0.1537  0.2920  1.1931    0.1241
.10   1.1559  0.3984  0.6954    9.6645    0.1924  0.2224  1.8529    0.0877

% PLExp1 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
%  fits:--------------------------------------------------------------------------
%       .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

 mdpr1_3=[PLE_exponent(par1_3(1,:),[T Ts]),PLE_exponent(par1_3(2,:),[T Ts]),PLE_exponent(par1_3(3,:),[T Ts])];
 cc=corrcoef([edpr1 mdpr1_3]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
R2 = [0.8466  0.7937  0.7374]   % variance accounted for
   % PLExp1 fits: R2 = [0.8715  0.8681  0.7639]


%%%%%%  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',edpr(:),foo,[T Ts]) ;
 [par1_9,rmse1_9] = PLE_curvefit('PLE_exponent',edpr(:),par1_9,[T Ts])
par1_9 =[2.6235 1.9116 1.1494, 0.4290 0.4247 0.4047, 11.0850 0.1496 1.4218]
rmse1_9 = 0.1207

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.6235  0.4290  1.4981   11.0850    0.1496  0.3925  1.4218    0.1421
.15   1.9116  0.4247  1.0998   11.0850    0.1496  0.2860  1.4218    0.1243
.10   1.1494  0.4047  0.6842   11.0850    0.1496  0.1719  1.4218    0.0898

%PLExp1 par1_9 = [2.3932 1.9023 1.2230 , 0.4382 0.4854 0.5276 , 10.1879  0.1835 1.1772]
%fits:  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,[T Ts]),36,3);foo=(edpr-mdpr1_9).^2;sqrt(mean(foo))
RMSE = [0.1421  0.1243  0.0898]     % PLExp1 fits: RMSE = [0.1291  0.1096  0.1098]
 cc=corrcoef([edpr mdpr1_9]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
36-point R2s= [0.8459  0.7929  0.7261]    % PLExp1: 32-point R2s = [0.8696  0.8630  0.7524]
 cc=corrcoef([edpr(:) mdpr1_9(:)]);R2r=cc(2).^2
108-point R2r = 0.9549                    % PLExp1: 96-point R2 = 0.9459

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

 F=((108-15-1)*(R2f-R2r))/((15-9)*(1-R2f)), 1-fcdf(F,6,92)
F(6,92) = 0.1849,  significance = 0.9804


%%%%%%  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',edpr(:),foo,[T Ts]) ;
 [par1_7,rmse1_7] = PLE_curvefit('PLE_exponent',edpr(:),par1_7,[T Ts])

par1_7 = [2.6185 1.9106 1.1560, 0.4251, 11.0210  0.1498  1.4216]
rmse1_7 = 0.1208
%PLExp1 par1_7 = [2.4145 1.8892 1.1974 , 0.4660 , 10.1551  0.1837  1.1784]
%fits:  rmse1_7 = 0.1181

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.6185  0.4251  1.5054   11.0210    0.1498 0.3922   1.4216    0.1422
.15   1.9106  0.4251  1.0984   11.0210    0.1498 0.2862   1.4216    0.1243
.10   1.1560  0.4251  0.6646   11.0210    0.1498 0.1732   1.4216    0.0900

%PLExp1 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
%fits: --------------------------------------------------------------------------
%      .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

 mdpr1_7=reshape(PLE_exponent(par1_7,[T Ts]),36,3);foo=(edpr-mdpr1_7).^2;sqrt(mean(foo))
RMSE = [0.1422  0.1243  0.0900]   % PLExp1 fits: RMSE = [0.1308  0.1102  0.1124]
 cc=corrcoef([edpr mdpr1_7]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
36-point R2s = [0.8459  0.7929  0.7261]  % PLExp1: 32-point R2s = [0.8672  0.8611  0.7485]
 cc=corrcoef([edpr(:) mdpr1_7(:)]);R2r=cc(2).^2
108-point R2r = 0.9549            % PLExp1 fits: 96-point R2r = 0.9444

%%  The 7-parameter model is obviously superior to the 9-parameter one
R2f = R2r = 0.9549


%%%  1000 bootstrap fits, resampling the 18 subjects:
%    (See PLExp2/var_estim.txt for details.)

    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   2.646    0.188      2.03    2.52    2.65    2.77    3.39  % D1
   1.929    0.122      1.58    1.85    1.93    2.01    2.49  % D2
   1.164    0.070      0.95    1.12    1.17    1.21    1.41  % D3
   0.434    0.049      0.25    0.40    0.43    0.47    0.55  % g
  11.621    4.178      1.03    8.96   10.73   13.15   46.60  % tau
   0.149    0.024      0.08    0.13    0.15    0.16    0.23  % s
   1.399    0.376      0.01    1.24    1.42    1.59    2.70  % tau_s
------------------------------------------------------------
    Mean  Std.dev       Min     Q25  Median     Q75     Max



******  SUMMARY TABLE : EMPIRICAL GROUP DATA : NO FEEDBACK : NOV 2003  *****

 Cp   F_asym    g     F_init    tau       s      s*D    tau_s     RMSE   R2
----------------------------------------------------------------------------
Average curve, Null hypothesis, 3 parameters
---   1.7921  0.4566  0.9739    9.4381   ---     ---     ---     .1405 .7105

Average curve, Switch-cost hypothesis, 5 parameters
---   1.8895  0.4227  1.0908   10.6962  0.1546  0.2922  1.4360   .1016 .8488

---------------- Switch-cost fits at separate contrast levels --------------
15-parameter model,R2=0.9555
.23   2.6309  0.4298  1.5001   11.3965  0.1403  0.3691  1.5013   .1417 .8466
.15   1.8997  0.4258  1.0908   10.8365  0.1537  0.2920  1.1931   .1241 .7937
.10   1.1559  0.3984  0.6954    9.6645  0.1924  0.2224  1.8529   .0877 .7374

9-parameter model, R2=0.9549
.23   2.6235  0.4290  1.4981   11.0850  0.1496  0.3925  1.4218   .1421 .8459
.15   1.9116  0.4247  1.0998   11.0850  0.1496  0.2860  1.4218   .1243 .7929
.10   1.1494  0.4047  0.6842   11.0850  0.1496  0.1719  1.4218   .0898 .7261

7-parameter model, R2=0.9549
.23   2.6185  0.4251  1.5054   11.0210  0.1498 0.3922   1.4216   .1422 .8459
.15   1.9106  0.4251  1.0984   11.0210  0.1498 0.2862   1.4216   .1243 .7929
.10   1.1560  0.4251  0.6646   11.0210  0.1498 0.1732   1.4216   .0900 .7261

5-parameter model, R2=0.9545   29 Nov 2004
Fix tau=10.1551, tau_s=1.1784
.23   2.5819  0.4296  1.4727   10.1551* 0.1496 0.3922   1.1784*  .1433 .8431
.15   1.8840  0.4296  1.0746   10.1551* 0.1496 0.2862   1.1784*  .1242 .7933
.10   1.1400  0.4296  0.6502   10.1551* 0.1496 0.1732   1.1784*  .0902 .7225
----------------------------------------------------------------------------
 Cp   F_asym    g     F_init    tau       s      s*D    tau_s     RMSE   R2



*************************************************************************************
******
******  Re-analysis of PLExp2 data with fixed time constants ---  Nov 29, 2004
******  The M-file PLExp2/PLE_EXPONENT2 has the correct TAUs hard-wired.
******
 clear all ; load('C:\work\PLExp2\estats18.mat')
 edpr18=zeros(36,3,18);for k=1:18 edpr(:,:,k)=estats18(k).dprime(1:36,:);end
 edpr=mean(edpr18,3);plot(edpr,'.-');axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);
 T=[1:36]';Ts=[Inf,1:8,1:8,1:8,1:8,1:3]';xdata=[T Ts];


%%%%%%  Constrain g across the three difficulty levels.
%%%%%%  Fix TAU = 10.1551 ; TAU_S = 1.1784  -- see PLE_EXPONENT2
%%%%%%  These are the values reported in Table 1 of the Psych Reveiw paper.
%%%%%%
%% 5 parameters total:  params=[D1 D2 D3, g, s]
%
 par0 = [2.62 1.91 1.15, 0.42, .15] ;   % from the 7-par fit above
 [par,rmse] = PLE_curvefit('PLE_exponent2',edpr(:),par0,xdata)
 [par,rmse] = PLE_curvefit('PLE_exponent2',edpr(:),par,xdata)

par = [2.5819  1.8840  1.1400  0.4296  0.1496]
rmse = 0.1213
% par1_7 = [2.6185 1.9106 1.1560, 0.4251, 11.0210  0.1498  1.4216]
% rmse1_7 = 0.1208

 par(1:3).*(1-par(4))
ans = [1.4727  1.0746  0.6502]   % F_init = D * (1-g)

 Cp   F_asym    g     F_init    tau         s      s*D    tau_s      RMSE
--------------------------------------------------------------------------
.23   2.5819  0.4296  1.4727   10.1551*   0.1496 0.3922   1.1784*   0.1433
.15   1.8840  0.4296  1.0746   10.1551*   0.1496 0.2862   1.1784*   0.1242
.10   1.1400  0.4296  0.6502   10.1551*   0.1496 0.1732   1.1784*   0.0902


 mdpr=reshape(PLE_exponent2(par,xdata),36,3);foo=(edpr-mdpr).^2;sqrt(mean(foo))
RMSE = [0.1433  0.1242  0.0902]  % 5-parameter PLExp2 fits
% 7-parameter PLExp2 fits: RMSE = [0.1422  0.1243  0.0900]
% 7-parameter PLExp1 fits: RMSE = [0.1308  0.1102  0.1124]

 cc=corrcoef([edpr mdpr]);cc=[cc(1,4) cc(2,5) cc(3,6)].^2
36-point R2s = [0.8431  0.7933  0.7225]
% 7-parameter PLExp2: 32-point R2s = [0.8459  0.7929  0.7261]
% 7-parameter PLExp1: 32-point R2s = [0.8672  0.8611  0.7485]

 cc=corrcoef([edpr(:) mdpr(:)]);R2r=cc(2).^2
108-point R2r = 0.9545
% 7-parameter PLExp2 fits: 96-point R2r = 0.9549
% 7-parameter PLExp1 fits: 96-point R2r = 0.9444

%%  Test whether the 5-parameter model is superior to the 7-parameter one
 R2f = 0.9549 ; R2r = 0.9545 ;

 F=((108-7-1)*(R2f-R2r))/((7-5)*(1-R2f)), 1-fcdf(F,2,100)
F(2,100) = 0.4435,  significance = 0.6431

 save('C:\work\PLExp2\PLE_exponent2.mat')



*******************************************************************************
*****  28 Nov 2004
*****  Plot the data by difficulty level, with 5-parameter regression curves.
*****  See Plexp1/dprime_fits.txt for analogous fits with feedback.
*****
% T=[1:36]';Ts=[Inf,1:8,1:8,1:8,1:8,1:3]';                  % [36x1]
 T1=[1, 2:.5:9, 10:.5:17, 18:.5:25, 26:.5:33, 34:.5:36]';  % [66x1]
 Ts1=[Inf, 1:.5:8, 1:.5:8, 1:.5:8, 1:.5:8, 1:.5:3]';       % [66x1]
 m=reshape(PLE_exponent2(par,[T1,Ts1]),66,3);

 plot(T,edpr(:,1),'b^',T,edpr(:,2),'ro',T,edpr(:,3),'gs',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:33,m(47:61,:),'k-',34:.5:36,m(62:66,:),'k-',T1(2:15:62),m(2:15:62,:),'kx');
 legend('Contrast 0.245','Contrast 0.160','Contrast 0.106','Switch cost fits',4);
 axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);xlabel('Block number  (300 trials/block, 4 blocks/day)');ylabel('d''');
% Set page size=[1 3.25 6.5 4.5]; Line width=1; Open symbols.


*******************************************************************************
***  08 March 2005
***  Combine the two d' figures into a single plot.

 cd('C:\work\PLExp2');load estats18
 edpr = mean(cat(3,estats18(:).dprime),3); edpr(37:48,:)=[];size(edpr)

 T=[1:36]';Ts=[Inf,1:8,1:8,1:8,1:8,1:3]';                  % [36x1]
 T1=[1, 2:.5:9, 10:.5:17, 18:.5:25, 26:.5:33, 34:.5:36]';  % [66x1]
 Ts1=[Inf, 1:.5:8, 1:.5:8, 1:.5:8, 1:.5:8, 1:.5:3]';       % [66x1]
 par = [2.5819  1.8840  1.1400  0.4296  0.1496] ;   % [D1 D2 D3 g s]
 m=reshape(PLE_exponent2(par,[T1,Ts1]),66,3);

% Load the corresponding feedback data from C:\work\PLExp1\dprime_fits.mat
% Change all feedback-related variable names to end on '_f'.

 T1_f=[1, 2:.5:9, 10:.5:17, 18:.5:25, 26:.5:31, 32]';  % [58x1]
 Ts1_f=[Inf, 1:.5:8, 1:.5:8, 1:.5:8, 1:.5:6, 1]';      % [58x1]
% par7_f = [2.4145  1.8892  1.1974  0.4660  10.1551  0.1837  1.1784];  % The 7-parameter fit reported in Table 1 of the Psych Review paper
 cd('C:\work\PLExp1');m_f=reshape(PLE_exponent(par7_f,[T1_f,Ts1_f]),58,3);cd('C:\work\PLExp2');

 whos
  Name         Size         Bytes  Class
-------------------------------------------------
  T           36x1            288  double array
  T1          66x1            528  double array
  Ts          36x1            288  double array
  Ts1         66x1            528  double array
  edpr        36x3            864  double array
  m           66x3           1584  double array
  par          1x5             40  double array

  T_f         32x1            256  double array
  T1_f        58x1            464  double array
  Ts_f        32x1            256  double array
  Ts1_f       58x1            464  double array
  edpr_f      32x3            768  double array
  m_f         58x3           1392  double array
  par7_f       1x7             56  double array
-------------------------------------------------
Grand total is 972 elements using 7776 bytes

 save('C:\work\PLPaper2\fig\edprimes.mat')


%%  No-feedback empirical d' curves -- top panel of Figure 2 in PLFeedback.doc
 subplot(2,1,1);plot(T,edpr(:,1),'b^',T,edpr(:,2),'ro',T,edpr(:,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:33,m(47:61,:),'k-',34:.5:36,m(62:66,:),'k-',T1(2:15:62),m(2:15:62,:),'kx');
 axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);xlabel('Block number  (300 trials/block, 4 blocks/day)');ylabel('d''');title('No feedback  (N=18)');

%%  Feedback empirical d' curves -- bottom panel of Figure 2 in PLFeedback.doc
%%  Re-print of Figure 4 in the Psych Review paper.
 subplot(2,1,2);plot(T_f,edpr_f(:,1),'b^',T_f,edpr_f(:,2),'ro',T_f,edpr_f(:,3),'bs',1,m_f(1,:),'kx-',2:.5:9,m_f(2:16,:),'k-',10:.5:17,m_f(17:31,:),'k-',18:.5:25,m_f(32:46,:),'k-',26:.5:31,m_f(47:57,:),'k-',T1_f([2 17 32 47 58]),m_f([2 17 32 47 58],:),'kx');
 axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);xlabel('Block number  (300 trials/block, 4 blocks/day)');ylabel('d''');title('Feedback on error  (N=13)');

% Set page size=[1 1 6.5 8 ]; % [Left Bottom Width Height]


*******************************************************************************
***  08 March 2005
***  z-Probability plots, no-feedback.

 cd('C:\work\PLExp2');load estats18
 zPcongr18=cat(3,estats18(:).zP_congr);zPcongr18(37:48,:,:)=[]; % [36x3x18]
 zPincon18=cat(3,estats18(:).zP_incon);zPincon18(37:48,:,:)=[]; % [36x3x18]
 zPcongr=mean(zPcongr18,3) ; zPincon=mean(zPincon18,3);         % [36x3]
% save C:\work\PLPaper2\fig\zPcongr18.mat
 foo=[mean(zPincon) ; mean(zPcongr)];[foo,mean(foo,2);mean(foo),mean(mean(foo))]

         Cp=.23    Cp=.15    Cp=.10    Total
----------------------------------------------
incon    0.9630    0.2839   -0.4317    0.2717
congr    1.2077    1.2998    1.3902    1.2993
----------------------------------------------
d'/2     1.0853    0.7919    0.4793    0.7855


%% Plot Incongruent on top
 z=zPincon';subplot(2,1,1);hold on;
 idx=[2:9]';h=plot(idx,z(1,idx),'b^-',idx,z(2,idx),'ro-',idx,z(3,idx),'gs-');set(h,'LineWidth',1);
 idx=[10:17]';h=plot(idx,z(1,idx),'b^-',idx,z(2,idx),'ro-',idx,z(3,idx),'gs-');set(h,'LineWidth',1);
 idx=[18:25]';h=plot(idx,z(1,idx),'b^-',idx,z(2,idx),'ro-',idx,z(3,idx),'gs-');set(h,'LineWidth',1);
 idx=[26:33]';h=plot(idx,z(1,idx),'b^-',idx,z(2,idx),'ro-',idx,z(3,idx),'gs-');set(h,'LineWidth',1);
 idx=[34:36]';h=plot(idx,z(1,idx),'b^-',idx,z(2,idx),'ro-',idx,z(3,idx),'gs-');set(h,'LineWidth',1);
 idx=[1]';plot(idx,z(1,idx),'b^',idx,z(2,idx),'ro',idx,z(3,idx),'gs');set(h,'LineWidth',1);
 title('Incongruent stimuli');axis([0 37 -1 2]);h=refline(0,0);set(h,'Color','k');hold off;
 set(gca,'box','on','xtick',[0:4:36],'ytick',[-1:.5:2]);ylabel('z-probability correct')

%% Plot Congruent at bottom
 z=zPcongr';subplot(2,1,2);h=plot(1:36,z(1,:),'b^',1:36,z(2,:),'ro',1:36,z(3,:),'gs');set(h,'LineWidth',1);
 title('Congruent stimuli');axis([0 37 -1 2]);h=refline(0,0);set(h,'Color','k');
 set(gca,'box','on','xtick',[0:4:36],'ytick',[-1:.5:2]);ylabel('z-probability correct')
 xlabel('Block number  (300 trials/block, 4 blocks/day)');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(zPincon) ; mean(zPcongr)],h=axes('position',[0.69 0.12 0.2 0.14]);plot(1:2,m(:,1),'b^-',1:2,m(:,2),'ro-',1:2,m(:,3),'gs-');
 axis([.5 2.5 -.75 1.75]);set(gca,'xtick',[1 2],'xticklabel','Incon|Congr','XAxisLocation','top','ytick',[-.5:.5:1.5],'yticklabel','|0||1|');

m = 0.9630    0.2839   -0.4317
    1.2077    1.2998    1.3902

%% See PLExp2\modelfits\Hebb_fits.txt for the corresponding model fit figures



********************  OLD OLD OLD OLD OLD  ***************************
%%%%%%%%%% Congruence effects   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
 a=cat(1,st(:).accZP1),describe(a)'

      -------- Congruent -------     ------- Incongruent ------
      Cp=.23    Cp=.15    Cp=.10     Cp=.23    Cp=.15    Cp=.10
---------------------------------------------------------------
 1    1.3014    1.3714    1.5024     1.5265    0.4454   -0.5937
 2    1.1074    1.0648    1.0930     0.3767   -0.0954   -0.5605
 3    1.3032    1.5993    1.7978     1.5695    0.5301   -0.6754
 4    0.8001    0.8958    0.8887     0.7368    0.3561   -0.0329
 5    1.1561    1.4181    1.5317     0.8987    0.0566   -0.7600
 6    1.0059    1.0921    1.2278     0.8155    0.3084   -0.1115
 7    0.9330    1.0347    1.1157     1.0147    0.3710   -0.3313
 8    0.9334    0.9099    0.8971     0.7362    0.4215   -0.0962
 9    0.9648    1.0316    1.1925     1.2521    0.6105   -0.1944
10    0.7714    0.8690    0.8748     0.3808    0.0204   -0.3147
11    1.2825    1.4142    1.4512     0.9481    0.4796   -0.1264
12    1.3846    1.3766    1.4050     1.0013    0.4244   -0.0742
13    1.2371    1.3490    1.4401     0.4556   -0.0536   -0.5597
14    1.4780    1.5939    1.6559     1.3601    0.4299   -0.5727
15    1.7884    1.7624    1.8604     0.3629   -0.4467   -1.2528
16    1.3637    1.4953    1.7117     1.0536    0.2331   -0.6036
17    1.1100    1.1973    1.3362     1.2257    0.5233   -0.3164
18    1.8176    1.9213    2.0422     1.6185    0.4957   -0.5938
---------------------------------------------------------------
mean  1.2077    1.2998    1.3902    0.9630    0.2839   -0.4317
std   0.2962    0.3075    0.3445    0.4106    0.2789    0.3129

min   0.7714    0.8690    0.8748    0.3629   -0.4467   -1.2528
Q25   0.9648    1.0347    1.1157    0.7362    0.0566   -0.5938
med   1.1966    1.3602    1.4226    0.9747    0.3963   -0.4455
Q75   1.3637    1.4953    1.6559    1.2521    0.4796   -0.1264
max   1.8176    1.9213    2.0422    1.6185    0.6105   -0.0329
---------------------------------------------------------------
      Cp=.23    Cp=.15    Cp=.10     Cp=.23    Cp=.15    Cp=.10
      -------- Congruent -------     ------- Incongruent ------


%%%%%%%%%%%% Sequential effects  %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
 a=[cat(1,st(:).seq_eff),cat(1,st(:).prSprR)],describe(a)'

     -- Main effects --   ----- Conditional probabilities -----
      S(t-1)    R(t-1)     sLrL      sLrR      sRrL      sRrR
-----------------------   -------------------------------------
 1   -0.0451    0.0619     0.5008    0.5677    0.4607    0.5176
 2   -0.0402    0.0766     0.5250    0.5887    0.4718    0.5613
 3   -0.0158    0.0546     0.4677    0.5308    0.4605    0.5066
 4   -0.0132   -0.0365     0.5660    0.5133    0.5367    0.5163
 5   -0.0353    0.2122     0.4376    0.6396    0.3921    0.6144
 6   -0.0246    0.1014     0.4497    0.5564    0.4304    0.5264
 7   -0.0193    0.2488     0.3928    0.6391    0.3710    0.6224
 8   -0.0123   -0.0486     0.5808    0.5400    0.5763    0.5198
 9   -0.0261    0.0146     0.5154    0.5180    0.4772    0.5039
10   -0.0384    0.2370     0.3827    0.6049    0.3295    0.5813
11   -0.0513    0.1216     0.4712    0.5788    0.4059    0.5414
12   -0.0352    0.0211     0.5406    0.5424    0.4860    0.5265
13   -0.0429    0.1282     0.4628    0.5868    0.4157    0.5481
14   -0.0277   -0.0005     0.5763    0.5817    0.5545    0.5481
15    0.0012    0.0209     0.5015    0.5301    0.5104    0.5236
16   -0.0657    0.1290     0.4652    0.6003    0.4056    0.5284
17   -0.0470    0.1353     0.4644    0.6164    0.4342    0.5527
18   -0.0872    0.1274     0.4927    0.6218    0.4072    0.5328
-----------------------   -------------------------------------
mean -0.0348    0.0892     0.4885    0.5754    0.4514    0.5429
std   0.0209    0.0875     0.0565    0.0404    0.0653    0.0337

min  -0.0872   -0.0486     0.3827    0.5133    0.3295    0.5039
Q25  -0.0451    0.0209     0.4628    0.5400    0.4059    0.5198
med  -0.0353    0.0890     0.4819    0.5802    0.4473    0.5306
Q75  -0.0193    0.1290     0.5250    0.6049    0.4860    0.5527
max   0.0012    0.2488     0.5808    0.6396    0.5763    0.6224
-----------------------   -------------------------------------
      S(t-1)    R(t-1)     sLrL      sLrR      sRrL      sRrR
     -- Main effects --   ----- Conditional probabilities -----
