*******  PLEXP2/MODELFITS/HEBB_FITS.TXT  -- March 31, 2004   ***********
***
***  Fit PLModel1 to the no-feedback PLExp2 data.
***  Compare with PLMODEL1/SIMUL/PARAMSEARCH.HEBB.TXT
***

%%%%%  SCHED36.MAT is the PLExp2 analog to SCHED.MAT
%
 a=1;a8=ones(1,8);a3=ones(1,3); b=2;b8=b.*a8;b3=b8(1:3); T=300;
 s1=[a,b8,a8,b8,a8,b3]'; s2=[b,a8,b8,a8,b8,a3]'; N=length(s1);  %N=36
 sched36_1=expand_block_spec(num2cell([zeros(N,2),s1,repmat(T,N,1)],2));
 sched36_2=expand_block_spec(num2cell([zeros(N,2),s2,repmat(T,N,1)],2));
 bnd36=[1:8:33,36].*300
bnd36 = [300 2700 5100 7500 9900 10800]

% saved as PLExp2/modelfits/sched36.mat


%%%%%%  HUMAN_ZCONGR for PLExp2, phase 1, 36 blocks per subject
%
%  HUMAN_ZCONGR is a 6-by-36 matrix of z-transformed probabilities
%  of correct response for each stimulus type across blocks.
%  The rows (or "conditions") in ZCONGR have the same indices (or IDs)
%  as in PLM_SUMSQFIT2 and PLM_PBLOCK1. The first-dimension index is:
%   1: congruent Gabor, contrast 0.23
%   2: incongruent Gabor, contrast 0.23
%   3: congruent Gabor, contrast 0.15
%   4: incongruent Gabor, contrast 0.15
%   5: congruent Gabor, contrast 0.10
%   6: incongruent Gabor, contrast 0.10
%
 load C:\work\PLExp2\modelfits\estats18 ; estats18(1)
ans =    sbj: 1
       group: 1
    zP_congr: [48x3 double]
    zP_incon: [48x3 double]
      dprime: [48x3 double]
      accZP1: [1.3014 1.3714 1.5024 1.5265 0.4454 -0.5937]
      accZP2: [1.2798 1.2915 1.5034 1.2829 0.3310 -0.5384]
      accZP3: [2.0587 2.1340 2.0743 0.8710 0.0242 -0.9389]
    RT_congr: [48x3 double]
    RT_incon: [48x3 double]
          RT: [48x3 double]
      accRT1: [565.8611 568.8136 569.0420 625.5889 662.6574 632.6160]
      accRT2: [570.9417 562.7056 553.8722 604.3944 633.0833 609.3417]
      accRT3: [462.2333 456.1167 467.5000 571.8556 573.2667 509.2667]
    seq_data: [16x10 double]
      prSprR: [0.5008 0.5677 0.4607 0.5176]
      deltaP: [-0.0109 0.0560 -0.0510 0.0058]
     seq_eff: [-0.0451 0.0619]

 zPc=cat(3,estats18(:).zP_congr);zPi=cat(3,estats18(:).zP_incon);       % [48x3x18]
 zPc=mean(zPc,3);zPi=mean(zPi,3);plot(1:48,zPc,'b.-',1:48,zPi,'r.-');   % [48x3] averaged across Ss
 zPc(37:48,:)=[];zPi(37:48,:)=[];                                       % [36x3], ignoring phase 2 and phase 3
 Z=zeros(36,6);Z(:,[1 3 5])=zPc;Z(:,[2 4 6])=zPi;                       % rearrange
 global HUMAN_ZCONGR;HUMAN_ZCONGR=Z';clear Z zPc zPi;describe(HUMAN_ZCONGR');

    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.208    0.190      0.79    1.08    1.20    1.37    1.50
   0.963    0.250      0.25    0.87    1.00    1.12    1.33
   1.300    0.200      0.66    1.17    1.28    1.45    1.63
   0.284    0.143     -0.20    0.23    0.28    0.38    0.54
   1.390    0.188      0.85    1.26    1.44    1.52    1.70
  -0.432    0.106     -0.68   -0.52   -0.45   -0.35   -0.26
------------------------------------------------------------
   0.785    0.179      0.28    0.68    0.79    0.91    1.07

 save('C:\work\PLExp2\modelfits\HUMAN_ZCONGR.mat','HUMAN_ZCONGR')


%%%%%%  Generate 10+10 stimulus sequences according to STIM36.MAT and using the old PLM_CACHE
%
 load('C:\work\PLExp2\modelfits\sched36.mat');
 load('C:\work\PLModel1\PLM_CACHE'); 
 descr=zeros(10800,6,20);   % descr(:,6,:)==0 indicating lack of feedback
 for k=1:10 d=block_from_cache(sched36_1,PLM_CACHE);descr(:,1:5,k)=cat(1,d{:});end;for k=11:20 d=block_from_cache(sched36_2,PLM_CACHE);descr(:,1:5,k)=cat(1,d{:});end;foo=squeeze(descr(1,:,:));openvar foo

%% Pack in MDATA suitable for passing as a single argument to PLM_HEBB1
 for k=1:20 Mdata(k).stim=descr(:,:,k);Mdata(k).chev=make_chev(10800);end;


%%%%%%  Use the model parameters optimized for the feedback data (PLExp1) as initial approximation for the no-feedback parameter search.
%
 Mparams = PLM_params
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.1950
      W_minmax: [-1 1]
        W_init: 0.1700
       W0_seed: [35x1 double]
    learn_rate: 0.0015
    runav_rate: 0.0200
     criterion: 2.2000
    fdbk_fract: 1
      blk_size: 300


%%%%%  Test-run to make sure the program runs
%
 md=Mdata(1);md.params=Mparams
md =  stim: [10800x6 double]
      chev: [10800x2 double]
    params: [1x1 struct]

 Z=PLM_Zcongr1(md);plot(Z,'.-');   % [6x36]
 PLM_sumsqfit2(Z)
ans = 82.3720

 ssq=zeros(20,1);for k=1:20 md=Mdata(k);md.params=Mparams;Z=PLM_Zcongr1(md);ssq(k)=PLM_sumsqfit2(Z);end;describe(ssq)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
  71.172    8.205     55.28   65.05   71.13   76.21   86.08

 reshape(ssq,5,4)
   --  Schedule A --   --  Schedule B --
   82.3720   80.1614   67.2745   64.4364
   83.3742   72.0646   68.2646   62.6160
   55.2831   59.0270   77.2790   73.0840
   74.2242   75.1421   70.2002   86.0824
   74.0960   65.3927   64.7053   68.3665


%%%%%%%%%%%%  SEARCH PARAMETERS are the same as for the feedback fits
%
 Sparams=PLM_search_params;Sparams.params=Mparams
Sparams = 
    lsfun_name: 'PLM_paramsearch'
    sumsq_name: 'PLM_sumsqfit2'
    model_name: 'PLM_Zcongr1'
        params: [1x1 struct]
     p2v_templ: {1x5 cell}
     v2p_templ: {1x5 cell}
        bounds: [2x5 double]
         optns: [1x1 struct]
       verbose: 0

 Sparams.p2v_templ'
