********  HRCR_MSTATS.TXT -- 28 Dec 2003 - 4 Jan 2004  ********
**

P  -->  PS  -->  PSC  -->  PSCR  -->   PSACR=ANCHOR
            \->  PSA  -->  PSAC  /

Model P:     Perception mechanism only; nearest-anchor responses.
Model PS:    Perception + anchor Selection mechanisms only.
Model PSC:   Perception + Selection + Correction mechanisms only.
Model PSA:   Perception + Selection + Activation, no corrections.
Model PSAC:  All except location learning.
Model PSCR:  All except activation learning.
Model PSACR: Full ANCHOR model.

The general methodology is always the same:
1. Fit the same model to the 40 individual data sets from the
   Category Rating Expmt 2. This produces 40 parameter sets.
   Which parameters are "active", which are shunted off, and
   which are fixed for all subjects depends on the particular model.
3. Run the model with each CR parametrization on each of the
   new CR stimulus sequences, for a total of 9600=40*240 runs.
4. Calculate a battery of model statistics via CR_STATS and
   compare them with the empirical ones.

This file carries out steps 2. and 3. above.
See HRCR_PSEARCH.TXT for details on the parameter search (step 1.)
The results from each parameter search are stored in HIER_OPT_PARAMS.MAT.
The summary statistics for each model are stored in HIER_SUMSTATS.MAT.
See also HRAI_PSEARCH.TXT AND HRAI_MSTATS.TXT for identification analogs.

WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING


*************************************************************
***                ******************************************
***   Stimuli      ******************************************
***                ******************************************
*************************************************************

%%%%%  Generate 9600=24*400 absolute identification sequences.
%%     Each 'pseudosubject' gets the same 400 sequences,
%%     repeated for each model on the hierarchy.
%%     200 sequences according to Schedule 1 and 200 to Schedule 2.
%%     These are codes '201' and '202' in GENERATE_STIMULI.M.

 clear all ; cd 'C:\work\anchor\finalsim\' ;
 AIstim1=reshape(generate_stimuli(struct('N_runs',24*200,'schedule',201)),[450,200,24]);
 AIstim2=reshape(generate_stimuli(struct('N_runs',24*200,'schedule',202)),[450,200,24]);

%%%%%  Generate 9600=40*240 category rating sequences.
%%     These are codes '101' and '102' in GENERATE_STIMULI.M.
 CRstim1=reshape(generate_stimuli(struct('N_runs',40*120,'schedule',101)),[450,120,40]);
 CRstim2=reshape(generate_stimuli(struct('N_runs',40*120,'schedule',102)),[450,120,40]);

 whos
  Name          Size         Bytes  Class
-------------------------------------------------
  AIstim1     450x200x24  17280000  double array
  AIstim2     450x200x24  17280000  double array
  CRstim1     450x120x40  17280000  double array
  CRstim2     450x120x40  17280000  double array
-------------------------------------------------
Grand total is 8640000 elements using 69120000 bytes


*************************************************************
***                ******************************************
***   Model P      ******************************************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
28 Dec 2003

Active parameters:  perc_k, L1
Shunted parameters: mem_k=0, temper=0, hist=0, cutoff=Inf, alpha=0
Fixed parameters:   perc_k=0.04, decay=0.5, L9=0.675, avail=[1:9]'

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_P, among others

%% Optimal CR params
 resCR=resCR40_P ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.perc_k = VAL / 10 ;'
      'PARAMS.anchors(:,1) = [VAL:(0.675-VAL)/8:0.675]' ;'

 p0=resCR.search_params.params;x=resAI.opt_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.210    0.252      0.83    0.97    1.18    1.43    1.67  % perc_k * 10
   0.150    0.063      0.05    0.10    0.15    0.20    0.25  % L1