ans = 
    'VAL = PARAMS.learn_rate *100 ;'
    'VAL = PARAMS.out_noise ;'
    'VAL = PARAMS.W_init ;'
    'VAL = PARAMS.rep_gain ;'
    'VAL = PARAMS.criterion ;'

 Sparams.bounds
ans =
    0.0500    0.1000    0.0500    0.3000    1.0000
    1.0000    0.5000    0.4000    1.0000    5.0000

 x0=params2vector(Mparams,Sparams.p2v_templ)
x0= 0.1500    0.1950    0.1700    0.8000    2.2000

% Test-run on 1 stimulus set
 tic;md=Sparams;md.model_data=Mdata(1);foo=PLM_paramsearch(x0,md);toc;foo
foo = 82.3720    % elapsed_time = 3.3950

% Test-run on 20 stimulus sets
 tic;md=Sparams;md.model_data=Mdata;foo=PLM_paramsearch(x0,md);toc;foo
 tic;md=Sparams;md.model_data=Mdata;foo=PLM_paramsearch(x0,md);toc;foo
foo = 71.1723    % elapsed_time = 61.7990


%%%%%  Now try the search itself  May 31, 2004
%
 tic;[Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams),toc

Sres =                       % elapsed_time = 1 hr 45 min, 124 function evaluations
  opt_par: [1x1 struct]
      ssq: 14.7507
    opt_X: [0.0500 0.2125 0.3999 0.4446 1.3316]
      det: [1x1 struct]

 Mparams=Sres.opt_par
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.4446
     rep_noise: 0.1000
     out_noise: 0.2125
      W_minmax: [-1 1]
        W_init: 0.3999
       W0_seed: [35x1 double]
    learn_rate: 5.0000e-004
    runav_rate: 0.0200
     criterion: 1.3316
    fdbk_fract: 1
      blk_size: 300


%%%%%%  Test the new parameters on 50+50 fresh stimulus sequences
%
 descr=zeros(10800,6,100);for k=1:50 d=block_from_cache(sched36_1,PLM_CACHE);descr(:,1:5,k)=cat(1,d{:});end;for k=51:100 d=block_from_cache(sched36_2,PLM_CACHE);descr(:,1:5,k)=cat(1,d{:});end;foo=squeeze(descr(1,:,:));openvar foo
 tic;Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Sres.opt_par);end;toc
elapsed_time = 207 sec on a Dell Optiplex GX400

 tic;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);toc
elapsed_time = 53 sec

 foo=reshape(sst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean
----------------------------------------
    1.2442    1.3091    1.3217    1.2917    % congruent
    1.0623    0.1191   -0.5242    0.2191    % incongruent
----------------------------------------
    2.3065    1.4282    0.7974    1.5107    % dprime

 nfsst = sst;plot_nfsst;sst05=sst;
 describe(nfsst.zP_congr)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.244    0.081      1.11    1.19    1.23    1.32    1.40
   1.309    0.070      1.20    1.25    1.30    1.35    1.46
   1.322    0.061      1.22    1.29    1.32    1.35    1.51
------------------------------------------------------------
   1.292    0.071      1.18    1.24    1.28    1.34    1.46

 describe(nfsst.zP_incon)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.062    0.118      0.80    0.97    1.06    1.15    1.25
   0.119    0.071     -0.02    0.06    0.12    0.17    0.26
  -0.524    0.048     -0.66   -0.55   -0.52   -0.49   -0.45
------------------------------------------------------------
   0.219    0.079      0.04    0.16    0.22    0.28    0.35

 describe(nfsst.dprime)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   2.307    0.111      1.99    2.22    2.32    2.38    2.49
   1.428    0.064      1.31    1.37    1.44    1.47    1.55
   0.797    0.039      0.71    0.77    0.79    0.83    0.87
------------------------------------------------------------
   1.511    0.071      1.34    1.46    1.52    1.56    1.64

 foo=Wh(:,1:10:end);
 for k=1:5 subplot(2,3,k);plot(foo([1:7]+7*k-7,:)','b');axis([0 1090 -.5 +.5]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.5:.1:+.5]);grid on;end
 subplot(2,3,6);plot(foo(36,:)','k');axis([0 1090 -.7 +.7]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.6:.2:+.6]);grid on

%%%  These fits are much better than the fits under the feedback-derived Mparams.
%%%  At least now there is no (or much weeker) unstable oscillation.
%%%  But that's in large part because the learning rate is so low (0.0005) that esentially there is no learning.
%%%  So, restrict the learning rate parameter and try again.

 Mparams = Sres.opt_par;Mparams.learn_rate=0.0015;  % same as with feedback
 Sparams.params=Mparams;Sparams.p2v_templ(1)=[];Sparams.v2p_templ(1)=[];Sparams.bounds(:,1)=[]
Sparams = 
    lsfun_name: 'PLM_paramsearch'
    sumsq_name: 'PLM_sumsqfit2'
    model_name: 'PLM_Zcongr1'
        params: [1x1 struct]
     p2v_templ: {1x4 cell}
     v2p_templ: {1x4 cell}
        bounds: [2x4 double]
         optns: [1x1 struct]
       verbose: 0

 Sparams.p2v_templ'
    'VAL = PARAMS.out_noise ;'
    'VAL = PARAMS.W_init ;'
    'VAL = PARAMS.rep_gain ;'
    'VAL = PARAMS.criterion ;'

 Sparams.bounds(2,1)=0.5;Sparams.bounds(2,2)=0.5;Sparams.bounds(2,4)=3;Sparams.bounds
    0.1000    0.0500    0.3000    1.0000
    0.5000    0.5000    1.0000    3.0000

 x0=params2vector(Sparams.params,Sparams.p2v_templ)
x0= 0.2125    0.3999    0.4446    1.3316

 [Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams);
Sres = 
    opt_par: [1x1 struct]
        ssq: 18.1222
      opt_X: [0.2261 0.4974 0.3867 1.6463]
        det: [1x1 struct]

 tic;Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Sres.opt_par);end;toc
elapsed_time = 250 sec on a Dell Optiplex GX400

 tic;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));sst=PLM_summary(st);toc
elapsed_time = 79 sec

 foo=reshape(sst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean
----------------------------------------
    1.2207    1.2638    1.2754    1.2533
    1.1665    0.2065   -0.4429    0.3100
----------------------------------------
    2.3872    1.4704    0.8325    1.5634

 nfsst = sst;plot_nfsst;sst15=sst;
 describe(nfsst.zP_congr)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.221    0.110      1.04    1.12    1.22    1.32    1.40
   1.264    0.083      1.12    1.19    1.26    1.33    1.47
   1.275    0.052      1.19    1.23    1.27    1.31    1.38
------------------------------------------------------------
   1.253    0.082      1.11    1.18    1.25    1.32    1.42

 describe(nfsst.zP_incon)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.166    0.158      0.91    1.03    1.14    1.32    1.46
   0.207    0.107      0.03    0.13    0.17    0.31    0.42
  -0.443    0.075     -0.58   -0.50   -0.46   -0.38   -0.31
------------------------------------------------------------
   0.310    0.113      0.12    0.22    0.29    0.42    0.52

 describe(nfsst.dprime)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   2.387    0.094      2.15    2.34    2.41    2.46    2.53
   1.470    0.066      1.36    1.43    1.47    1.53    1.59
   0.833    0.051      0.70    0.80    0.83    0.86    0.94
------------------------------------------------------------
   1.563    0.070      1.40    1.52    1.57    1.62    1.68

 foo=Wh(:,1:10:end);for k=1:5 subplot(2,3,k);plot(foo([1:7]+7*k-7,:)','b');axis([0 1090 -.5 +.5]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.5:.1:+.5]);grid on;end
 subplot(2,3,6);plot(foo(36,:)','k');axis([0 1090 -.7 +.7]);grid on



%%%%%  Try again, with initial W_init=0.20;
 Mparams=Sres.opt_par;Mparams.W_init=0.20,Sparams.params=Mparams;
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.3867
     rep_noise: 0.1000
     out_noise: 0.2261
      W_minmax: [-1 1]
        W_init: 0.2000
       W0_seed: [35x1 double]
    learn_rate: 0.0015
    runav_rate: 0.0200
     criterion: 1.6463
    fdbk_fract: 1
      blk_size: 300

 [Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams);
Sres =                   % MaxiFunEvals (120) exceeded.
    opt_par: [1x1 struct]
        ssq: 24.6730
      opt_X: [0.2263 0.3322 0.5400 1.8758]
        det: [1x1 struct]



*****************************************************************
*****************************************************************
*****************************************************************
*** April 2, 2004
***
*** Change the learning rule so that the self-teaching signal
*** has a mean of zero. Line 172 in PLM_HEBB1.M now reads:
***   dW = (learn_rate*(curr_fdbk-runav)) .* rep' ;
*** The old (PLExp1) version took advantage that the mean of the EXTERNAL
*** feedback is always zero, and didn't subtract anything from CURR_FDBK.
***
*****************************************************************
*****************************************************************
*****************************************************************

%%%  Now, re-do the fits with the mean-centered learning rule:
%%
Mparams =                   % Initial approximation informed by earlier searches
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.5889      % that's too low; fix at .80 to get the reversal for congruent stimuli
     rep_noise: 0.1000
     out_noise: 0.1810
      W_minmax: [-1 1]
        W_init: 0.2599
       W0_seed: [35x1 double]
    learn_rate: 0.0025      % <-- fixed
    runav_rate: 0.0200
     criterion: 1           % <-- at low bound
    fdbk_fract: 1
      blk_size: 300

 Mparams.rep_gain=0.8;Sparams.params=Mparams
Sparams = 
    lsfun_name: 'PLM_paramsearch'
    sumsq_name: 'PLM_sumsqfit2'
    model_name: 'PLM_Zcongr1'
        params: [1x1 struct]
     p2v_templ: {1x4 cell}
     v2p_templ: {1x4 cell}
        bounds: [2x4 double]
         optns: [1x1 struct]
       verbose: 0

 Sparams.p2v_templ(3)=[];Sparams.v2p_templ(3)=[];Sparams.bounds(:,3)=[];
 Sparams.p2v_templ'
    'VAL = PARAMS.out_noise ;'
    'VAL = PARAMS.W_init ;'
    'VAL = PARAMS.criterion ;'

 Sparams.bounds(1,3)=0.1;Sparams.bounds
    0.1000    0.0500    0.1000
    0.5000    0.4000    3.0000

 x0=params2vector(Sparams.params,Sparams.p2v_templ)
x0= 0.1810    0.2599    1.0000

 [Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams);
Sres = 
    opt_par: [1x1 struct]
        ssq: 13.6502
      opt_X: [0.1782 0.1996 1.0740]
        det: [1x1 struct]

 Sparams.params=Sres.opt_par;Sparams.params.W_init=0.18;
 [Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams);
Sres = 
    opt_par: [1x1 struct]
        ssq: 12.7072
      opt_X: [0.1768 0.1962 0.9724]
        det: [1x1 struct]

 Mparams=Sres.opt_par
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.1768
      W_minmax: [-1 1]
        W_init: 0.1962
       W0_seed: [35x1 double]
    learn_rate: 0.0025
    runav_rate: 0.0200
     criterion: 0.9724
    fdbk_fract: 1
      blk_size: 300

 tic;Nruns=20;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;toc
elapsed_time =  43 sec on a Dell OptiPlex GX400

 tic;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));nfsst=PLM_summary(st);toc
elapsed_time = 10.5950

 plot_nfsst ; sst25=nfsst;Sres25=Sres;    % since learn_rate=0.0025
 foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % learn_rate=0.0025
----------------------------------------
    1.1920    1.2767    1.3214    1.2634   % see 100-run-based estimates below
    1.0967    0.2310   -0.4488    0.2929
----------------------------------------
    2.2887    1.5077    0.8726    1.5563

 describe(nfsst.zP_congr)
    Mean  Std.dev       Min     Q25  Median     Q75     Max   % learn_rate=0.0025
------------------------------------------------------------
   1.192    0.221      0.88    1.03    1.16    1.36    1.67
   1.277    0.206      0.90    1.14    1.25    1.44    1.71
   1.321    0.145      1.04    1.24    1.32    1.39    1.72
------------------------------------------------------------
   1.263    0.191      0.94    1.13    1.24    1.40    1.70

 describe(nfsst.zP_incon)
    Mean  Std.dev       Min     Q25  Median     Q75     Max   % learn_rate=0.0025
------------------------------------------------------------
   1.097    0.309      0.36    0.93    1.11    1.33    1.63
   0.231    0.244     -0.29    0.09    0.25    0.43    0.61
  -0.449    0.193     -0.93   -0.55   -0.44   -0.29   -0.19
------------------------------------------------------------
   0.293    0.249     -0.29    0.15    0.31    0.49    0.69

 describe(nfsst.dprime)
    Mean  Std.dev       Min     Q25  Median     Q75     Max   % learn_rate=0.0025
------------------------------------------------------------
   2.289    0.253      1.61    2.17    2.31    2.48    2.66
   1.508    0.180      0.97    1.41    1.50    1.64    1.94
   0.873    0.127      0.62    0.81    0.88    0.97    1.08
------------------------------------------------------------
   1.556    0.187      1.06    1.46    1.56    1.70    1.89

 foo=Wh(:,1:10:end);for k=1:5 subplot(2,3,k);plot(foo([1:7]+7*k-7,:)','b');axis([0 1090 -.5 +.5]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.5:.1:+.5]);grid on;end;subplot(2,3,6);plot(foo(36,:)','k');axis([0 1090 -.7 +.7]);grid on


%%%%%%  LEARN_RATE=0.0020
 Sparams20=Sparams;Sparams20.params.learn_rate=0.0020;
 [Sres20.opt_par,Sres20.ssq,Sres20.opt_X,Sres20.det]=paramsearch(Mdata,Sparams20);disp(Sres20);
    opt_par: [1x1 struct]
        ssq: 10.3240
      opt_X: [0.1602 0.1532 0.8323]
        det: [1x1 struct]


%%%%%%  LEARN_RATE=0.0025 -- the same as in the feedback condition (PLExp1)
 Sparams15=Sparams;Sparams15.params.learn_rate=0.0015;
[Sres15.opt_par,Sres15.ssq,Sres15.opt_X,Sres15.det]=paramsearch(Mdata,Sparams15);disp(Sres15);
    opt_par: [1x1 struct]
        ssq: 10.1589
      opt_X: [0.1661 0.1736 0.8528]
        det: [1x1 struct]


%%%%%%  Test the new parameters on 50+50 fresh stimulus sequences
%
 descr=zeros(10800,6,100);for k=1:50 d=block_from_cache(sched36_1,PLM_CACHE);descr(:,1:5,k)=cat(1,d{:});end;for k=51:100 d=block_from_cache(sched36_2,PLM_CACHE);descr(:,1:5,k)=cat(1,d{:});end;foo=squeeze(descr(1,:,:));openvar foo
 tic;Mparams=Sres15.opt_par;Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));nfsst=PLM_summary(st);toc;figure(1);plot_nfsst