% Due to an unfortunate programming shortcut, ANCHOR2.M multiplies
% all anchor magnitudes by 0.0 when mem_k=0. Setting it to 1e-10
% achieves the desired effect.
% Letting T=0 does not cause division by zero -- see SOFTMAX1.M.

 p0.mem_k=1e-10 ;
 for k=1:40 CR_params(k)=vector2params(x(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0867
       mem_k: 1.0000e-010
       avail: [9x1 double]
     anchors: [9x6 double]
     cutoffs: [-Inf -Inf Inf Inf]
      temper: 0
     history: 0
       alpha: 0
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.1066 0.1777 0.2487 0.3198 0.3908 0.4619 0.5329 0.6040 0.6750]

%% Run the model
 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 850.5930 sec (15 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_P=CR_stats(s,mr);toc

% Elapsed_time = 73 sec on a Dell Optiplex GX400.

 sstCR40_P=CR_summary(stCR_P)
sstCR40_P =
    resp_mean: [5.5580 0.5180 5.8355 0.3756 4.8311 5.8422 6.6044]
     resp_std: [1.7720 0.2440 2.0810 0.2365 1.6100 2.0767 2.6963]
           R2: [0.7750 0.0820 0.8240 0.0560 0.6439 0.8290 0.9252]
       deltaS: [0.5533 0.3520 8.1430e-4 0.1462 -0.5347 1.4943e-4 0.5580]
      deltaR2: [-0.0160 0.1000 5.1584e-4 0.0495 -0.2627 8.2064e-5 0.2552]
     deltaARL: [0.4870 0.5830 0.0013 0.1322 -0.4929 5.3270e-4 0.5769]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 -0.0035 0.2424 -1.4123 5.8337e-4 1.0775]
           d1: [0.4510 0.3440 -4.1078e-4 0.2408 -1.4123 0.0053 0.9447]
           d2: [-0.0350 0.5120 -0.0065 0.2439 -0.9423 -0.0036 1.0775]
     resp_acf: [0.1950 0.1000 0.1073 0.0499 -0.1185 0.1078 0.2860]
    resid_acf: [0.3430 0.1160 -0.0027 0.0471 -0.2106 -0.0020 0.1909]
    resid_std: [0.8090 0.1240 0.8590 0.1590 0.5531 0.8403 1.3573]
    mn_blk_R2: [0.8238 0.7585 0.8238 0.7580 0.8243]
    mn_blkstd: [2.0691 1.7607 2.0683 1.7604 2.0683]
       mn_ARL: [5.836 5.837 5.880 5.881 5.836 5.834 5.880 5.881 5.837 5.840;
                5.835 5.837 5.877 5.877 5.835 5.837 5.882 5.883 5.837 5.836]

 a=sstCR40_P.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_P sstCR40_P


*************************************************************
***                ******************************************
***   Model PS     ******************************************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
29 Dec 2003

Active parameters:  mem_k, temper, L1
Shunted parameters: hist=0, cutoff=Inf, alpha=0
Fixed parameters:   perc_k=0.04, decay=0.5, L9=0.675, avail=[1:9]'

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_PS, among others

%% Optimal CR params
 resCR=resCR40_PS ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.mem_k = VAL / 10 ;'
      'PARAMS.temper = VAL / 10 ;'
      'PARAMS.anchors(:,1) = [VAL:(0.675-VAL)/8:0.675]' ;'

 p0=resCR.search_params.params;x=resCR.opt_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   1.433    0.502      0.20    1.07    1.48    1.96    2.00  % mem_k * 10
   0.271    0.116      0.10    0.18    0.27    0.36    0.50  % temper* 10
   0.155    0.065      0.05    0.10    0.17    0.21    0.26  % L1

 for k=1:40 CR_params(k)=vector2params(x(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0400
       mem_k: 0.0967
       avail: [9x1 double]
     anchors: [9x6 double]
     cutoffs: [-Inf -Inf Inf Inf]
      temper: 0.0192
     history: 0
       alpha: 0
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.1097 0.1803 0.2510 0.3217 0.3923 0.4630 0.5337 0.6043 0.6750]

%% Run the model
 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 4.518e+003 sec (75 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_PS=CR_stats(s,mr);toc

% Elapsed_time = 503.1330 sec (9 min) on a Dell Optiplex GX400.

 sstCR40_PS=CR_summary(stCR_PS)
sstCR40_PS =
    resp_mean: [5.5580 0.5180 5.8350 0.3861 4.8289 5.8400 6.6689]
     resp_std: [1.7720 0.2440 2.0606 0.2212 1.5976 2.0608 2.5735]
           R2: [0.7750 0.0820 0.6848 0.1085 0.3677 0.6871 0.9104]
       deltaS: [0.5533 0.3520 2.3482e-4 0.1541 -0.6885 0.0013 0.5976]
      deltaR2: [-0.0160 0.1000 7.2628e-6 0.0821 -0.4305 6.6577e-4 0.3736]
     deltaARL: [0.4870 0.5830 -8.8815e-4 0.1748 -0.7238 1.3612e-4 0.8465]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 -0.0130 0.3216 -1.4289 -0.0063 1.2621]
           d1: [0.4510 0.3440 -0.0150 0.3197 -1.4289 -0.0097 1.1211]
           d2: [-0.0350 0.5120 -0.0110 0.3235 -1.4096 -0.0035 1.2621]
     resp_acf: [0.1950 0.1000 0.0887 0.0515 -0.1231 0.0896 0.2653]
    resid_acf: [0.3430 0.1160 -0.0035 0.0471 -0.1766 -0.0042 0.1675]
    resid_std: [0.8090 0.1240 1.1390 0.2482 0.6871 1.1133 1.8489]
    mn_blk_R2: [0.6854 0.5974 0.6848 0.5974 0.6854]
    mn_blkstd: [2.0490 1.7960 2.0468 1.7955 2.0487]
       mn_ARL: [5.825 5.829 5.909 5.917 5.827 5.823 5.905 5.902 5.824 5.824;
                5.823 5.822 5.908 5.903 5.819 5.827 5.911 5.914 5.825 5.821]

 a=sstCR40_PS.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_PS sstCR40_PS


*************************************************************
***                ******************************************
***   Model PSC    ******************************************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
31 Dec 2003

Active parameters:  mem_k, temper, cutoff, L1
Shunted parameters: hist=0, alpha=0
Fixed parameters:   perc_k=0.04, decay=0.5, L9=0.675, avail=[1:9]',
                    cutoff multipliers=[-3 -1 +0.9 +2.7]

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_PSC, among others

%% Optimal CR params
 resCR=resCR40_PSC ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.mem_k = VAL / 10 ;'
      'PARAMS.temper = VAL / 10 ;'
      'PARAMS.cutoffs = VAL*[-3 -1 +0.9 +2.7] ;'
      'PARAMS.anchors(:,1) = [VAL:(0.675-VAL)/8:0.675]' ;'

 p0=resCR.search_params.params;x=resCR.opt_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   0.920    0.207      0.22    1.00    1.00    1.00    1.00  % mem_k * 10
   0.619    0.250      0.20    0.39    0.80    0.80    0.80  % temper* 10
   0.767    0.398      0.50    0.50    0.50    1.07    1.50  % cutoff
   0.164    0.068      0.05    0.11    0.17    0.22    0.26  % L1

 for k=1:40 CR_params(k)=vector2params(x(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0400
       mem_k: 0.1000
       avail: [9x1 double]
     anchors: [9x6 double]
     cutoffs: [-1.9774 -0.6591 0.5932 1.7796]
      temper: 0.0200
     history: 0
       alpha: 0
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.1169 0.1867 0.2564 0.3262 0.3959 0.4657 0.5355 0.6052 0.6750]

%% Run the model
 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 874.8780 sec (15 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_PSC=CR_stats(s,mr);toc

% Elapsed_time = 72.8750 sec (2 min) on a Dell Optiplex GX400.

 sstCR40_PSC=CR_summary(stCR_PSC)
sstCR40_PSC =
    resp_mean: [5.5580 0.5180 5.8122 0.4213 4.7756 5.8089 6.6778]
     resp_std: [1.7720 0.2440 2.0484 0.2160 1.5912 2.0552 2.5808]
           R2: [0.7750 0.0820 0.7992 0.0480 0.6168 0.7997 0.9261]
       deltaS: [0.5533 0.3520 0.0018 0.1484 -0.5895 0.0045 0.5668]
      deltaR2: [-0.0160 0.1000 -0.0005 0.0587 -0.2657 -0.0002 0.2585]
     deltaARL: [0.4870 0.5830 0.0018 0.1373 -0.6042 0.0022 0.7447]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 -0.0115 0.2522 -1.4866 -0.0102 1.0544]
           d1: [0.4510 0.3440 -0.0084 0.2523 -1.4866 -0.0069 1.0544]
           d2: [-0.0350 0.5120 -0.0145 0.2522 -1.2914 -0.0149 1.0130]
     resp_acf: [0.1950 0.1000 0.1039 0.0498 -0.0948 0.1047 0.2708]
    resid_acf: [0.3430 0.1160 -0.0028 0.0472 -0.1881 -0.0032 0.1879]
    resid_std: [0.8090 0.1240 0.9120 0.1559 0.6348 0.8960 1.3804]
    mn_blk_R2: [0.7989 0.7272 0.7994 0.7266 0.7985]
    mn_blkstd: [2.0371 1.7420 2.0360 1.7426 2.0353]
       mn_ARL: [5.808 5.809 5.861 5.863 5.807 5.806 5.853 5.855 5.810 5.808;
                5.807 5.805 5.858 5.848 5.809 5.808 5.861 5.862 5.808 5.810]

 a=sstCR40_PSC.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_PSC sstCR40_PSC


*************************************************************
***                ******************************************
***   Model PSA    ******************************************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
31 Dec 2003

Active parameters:  mem_k, temper, hist, L1
Shunted parameters: cutoff=Inf, alpha=0
Fixed parameters:   perc_k=0.04, decay=0.5, L9=0.675, avail=[1:9]',
                    cutoff multipliers=[-3 -1 +0.9 +2.7]

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_PSA, among others

%% Optimal CR params
 resCR=resCR40_PSA ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.mem_k = VAL / 10 ;'
      'PARAMS.temper = VAL / 10 ;'
      'PARAMS.history = VAL / 5 ;'
      'PARAMS.anchors(:,1) = [VAL:(0.675-VAL)/8:0.675]' ;'

 p0=resCR.search_params.params;x=resCR.opt_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   0.910    0.189      0.20    0.92    1.00    1.00    1.00  % mem_k * 10
   0.320    0.089      0.20    0.24    0.31    0.38    0.49  % temper* 10
   0.227    0.097      0.00    0.17    0.24    0.30    0.40  % history* 5
   0.179    0.067      0.05    0.13    0.19    0.23    0.29  % L1

 for k=1:40 CR_params(k)=vector2params(x0(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0400
       mem_k: 0.0926
       avail: [9x1 double]
     anchors: [9x6 double]
     cutoffs: [-Inf -Inf Inf Inf]
      temper: 0.0200
     history: 0.0200
       alpha: 0
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.1158 0.1857 0.2556 0.3255 0.3954 0.4653 0.5352 0.6051 0.6750]

%% Run the model
%  Note that without corrections, model PSA is prone to unstable
%  runaway activations in the absence of feedback. These runs
%  are likely to be very lopsided.

 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 936.7570 sec (16 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_PSA=CR_stats(s,mr);toc

% Elapsed_time = 74.1170 sec (2 min) on a Dell Optiplex GX400.
% Numerous "Divide by zero" warnings:
%> In C:\MATLABR11\toolbox\matlab\datafun\corrcoef.m at line 31
%  In C:\work\anchor\finalsim\CR_stats.m at line 71  % or 69, 70, 72
% This happens when CR_STATS attempts to calculate the S-R
% correlations and the response variance is zero.

 sstCR40_PSA=CR_summary(stCR_PSA)
sstCR40_PSA =
    resp_mean: [5.5580 0.5180 5.7045 0.4651 3.7956 5.7089 6.8356]
     resp_std: [1.7720 0.2440 1.9640 0.3510 0 1.9823 3.9947]
           R2: [0.7750 0.0820 NaN NaN 0.0095 NaN 0.9078]
       deltaS: [0.5533 0.3520 0.0515 0.1847 -0.9065 0.0482 1.0419]
      deltaR2: [-0.0160 0.1000 NaN NaN -0.4050 NaN 0.4316]
     deltaARL: [0.4870 0.5830 -0.0022 0.1940 -0.9414 -4.6375e-004 0.8501]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 -0.1575 0.3243 -2.0746 -0.1360 1.1340]
           d1: [0.4510 0.3440 -0.1573 0.3245 -2.0106 -0.1321 1.1340]
           d2: [-0.0350 0.5120 -0.1576 0.3241 -2.0746 -0.1387 0.9831]
     resp_acf: [0.1950 0.1000 NaN NaN -0.0815 NaN 0.3040]
    resid_acf: [0.3430 0.1160 NaN NaN -0.1698 NaN 0.2797]
    resid_std: [0.8090 0.1240 NaN NaN 0.2287 NaN 2.2796]
    mn_blk_R2: [NaN NaN NaN NaN NaN]
    mn_blkstd: [1.9871 1.6551 1.9329 1.6903 1.9356]
       mn_ARL: [5.696 5.698 5.825 5.897 5.817 5.800 5.787 5.739 5.715 5.712;
                5.698 5.702 5.640 5.563 5.579 5.597 5.674 5.721 5.682 5.677]

 a=sstCR40_PSA.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_PSA sstCR40_PSA


*************************************************************
***                ******************************************
***   Model PSAC   ******************************************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
4 Jan 2004

Active parameters:  mem_k, temper, hist, cutoff, L1
Shunted parameters: alpha=0
Fixed parameters:   perc_k=0.04, decay=0.5, L9=0.675, avail=[1:9]'
                    cutoff multipliers=[-3 -1 +0.9 +2.7]

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_PSAC, among others

%% Optimal CR params
 resCR=resCR40_PSAC ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.mem_k = VAL / 10 ;'
      'PARAMS.temper = VAL / 10 ;'
      'PARAMS.history = VAL / 5 ;'
      'PARAMS.cutoffs = VAL*[-3 -1 +0.9 +2.7] ;'
      'PARAMS.anchors(:,1) = [VAL:(0.675-VAL)/8:0.675]' ;'

 p0=resCR.search_params.params;x=resCR.corr_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   0.855    0.150      0.22    0.88    0.92    0.92    0.92  % mem_k * 10
   0.542    0.082      0.40    0.47    0.60    0.60    0.60  % temper* 10
   0.751    0.221      0.25    0.64    0.86    0.94    0.94  % history* 5
   1.020    0.109      0.50    1.04    1.07    1.07    1.07  % cutoff
   0.200    0.068      0.05    0.15    0.21    0.25    0.32  % L1

 for k=1:40 CR_params(k)=vector2params(x(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0400
       mem_k: 0.0839
       avail: [9x1 double]
     anchors: [9x6 double]
     cutoffs: [-2.7684 -0.9228 0.8305 2.4916]
      temper: 0.0400
     history: 0.1107
       alpha: 0
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.1323 0.2002 0.2680 0.3358 0.4037 0.4715 0.5393 0.6072 0.6750]

%% Run the model
 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 873.9660 sec (15 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_PSAC=CR_stats(s,mr);toc

% Elapsed_time = 81.4870 sec (2 min) on a Dell Optiplex GX400.

 sstCR40_PSAC=CR_summary(stCR_PSAC)
sstCR40_PSAC = 
    resp_mean: [5.5580 0.5180 5.6102 0.5207 3.2289 5.6178 6.8422]
     resp_std: [1.7720 0.2440 1.9104 0.2750 1.3198 1.8871 3.0789]
           R2: [0.7750 0.0820 0.8184 0.0408 0.6260 0.8182 0.9169]
       deltaS: [0.5533 0.3520 -0.0129 0.2092 -0.8277 -0.0150 1.1483]
      deltaR2: [-0.0160 0.1000 0.0154 0.0520 -0.2528 0.0135 0.2412]
     deltaARL: [0.4870 0.5830 -0.0061 0.2141 -1.6827 -0.0046 1.5550]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 -0.2546 0.2723 -1.7293 -0.2369 0.6989]
           d1: [0.4510 0.3440 -0.2606 0.2731 -1.7293 -0.2419 0.6989]
           d2: [-0.0350 0.5120 -0.2486 0.2714 -1.4776 -0.2319 0.5813]
     resp_acf: [0.1950 0.1000 0.1418 0.0525 -0.0397 0.1424 0.3143]
    resid_acf: [0.3430 0.1160 0.0400 0.0559 -0.1874 0.0391 0.3847]
    resid_std: [0.8090 0.1240 0.8051 0.1292 0.5529 0.7927 1.4949]
    mn_blk_R2: [0.8101 0.7586 0.8147 0.7522 0.8255]
    mn_blkstd: [1.8932 1.5606 1.8628 1.6118 1.9061]
       mn_ARL: [5.605 5.605 5.777 5.902 5.834 5.799 5.721 5.641 5.637 5.638;
                5.605 5.604 5.480 5.351 5.376 5.407 5.519 5.600 5.562 5.558]

 a=sstCR40_PSAC.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_PSAC sstCR40_PSAC


*************************************************************
***                ******************************************
***   Model PSCR   ******************************************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
2 Jan 2004

Active parameters:  mem_k, temper, cutoff
Shunted parameters: hist=0
Fixed parameters:   perc_k=0.04, alpha=0.3, decay=0.5,
                    L1=0.275, L9=0.675, avail=[1 ; 9] only
                    cutoff multipliers=[-3 -1 +0.9 +2.7]

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_PSCR, among others

%% Optimal CR params
 resCR=resCR40_PSCR ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.mem_k = VAL / 10 ;'
      'PARAMS.temper = VAL / 10 ;'
      'PARAMS.cutoffs = VAL*[-3 -1 +0.9 +2.7] ;'

% Convert the CORRECTED parameters. See the post-optimization
% correction for Models PSCR and PSACR in HRCR_PSEARCH.TXT.
 p0=resCR.search_params.params;x=resCR.corr_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   0.750    0.172      0.34    0.61    0.80    0.92    0.92  % mem_k * 10
   0.517    0.086      0.40    0.40    0.53    0.60    0.60  % temper* 10
   0.873    0.189      0.50    0.76    0.89    1.07    1.07  % cutoff

 for k=1:40 CR_params(k)=vector2params(x(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0400
       mem_k: 0.0800
       avail: [2x1 double]
     anchors: [9x6 double]
     cutoffs: [-2.4052 -0.8017 0.7216 2.1647]
      temper: 0.0552
     history: 0
       alpha: 0.3000
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.2750 0.3250 0.3750 0.4250 0.4750 0.5250 0.5750 0.6250 0.6750]

%% Run the model
 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 898.0510 sec (15 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_PSCR=CR_stats(s,mr);toc

% Elapsed_time = 64.3120 sec (1 min) on a Dell Optiplex GX400.

 sstCR40_PSCR=CR_summary(stCR_PSCR)
sstCR40_PSCR = 
    resp_mean: [5.5580 0.5180 5.0782 0.3151 3.6156 5.0711 6.4667]
     resp_std: [1.7720 0.2440 2.5653 0.0933 2.2325 2.5637 2.9350]
           R2: [0.7750 0.0820 0.7341 0.1178 0.0844 0.7604 0.9268]
       deltaS: [0.5533 0.3520 0.1203 0.2438 -0.9873 0.1227 1.1155]
      deltaR2: [-0.0160 0.1000 0.1128 0.1434 -0.3902 0.0839 0.7404]
     deltaARL: [0.4870 0.5830 0.0816 0.6914 -2.8977 0.1016 2.6891]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 0.6211 0.5548 -1.7616 0.6107 2.9402]
           d1: [0.4510 0.3440 0.6631 0.5490 -1.5826 0.6637 2.9392]
           d2: [-0.0350 0.5120 0.5791 0.5575 -1.7616 0.5682 2.9402]
     resp_acf: [0.1950 0.1000 0.0945 0.0576 -0.1015 0.0938 0.3600]
    resid_acf: [0.3430 0.1160 0.0828 0.0695 -0.1567 0.0789 0.4116]
    resid_std: [0.8090 0.1240 1.2964 0.2939 0.6686 1.2497 2.5637]
    mn_blk_R2: [0.6961 0.6457 0.7698 0.7281 0.8088]
    mn_blkstd: [2.6252 2.2713 2.5384 2.1711 2.5049]
       mn_ARL: [5.033 5.025 4.933 4.701 4.673 4.793 5.052 5.364 5.358 5.282;
                5.028 5.028 5.237 5.519 5.465 5.352 5.211 4.940 4.864 4.936]

 a=sstCR40_PSCR.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_PSCR sstCR40_PSCR


*************************************************************
***                ******************************************
***   Model PSACR  = Full ANCHOR Model **********************
***   Cat. Rat     ******************************************
***                ******************************************
*************************************************************
2 Jan 2004

Active parameters:  mem_k, temper, hist, cutoff, L1
Shunted parameters: alpha=0
Fixed parameters:   perc_k=0.04, decay=0.5, L9=0.675, avail=[1:9]'
                    cutoff multipliers=[-3 -1 +0.9 +2.7]

 cd 'C:\work\anchor\finalsim\' ; load hier_stimuli ;
 load hier_opt_params ;      % resCR40_PSACR, among others

%% Optimal CR params
 resCR=resCR40_PSACR ; v2p=resCR.search_params.v2p_templ;v2p'
ans = 'PARAMS.mem_k = VAL / 10 ;'
      'PARAMS.temper = VAL / 10 ;'
      'PARAMS.history = VAL / 5 ;'
      'PARAMS.cutoffs = VAL*[-3 -1 +0.9 +2.7] ;'

% Convert the CORRECTED parameters. See the post-optimization
% correction for Models PSCR and PSACR in HRCR_PSEARCH.TXT.
 p0=resCR.search_params.params;x=resCR.corr_X;describe(x)
    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
   0.762    0.159      0.35    0.65    0.76    0.92    0.92  % mem_k * 10
   0.495    0.081      0.40    0.42    0.49    0.60    0.60  % temper* 10
   0.534    0.187      0.19    0.41    0.49    0.67    0.94  % history* 5
   0.864    0.153      0.50    0.75    0.83    1.03    1.07  % cutoff

 for k=1:40 CR_params(k)=vector2params(x(k,:),p0,v2p);end;CR_params(40)
CR_params(40) =
       scale: 'LINEAR'
       N_cat: 9
     SM_conv: 0.0010
      cat_sz: 0.0500
      perc_k: 0.0400
       mem_k: 0.0828
       avail: [2x1 double]
     anchors: [9x6 double]
     cutoffs: [-2.4298 -0.8099 0.7290 2.1869]
      temper: 0.0423
     history: 0.0945
       alpha: 0.3000
       decay: 0.5000
         ITI: 4
    M_raster: 7
    A_raster: [5 5 3 3]
    mnfieldp: 1

 CR_params(40).anchors(:,1)'
ans = [0.2750 0.3250 0.3750 0.4250 0.4750 0.5250 0.5750 0.6250 0.6750]

%% Run the model
 CRmresp1=zeros(size(CRstim1));CRmresp2=zeros(size(CRstim2));
 tic; for k=1:40;fprintf('%d ',k);p=CR_params(k);
    CRmresp1(:,:,k)=anchor2(CRstim1(:,:,k),p);
    CRmresp2(:,:,k)=anchor2(CRstim2(:,:,k),p); end;toc

% Elapsed_time = 880.0760 sec (15 min) on a Dell Optiplex GX400.

%% Calculate statistics
 s=[reshape(CRstim1,450,120*40),reshape(CRstim2,450,120*40)];
 mr=[reshape(CRmresp1,450,120*40),reshape(CRmresp2,450,120*40)];
 tic;stCR_PSACR=CR_stats(s,mr);toc

% Elapsed_time = 67.4670 sec (1 min) on a Dell Optiplex GX400.

 sstCR40_PSACR=CR_summary(stCR_PSACR)
sstCR40_PSACR = 
    resp_mean: [5.5580 0.5180 5.3065 1.4930 2.2000 5.3056 8.0356]
     resp_std: [1.7720 0.2440 1.9301 0.5466 0.9150 1.7625 3.4519]
           R2: [0.7750 0.0820 0.7666 0.0885 0.3118 0.7796 0.9324]
       deltaS: [0.5533 0.3520 0.2070 0.3737 -0.8092 0.1375 2.2110]
      deltaR2: [-0.0160 0.1000 0.0318 0.1002 -0.5039 0.0240 0.7189]
     deltaARL: [0.4870 0.5830 0.2703 0.6421 -3.6024 0.2345 3.0484]
WARNING: The sign of the context measure d is reversed in the manuscript.  WARNING
            d: [0.2080 0.4960 0.5256 0.5113 -2.4245 0.5194 2.7662]
           d1: [0.4510 0.3440 0.6447 0.4981 -2.4245 0.6284 2.7662]
           d2: [-0.0350 0.5120 0.4064 0.4963 -1.4506 0.4053 2.6814]
     resp_acf: [0.1950 0.1000 0.1047 0.0593 -0.0974 0.1034 0.4647]
    resid_acf: [0.3430 0.1160 0.1684 0.1131 -0.1406 0.1516 0.7552]
    resid_std: [0.8090 0.1240 0.9129 0.3133 0.4073 0.8177 2.3398]
    mn_blk_R2: [0.7831 0.7503 0.8099 0.7501 0.8149]
    mn_blkstd: [2.0513 1.6583 1.8899 1.6224 1.8443]
       mn_ARL: [5.092 5.166 5.202 5.163 5.191 5.326 5.567 5.808 5.770 5.696;
                5.117 5.193 5.337 5.479 5.418 5.333 5.246 5.073 5.045 5.138]

 a=sstCR40_PSACR.mn_ARL;fprintf('%.3f ',a(1,:));fprintf('\n');fprintf('%.3f ',a(2,:));
 save sstCR40_PSACR sstCR40_PSACR

%%
%%%%%%%%%%%  End of file HRCR_MSTATS.TXT