elapsed_time = 319 sec on a Dell Optiplex GX400

 sst15=nfsst;foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % learn_rate=0.0015
----------------------------------------
    1.1809    1.2690    1.3225    1.2575
    0.9954    0.1535   -0.4771    0.2239
----------------------------------------
    2.1762    1.4225    0.8454    1.4814

 foo=Wh(:,1:10:end);figure(2);for k=1:5 subplot(2,3,k);plot(foo([1:7]+7*k-7,:)','b');axis([0 1090 -.5 +.5]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.5:.1:+.5]);grid on;end;subplot(2,3,6);plot(foo(36,:)');axis([0 1090 -.7 +.7]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.6:.2:+.6]);grid on

%%%%%%%%%%%%%%%%
 tic;Mparams=Sres20.opt_par;Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));nfsst=PLM_summary(st);toc;figure(1);plot_nfsst
elapsed_time = 316 sec on a Dell Optiplex GX400

 sst20=nfsst;foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % learn_rate=0.0020
----------------------------------------
    1.2032    1.2824    1.3208    1.2688
    1.0477    0.2065   -0.4513    0.2677
----------------------------------------
    2.2509    1.4889    0.8696    1.5365

%%%%%%%%%%%%%%%%
 tic;Mparams=Sres25.opt_par;Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));nfsst=PLM_summary(st);toc;figure(1);plot_nfsst
elapsed_time = 316 sec on a Dell Optiplex GX400

 sst25=nfsst;foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % learn_rate=0.0025
----------------------------------------
    1.2043    1.2957    1.3327    1.2776
    1.1154    0.2303   -0.4378    0.3026
----------------------------------------
    2.3197    1.5260    0.8948    1.5802


%%%%%%%%%%%%%%%%   Feedback-derived parameters, see C:\work\PLModel1\simul\paramsearch.Hebb.txt
 tic;Mparams=PLM_params;Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));nfsst=PLM_summary(st);toc;figure(1);plot_nfsst
elapsed_time = 307 sec on a Dell Optiplex GX400

 sst_fdbk=nfsst;foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % learn_rate=0.0015, feedback-derived parameters
----------------------------------------
    0.6195    0.6881    0.7221    0.6765
    0.9190    0.2999   -0.1527    0.3554
----------------------------------------
    1.5385    0.9880    0.5694    1.0320


***************************************************************************
****
****  April 6, 2004  --  Final no-feedback fits (hopefully...)
****

% Fix the W_init to the feedback-derived value = 0.17
% This leaves only two free parameters: out_noise and criterion
 Sparams=Sres15.Sparams;Sparams.params.W_init=0.17;Sparams.params.rep_gain=0.80;
 Sparams.params.out_noise=0.17;Sparams.params.criterion=0.9;Sparams.params
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.1700
      W_minmax: [-1 1]
        W_init: 0.1700
       W0_seed: [35x1 double]
    learn_rate: 0.0015
    runav_rate: 0.0200
     criterion: 0.9000
    fdbk_fract: 1
      blk_size: 300

 Sparams.p2v_templ(2)=[];Sparams.v2p_templ(2)=[];Sparams.bounds(:,2)=[]
Sparams = 
    lsfun_name: 'PLM_paramsearch'
    sumsq_name: 'PLM_sumsqfit2'
    model_name: 'PLM_Zcongr1'
        params: [1x1 struct]
     p2v_templ: {'VAL = PARAMS.out_noise ;'  'VAL = PARAMS.criterion ;'}
     v2p_templ: {'PARAMS.out_noise = VAL ;'  'PARAMS.criterion = VAL ;'}
        bounds: [2x2 double]
         optns: [1x1 struct]
       verbose: 0

 Sparams.bounds
    0.1000    0.1000
    0.5000    3.0000

 [Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams);disp(Sres);
    opt_par: [1x1 struct]         
        ssq: 9.0005               % pass 1
      opt_X: [0.1580 0.8721]
        det: [1x1 struct]
          exitflag: 1
           N_evals: 29
           N_iters: 4
              time: 1.4293e+003

 Sparams.params=Sres.opt_par;[Sres.opt_par,Sres.ssq,Sres.opt_X,Sres.det]=paramsearch(Mdata,Sparams);disp(Sres);
    opt_par: [1x1 struct]
        ssq: 8.9745               % pass 2, final
      opt_X: [0.1578 0.8679]
        det: [1x1 struct]
          exitflag: 1
           N_evals: 39
           N_iters: 3
              time: 2.8556e+003

%% The summed squared error is SSQ=8.9745=9*0.997 for the no-feedback condition.
%% For comparison, in the feedback condition it is SSQ=7.7384=8*0.9673. See C:\work\PLModel1\simul\paramsearch.Hebb.txt.

 Sres_final=Sres;clear Sres;save Sres Sres_final Sres15 Sres20 Sres25
 Mparams=Sres_final.opt_par
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.1578
      W_minmax: [-1 1]
        W_init: 0.1700
       W0_seed: [35x1 double]
    learn_rate: 0.0015
    runav_rate: 0.0200
     criterion: 0.8679
    fdbk_fract: 1
      blk_size: 300

 H=Hebb_stats(Mparams);print_Hebb_stats(H);     % cf. C:\work\PLPaper1\scripts\HomogVariance.txt

Condition   dprime   bias  : mnYcon  mnYinc  :   sdY   SDcon  SDinc 
-------------------------------------------------------------------
Wsame
  Cp=.23    2.002  -0.054  :  0.108  -0.215  :  0.161  0.026  0.035
  Cp=.15    1.444   0.003  :  0.119  -0.114  :  0.162  0.026  0.042
  Cp=.10    0.885   0.055  :  0.126  -0.017  :  0.161  0.026  0.041
  Mean  +   1.444   0.001  +  0.118  -0.115  +  0.161  0.026  0.039
  Diff     -1.118          :  0.018   0.198  :

Wother
  Cp=.23    1.563   0.265  :  0.392   0.139  :  0.162  0.025  0.045
  Cp=.15    1.038   0.303  :  0.387   0.219  :  0.162  0.028  0.045
  Cp=.10    0.595   0.331  :  0.379   0.283  :  0.162  0.030  0.041
  Mean  +   1.065   0.300  +  0.386   0.213  +  0.162  0.028  0.044
  Diff     -0.969          : -0.012   0.144  :

Wmix
  Cp=.23    1.783   0.106  :  0.250  -0.038  :  0.161  0.025  0.040
  Cp=.15    1.242   0.153  :  0.253   0.052  :  0.162  0.026  0.043
  Cp=.10    0.740   0.193  :  0.253   0.133  :  0.161  0.027  0.040
  Mean  +   1.255   0.151  +  0.252   0.049  +  0.162  0.026  0.041
  Diff     -1.043          :  0.003   0.171  :

W0step
  Cp=.23    2.377   0.334  :  0.536   0.132  :  0.170  0.053  0.072
  Cp=.15    1.546   0.426  :  0.558   0.294  :  0.170  0.054  0.073
  Cp=.10    0.881   0.493  :  0.568   0.418  :  0.170  0.054  0.069
  Mean  +   1.601   0.418  +  0.554   0.281  +  0.170  0.054  0.072
  Diff     -1.496          :  0.032   0.287  :

W0lin
  Cp=.23    1.778   0.264  :  0.412   0.117  :  0.166  0.047  0.056
  Cp=.15    1.118   0.336  :  0.429   0.244  :  0.166  0.048  0.057
  Cp=.10    0.618   0.387  :  0.439   0.336  :  0.166  0.047  0.055
  Mean  +   1.171   0.329  +  0.427   0.232  +  0.166  0.047  0.056
  Diff     -1.161          :  0.027   0.220  :
-------------------------------------------------------------------
Condition   dprime   bias  : mnYcon  mnYinc  :   sdY   SDcon  SDinc 


***************************************************************************
****
****  April 6, 2004  --  Test the model with the final parameters, 50+50 runs
****

 Nruns=100;o=zeros(10800,5,Nruns);for k=1:Nruns [o(:,:,k),Wh]=PLM_Hebb1(descr(:,:,k),Mparams);end;st=PLM_stats2(descr(:,:,1:Nruns),o(:,1,:),o(:,2,:),o(:,3,:));nfsst=PLM_summary(st);figure(1);plot_nfsst
 sst_final=nfsst;foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % final parameters, learn_rate=0.0015; rep_gain=0.80; W_init=0.17;out_noise=0.1578;criterion=0.8679
----------------------------------------
    1.2280    1.3143    1.3712    1.3045   % congruent stimuli, z(Pcorrect)
    1.0703    0.1880   -0.4738    0.2615   % incongruent stimuli, z(Pcorrect)
----------------------------------------
    2.2984    1.5022    0.8974    1.5660   % dprime = zP_congr + zP_incongr

 describe(nfsst.zP_congr)
    Mean  Std.dev       Min     Q25  Median     Q75     Max   % final params
------------------------------------------------------------
   1.228    0.192      0.90    1.11    1.23    1.35    1.64   % Cpeak=.23
   1.314    0.175      1.02    1.22    1.32    1.46    1.69   % Cpeak=.15
   1.371    0.161      1.11    1.29    1.38    1.49    1.72   % Cpeak=.10
------------------------------------------------------------
   1.305    0.176      1.01    1.20    1.31    1.43    1.68

 describe(nfsst.zP_incon)
    Mean  Std.dev       Min     Q25  Median     Q75     Max   % final params
------------------------------------------------------------
   1.070    0.250      0.48    0.89    1.10    1.26    1.50   % Cpeak=.23
   0.188    0.183     -0.26    0.07    0.21    0.32    0.49   % Cpeak=.15
  -0.474    0.129     -0.79   -0.56   -0.46   -0.37   -0.27   % Cpeak=.10
------------------------------------------------------------
   0.262    0.187     -0.19    0.14    0.28    0.40    0.58

 describe(nfsst.dprime)
    Mean  Std.dev       Min     Q25  Median     Q75     Max   % final params
------------------------------------------------------------
   2.298    0.274      1.71    2.12    2.32    2.52    2.72   % Cpeak=.23
   1.502    0.204      1.10    1.34    1.52    1.66    1.84   % Cpeak=.15
   0.897    0.123      0.65    0.78    0.92    1.01    1.09   % Cpeak=.10
------------------------------------------------------------
   1.566    0.201      1.15    1.41    1.59    1.73    1.89

 foo=Wh(:,1:10:end);figure(2);for k=1:5 subplot(2,3,k);plot(foo([1:7]+7*k-7,:)','b');axis([0 1090 -.5 +.5]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.5:.1:+.5]);grid on;end;subplot(2,3,6);plot(foo(36,:)');axis([0 1090 -.7 +.7]);set(gca,'xtick',[0:120:1080],'xticklabel',[],'ytick',[-.6:.2:+.6]);grid on



***************************************************************************
****
****  April 7, 2004  --  Big model test, 500+500 runs
****

 clear all ; load('C:\work\PLModel1\PLM_CACHE.mat') ; load('C:\work\PLExp2\modelfits\HUMAN_ZCONGR.mat') ; load('C:\work\PLExp2\modelfits\sched36.mat') ;

% A brute-force approach doesn't work. Workspace size gets so big that it spills into virtual memory, which takes forever... Divide-and-conquer is needed.
 tic;descr1=zeros(10800,6,500);descr2=descr1;for k=1:500 d=block_from_cache(sched36_1,PLM_CACHE);descr1(:,1:5,k)=cat(1,d{:});end;for k=1:500 d=block_from_cache(sched36_2,PLM_CACHE);descr2(:,1:5,k)=cat(1,d{:});end;toc;
elapsed_time = 58 minutes on a Dell Optiplex GX400

 Mparams=PLM_params;Mparams.out_noise=.1578;Mparams.criterion=.8679
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.1578
      W_minmax: [-1 1]
        W_init: 0.1700
       W0_seed: [35x1 double]
    learn_rate: 0.0015
    runav_rate: 0.0200
     criterion: 0.8679
    fdbk_fract: 1
      blk_size: 300

 tic;Nruns=500;o1=zeros(10800,5,Nruns);for k=1:Nruns [o1(:,:,k),Wh]=PLM_Hebb1(descr1(:,:,k),Mparams);if (mod(k,10)==5) fprintf('.');end;end;toc
 tic;st1=PLM_stats2(descr1,o1(:,1,:),o1(:,2,:),o1(:,3,:));toc     % elapsed_time = 315 sec
 save st1 st1 ; clear all ; load descr2 ;  load('C:\work\PLModel1\PLM_CACHE.mat') ;

 tic;Nruns=500;o2=zeros(10800,5,Nruns);for k=1:Nruns [o2(:,:,k),Wh]=PLM_Hebb1(descr2(:,:,k),Mparams);if (mod(k,10)==5) fprintf('.');end;end;toc     % elapsed_time = 19 min
 tic;st2=PLM_stats2(descr2,o2(:,1,:),o2(:,2,:),o2(:,3,:));toc     % elapsed_time = 392 sec
 load st1; model_stats=merge_PLM_stats(st1,st2);
model_stats = 
        group: [1000x1 double]
     zP_congr: [36x3x1000 double]
     zP_incon: [36x3x1000 double]
       dprime: [36x3x1000 double]
        accZP: [1000x6 double]
    mnA_congr: [36x3x1000 double]
    mnA_incon: [36x3x1000 double]
       Aprime: [36x3x1000 double]
         accA: [1000x6 double]
    mnI_congr: [36x3x1000 double]
    mnI_incon: [36x3x1000 double]
       Iprime: [36x3x1000 double]
         accI: [1000x6 double]

 nfsst=model_sumstats;plot_nfsst;foo=reshape(nfsst.accZP,3,2)';foo=[foo mean(foo,2)];foo=[foo;sum(foo)]
     .23       .15       .10       mean    % final parameters, 500+500 runs
----------------------------------------
    1.2233    1.3194    1.3642    1.3023
    1.0694    0.1874   -0.4775    0.2598
----------------------------------------
    2.2927    1.5068    0.8867    1.5621

 mZ=zeros(36,6);mZ(:,[1 3 5])=model_sumstats.zP_congr;mZ(:,[2 4 6])=model_sumstats.zP_incon;
 global MODEL_ZCONGR;MODEL_ZCONGR=mZ';describe(MODEL_ZCONGR');


******************************************
***
***  April 8, 2004 -- Plot the fits
***

 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


************  D-prime plot
subplot(3,1,3);hold on;idx=[1:36]';h=plot(idx,edpr(:,1),'b^',idx,edpr(:,2),'ro',idx,edpr(:,3),'gs');set(h,'MarkerSize',4);
idx={[2:9]',[10:17]',[18:25]',[26:33]',[34:36]',[1]'};for k=1:length(idx);I=idx{k};h=plot(I,mdpr(I,1),'b.-',I,mdpr(I,2),'r.-',I,mdpr(I,3),'g.-');end
hold off;title('D-prime');axis([0 37 0 3]);set(gca,'box','on','xtick',[0:4:36],'ytick',[0:.5:3]);ylabel('d''');

************  Z-incongr plot
subplot(3,1,1);hold on;Z=eZ(:,[2 4 6]);I=[1:36]';h=plot(I,Z(:,1),'b^',I,Z(:,2),'ro',I,Z(:,3),'gs');set(h,'MarkerSize',4);
Z=mZ(:,[2 4 6]);idx={[2:9]',[10:17]',[18:25]',[26:33]',[34:36]',[1]'};for k=1:length(idx);I=idx{k};h=plot(I,Z(I,1),'b.-',I,Z(I,2),'r.-',I,Z(I,3),'g.-');end
hold off;title('Incongruent stimuli');axis([0 37 -1 2.5]);set(gca,'box','on','xtick',[0:4:36],'ytick',[-1:.5:2.5]);ylabel('z-probability correct');

************  Z-congruent plot
subplot(3,1,2);hold on;Z=eZ(:,[1 3 5]);I=[1:36]';h=plot(I,Z(:,1),'b^',I,Z(:,2),'ro',I,Z(:,3),'gs');set(h,'MarkerSize',4);
Z=mZ(:,[1 3 5]);idx={[2:9]',[10:17]',[18:25]',[26:33]',[34:36]',[1]'};for k=1:length(idx);I=idx{k};h=plot(I,Z(I,1),'b.-',I,Z(I,2),'r.-',I,Z(I,3),'g.-');end
hold off;title('Congruent stimuli');axis([0 37 -1 2.5]);set(gca,'box','on','xtick',[0:4:36],'ytick',[-1:.5:2.5]);ylabel('z-probability correct');



****************************************************************************
****************************************************************************
****************************************************************************
***
*** May 2, 2005 -- Replicate the feedback simulations from the Psych Review paper
*** using the new zero-centered learning rule. (See the entry of 2 April 2004 above.)
*** This requires a change to PLM_HEBB1.M, defining a new variable TEACH_RUNAV.
*** The learning rule (line 179) now reads:
***   runav = runav + runav_rate*(sign(resp-1.5)-teach_runav) ;
***
*** When there is feedback, the weight histories under the two learning rules
*** differ by less than .001, even though teach_runav oscillates b/n -.3 and +.3
*** with a standard deviation of .094 -- comparable to that of the response runav!
***
****************************************************************************
****************************************************************************
****************************************************************************

 clear all;load('C:\work\PLModel1\PLM_CACHE.mat');load('C:\work\PLModel1\sched.mat');

% A brute-force approach doesn't work. Workspace size gets so big that it spills into virtual memory, which takes forever... Divide-and-conquer is needed.
 tic;descr1=zeros(9600,5,500);descr2=descr1;for k=1:500 d=block_from_cache(sched1,PLM_CACHE);descr1(:,:,k)=cat(1,d{:});end;for k=1:500 d=block_from_cache(sched2,PLM_CACHE);descr2(:,:,k)=cat(1,d{:});end;toc;
elapsed_time = 27 minutes on a Dell Optiplex GX400
 cd('C:\work\PLPaper2\simul');descr1=uint16(descr1);save descr1 descr1;descr2=uint16(descr2);save descr2 descr2;  % 47MB each

 clear all;load('C:\work\PLModel1\PLM_CACHE.mat');load('C:\work\PLPaper2\simul\descr1');descr1=double(descr1);
 Mparams=PLM_params
Mparams = 
      rep_size: [7 5 1]
       max_act: 0.5000
      rep_gain: 0.8000
     rep_noise: 0.1000
     out_noise: 0.1950
      W_minmax: [-1 1]
        W_init: 0.1700
       W0_seed: [35x1 double]
    learn_rate: 0.0015
    runav_rate: 0.0200
     criterion: 2.2000
    fdbk_fract: 1
      blk_size: 300

 tic;Nruns=500;o1=zeros(9600,5,Nruns);for k=1:Nruns [o1(:,:,k),Wh]=PLM_Hebb1(descr1(:,:,k),Mparams);if (mod(k,10)==5) fprintf('.');end;end;toc  % elapsed time 17 min
 cd('C:\work\PLExp2\modelfits');tic;st1=PLM_stats2(descr1,o1(:,1,:),o1(:,2,:),o1(:,3,:));toc     % elapsed_time = 6 min
 cd('C:\work\PLPaper2\simul');save st1 st1 ; 

 clear all;load('C:\work\PLModel1\PLM_CACHE.mat');load('C:\work\PLPaper2\simul\descr2');descr2=double(descr2);Mparams=PLM_params;
 tic;Nruns=500;o1=zeros(9600,5,Nruns);for k=1:Nruns [o2(:,:,k),Wh]=PLM_Hebb1(descr2(:,:,k),Mparams);if (mod(k,10)==5) fprintf('.');end;end;toc  % elapsed time 22 min
 cd('C:\work\PLExp2\modelfits');tic;st2=PLM_stats2(descr2,o2(:,1,:),o2(:,2,:),o2(:,3,:));toc     % elapsed_time = 6 min
 cd('C:\work\PLPaper2\simul');save st2 st2 ; 
 M=0.6;for k=1:5;subplot(3,2,k);plot(Wh([1:7]+7*k-7,:)');axis([0 9600 -M +M]);set(gca,'xtick',bnd,'xticklabel',[]);grid on;end;
 b=10;subplot(3,2,6);plot(1:9600/b,Wh(36,1:b:9600),'b.',1:9600/b,Wh(37,1:b:9600),'r.');axis([0 9600/b -.5 +.5]);set(gca,'xtick',bnd/b,'xticklabel',[]);grid on;

 clear all ; load st1;load st2;load('C:\work\PLModel1\HUMAN_ZCONGR.mat');
 fdbk_stats=merge_PLM_stats(st1,st2)
fdbk_stats = 
        group: [1000x1 double]
     zP_congr: [32x3x1000 double]
     zP_incon: [32x3x1000 double]
       dprime: [32x3x1000 double]
        accZP: [1000x6 double]
    mnA_congr: [32x3x1000 double]
    mnA_incon: [32x3x1000 double]
       Aprime: [32x3x1000 double]
         accA: [1000x6 double]
    mnI_congr: [32x3x1000 double]
    mnI_incon: [32x3x1000 double]
       Iprime: [32x3x1000 double]
         accI: [1000x6 double]

 fdbk_sst=PLM_summary(fdbk_stats)
fdbk_sst = 
       N_runs: 1000
     zP_congr: [32x3 double]
     zP_incon: [32x3 double]
       dprime: [32x3 double]
        accZP: [0.8776 0.9460 0.9846 1.2854 0.4847 -0.1380]
    mnA_congr: [32x3 double]
    mnA_incon: [32x3 double]
       Aprime: [32x3 double]
         accA: [0.4968 0.5122 0.5199 -0.0260 -0.2008 -0.3392]
    mnI_congr: [32x3 double]
    mnI_incon: [32x3 double]
       Iprime: [32x3 double]
         accI: [0.4968 0.5125 0.5197 -0.0261 -0.2010 -0.3390]

 figure(1);plot_PLM_fit({fdbk_sst.zP_congr',fdbk_sst.zP_incon'});
 figure(2);plot_PLM_dynamics(fdbk_sst);


***********************************************************
***
***  May 6, 2005 -- Consolidate all data
***

%% Retrieve z-Probability data from PLExp1, PLExp2, and the recent simulations
%
 clear all;load('C:\work\PLExp2\modelfits\HUMAN_ZCONGR.mat');load('C:\work\PLExp2\modelfits\MODEL_ZCONGR.mat');
 eZnf36=HUMAN_ZCONGR' ; mZnf36=MODEL_ZCONGR' ;  % [36x6] no feedback
 clear HUMAN_ZCONGR MODEL_ZCONGR model_sumstats
 load('C:\work\PLModel1\HUMAN_ZCONGR.mat');eZfk32=HUMAN_ZCONGR';clear HUMAN_ZCONGR;  % [32x6] empirical data with feedback
 load('C:\work\PLPaper2\simul\fdbk_stats.mat');
 mZfk32=zeros(32,6);mZfk32(:,[1 3 5])=fdbk_sst.zP_congr;mZfk32(:,[2 4 6])=fdbk_sst.zP_incon;
 clear fdbk_stats fdbk_sst

%% Empirical and model d-primes
%
 edprnf36=eZnf36(:,[1 3 5])+eZnf36(:,[2 4 6]);  % empirical, no feedback, 9+9 Ss
 mdprnf36=mZnf36(:,[1 3 5])+mZnf36(:,[2 4 6]);  % model, no feedback, 500+500 runs
 edprfk32=eZfk32(:,[1 3 5])+eZfk32(:,[2 4 6]);  % empirical, feedback, 13 Ss
 mdprfk32=mZfk32(:,[1 3 5])+mZfk32(:,[2 4 6]);  % model, feedback, 500+500 runs with the new learning rule (see the entry of May 2, 2005 in Hebb_fits.txt)

 save('C:\work\PLPaper2\simul\all_zP_data.mat');whos
  Name           Size         Bytes  Class
--------------------------------------------------
  eZfk32        32x6           1536  double array
  eZnf36        36x6           1728  double array
  edprfk32      32x3            768  double array
  edprnf36      36x3            864  double array
  mZfk32        32x6           1536  double array
  mZnf36        36x6           1728  double array
  mdprfk32      32x3            768  double array
  mdprnf36      36x3            864  double array
--------------------------------------------------


%% Descriptive statistics -- verify Table 3 in PLFeedback.doc
%
 PLM_sumsqfit2(mZnf36',eZnf36')
sumsq2 = 7.5112 = 0.8346 * 9 sessions, no feedback

 PLM_sumsqfit2(mZfk32',eZfk32')
sumsq2 = 6.0037 = 0.7505 * 8 sessions, with feedback

 describe(eZnf36);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % empirical, 9+9 observers
------------------------------------------------------------   % no feedback
   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

 describe(edprnf36./2);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % empirical, 9+9 observers
------------------------------------------------------------   % no feedback
   1.085    0.184      0.67    0.97    1.12    1.22    1.35    % Cp=.23
   0.792    0.139      0.39    0.72    0.83    0.89    0.98    % Cp=.15
   0.479    0.087      0.29    0.42    0.49    0.56    0.61    % Cp=.10
------------------------------------------------------------
   0.785    0.136      0.45    0.71    0.81    0.89    0.98

 describe(mZnf36);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % model, 500+500 runs
------------------------------------------------------------   % no feedback
   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

 describe(mdprnf36./2);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % model, 500+500 runs
------------------------------------------------------------   % no feedback
   1.146    0.137      0.86    1.05    1.15    1.26    1.34    % Cp=.23
   0.753    0.095      0.56    0.68    0.76    0.83    0.89    % Cp=.15
   0.443    0.061      0.32    0.39    0.45    0.49    0.55    % Cp=.10
------------------------------------------------------------
   0.781    0.098      0.58    0.71    0.79    0.86    0.93


 describe(eZfk32);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % empirical, 13 observers
------------------------------------------------------------   % feedback
   0.821    0.150      0.47    0.70    0.85    0.93    1.05    % congr, Cp=.23
   1.126    0.242      0.48    0.98    1.14    1.30    1.50    % incon, Cp=.23
   0.956    0.179      0.60    0.87    0.98    1.07    1.46    % congr, Cp=.15
   0.565    0.151      0.16    0.48    0.57    0.67    0.81    % incon, Cp=.15
   1.031    0.207      0.63    0.90    1.04    1.20    1.34    % congr, Cp=.10
  -0.070    0.078     -0.21   -0.13   -0.08   -0.00    0.12    % incon, Cp=.10
------------------------------------------------------------
   0.738    0.168      0.35    0.63    0.75    0.86    1.05

 describe(edprfk32./2);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % empirical, 13 observers
------------------------------------------------------------   % feedback
   0.974    0.181      0.47    0.88    1.01    1.12    1.25    % Cp=.23
   0.761    0.150      0.40    0.69    0.80    0.86    1.11    % Cp=.15
   0.481    0.112      0.26    0.42    0.48    0.59    0.66    % Cp=.10
------------------------------------------------------------
   0.738    0.148      0.38    0.67    0.76    0.86    1.01

 describe(mZfk32);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % model, 500+500 runs under runav-corrected learning rule
------------------------------------------------------------   % feedback
   0.878    0.147      0.60    0.72    0.90    0.98    1.13    % congr, Cp=.23
   1.285    0.233      0.80    1.11    1.33    1.45    1.69    % incon, Cp=.23
   0.946    0.133      0.68    0.81    0.97    1.04    1.15    % congr, Cp=.15
   0.485    0.144      0.19    0.38    0.50    0.60    0.74    % incon, Cp=.15
   0.985    0.122      0.73    0.88    1.00    1.06    1.16    % congr, Cp=.10
  -0.138    0.078     -0.33   -0.19   -0.12   -0.08   -0.00    % incon, Cp=.10
------------------------------------------------------------
   0.740    0.143      0.44    0.62    0.76    0.84    0.98

 describe(mdprfk32./2);
    Mean  Std.dev       Min     Q25  Median     Q75     Max    % model, 500+500 runs under runav-corrected learning rule
------------------------------------------------------------   % feedback
   1.082    0.164      0.72    0.99    1.12    1.20    1.30    % Cp=.23
   0.715    0.115      0.45    0.64    0.74    0.80    0.87    % Cp=.15
   0.423    0.074      0.25    0.38    0.44    0.48    0.53    % Cp=.10
------------------------------------------------------------
   0.740    0.118      0.48    0.67    0.77    0.83    0.90


%%% Correlation coefficients and RMSEs
 e=edprfk32;m=mdprfk32;c=corrcoef(e(:),m(:));R2=c(2)^2 , rmse=sqrt(mean((m(:)-e(:)).^2))
R2 = 0.8822 ,  rmse = 0.2097     % with feedback, 96 data points total

 e=edprnf36;m=mdprnf36;c=corrcoef(e(:),m(:));R2=c(2)^2 , rmse=sqrt(mean((m(:)-e(:)).^2))
R2 = 0.9261 ,  rmse = 0.1670     % without feedback, 108 data points total



***********************************************************
***
***  May 6, 2005 -- Make Figures 8 and 9 (model fits) in PLFeedback.doc
***

%%  No-feedback d' fits -- top panel of Figure 8 in PLFeedback.doc
 subplot(2,1,1);e=edprnf36;m=mdprnf36;t={[1],[2:9]',[10:17]',[18:25]',[26:33]',[34:36]'};T=[1:36]';
 h=plot(T,e(:,1),'b^',T,e(:,2),'ro',T,e(:,3),'bs',T,m(:,1),'k^',T,m(:,2),'ko',T,m(:,3),'ks',t{2},m(t{2},:),'k-',t{3},m(t{3},:),'k-',t{4},m(t{4},:),'k-',t{5},m(t{5},:),'k-',t{6},m(t{6},:),'k-');
 set(h,'LineWidth',1);set(h(4:6),'MarkerSize',4);   % keep open symbols
 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  (1000 runs)');

%%  New feedback d' fits with the runav-centered learning rule -- bottom panel of Figure 8 in PLFeedback.doc
%%  Re-print of Figure 4 in the Psych Review paper.
 subplot(2,1,2);e=edprfk32;m=mdprfk32;t={[1],[2:9]',[10:17]',[18:25]',[26:31]',[32]};T=[1:32]';
 h=plot(T,e(:,1),'b^',T,e(:,2),'ro',T,e(:,3),'bs',T,m(:,1),'k^',T,m(:,2),'ko',T,m(:,3),'ks',t{2},m(t{2},:),'k-',t{3},m(t{3},:),'k-',t{4},m(t{4},:),'k-',t{5},m(t{5},:),'k-');
 set(h,'LineWidth',1);set(h(4:6),'MarkerSize',4);   % keep open symbols
 axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);xlabel('Block number  (300 trials/block, 4 blocks/day)');ylabel('d''');title('Trial-by-trial feedback  (1000 runs)');

% Set page size=[1 1 6.5 8 ]; % [Left Bottom Width Height]
% See PLExp2\PLE2_results.txt for the corresponding regression figures

%%  No-feedback incongruent fits -- top panel of Figure 9 in PLFeedback.doc
 subplot(2,1,1);e=eZnf36(:,[2 4 6]);m=mZnf36(:,[2 4 6]);t={[1],[2:9]',[10:17]',[18:25]',[26:33]',[34:36]'};T=[1:36]';
 h=plot(T,e(:,1),'b^',T,e(:,2),'ro',T,e(:,3),'bs',t{2},m(t{2},:),'k-',T,m(:,1),'k^',T,m(:,2),'ko',T,m(:,3),'ks',t{3},m(t{3},:),'k-',t{4},m(t{4},:),'k-',t{5},m(t{5},:),'k-',t{6},m(t{6},:),'k-');
 set(h,'LineWidth',1);set(h(7:9),'MarkerSize',4);   % keep open symbols
 axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);xlabel('Block number  (300 trials/block, 4 blocks/day)');ylabel('z-probability correct');
 title('Incongruent 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]);

%%  No-feedback congruent fits -- bottom panel of Figure 9 in PLFeedback.doc
 subplot(2,1,2);e=eZnf36(:,[1 3 5]);m=mZnf36(:,[1 3 5]);t={[1],[2:9]',[10:17]',[18:25]',[26:33]',[34:36]'};T=[1:36]';
 h=plot(T,e(:,1),'b^',T,e(:,2),'ro',T,e(:,3),'bs',t{2},m(t{2},:),'k-',T,m(:,1),'k^',T,m(:,2),'ko',T,m(:,3),'ks',t{3},m(t{3},:),'k-',t{4},m(t{4},:),'k-',t{5},m(t{5},:),'k-',t{6},m(t{6},:),'k-');
 set(h,'LineWidth',1);set(h(7:9),'MarkerSize',4);   % keep open symbols
 axis([0 37 0 3]);set(gca,'xtick',[0:4:36]);xlabel('Block number  (300 trials/block, 4 blocks/day)');ylabel('z-probability correct');
 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]);legend('Data, c=0.245','Data, c=0.160','Data, c=0.106','Model fits',3);

%% Now the fun part --- the inset
 e=[mean(eZnf36(:,[2 4 6])) ; mean(eZnf36(:,[1 3 5]))]
e = 0.9630    0.2839   -0.4317
    1.2077    1.2998    1.3902
 m=[mean(mZnf36(:,[2 4 6])) ; mean(mZnf36(:,[1 3 5]))]
m = 1.0694    0.1874   -0.4775
    1.2233    1.3194    1.3642

 ah=axes('position',[0.67 0.12 0.22 0.14]);t=[.75 2.25]';T=[1 2];
 h=plot(t,e(:,1),'b^',t,e(:,2),'ro',t,e(:,3),'bs',T,m(:,1),'k^-',T,m(:,2),'ko-',T,m(:,3),'ks-');
 set(h,'LineWidth',1);set(h(4:6),'MarkerSize',4);   % keep open symbols
 axis([.25 2.75 -.75 1.75]);set(gca,'xtick',[1 2],'xticklabel','Incon|Congr','XAxisLocation','top','ytick',[-.5:.5:1.5],'yticklabel','|0||1|');

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


************  2005-07-05  New fits under PLM_Hebb2 ver 2.1
See PLPaper2/simul/FinalSumSq.txt and C:\work\PLModel1\newsimul\Hebb2m_fits.txt for details
