fit_INST_ARL2_part1.m

File: work/CMUExper/fdbk1/data/fit_INST_ARL2_part1.m
Date: 2007-08-16
Alexander Petrov, http://alexpetrov.com

Contents

Load the data

clear all
cd(fullfile(work_pathstr,'CMUExper','fdbk1','data')) ;
load('S2.mat') ;
ARL2 = [S2(:).ARL2] ;    % 9x55
gr=[S2(:).group]' ;
ugr = unique(gr) ;
N_groups = length(ugr) ;
ugr_descr = {'U1' 'L1' 'H1', 'L2' 'H2' } ;   % [1 2 3, 5 6]

Collapse across groups

grARL2 = zeros(size(ARL2,1),N_groups) ;   % 9x5
gr_idx = cell(1,N_groups) ;
for k = 1:N_groups
    gr_idx{k} = find(gr==ugr(k))' ;
    grARL2(:,k) = mean(ARL2(:,gr_idx{k}),2) ;
end
grARL2

grARL2 =
    4.3358    4.1354    4.2075    4.1789    4.2197
    4.0525    4.0594    4.1666    4.6730    4.1990
    4.1878    3.9444    4.2666    4.7687    4.0890
    4.4236    4.4190    4.4370    3.9850    4.3336
    4.3422    4.2662    4.4461    3.9343    4.3488
    4.2404    3.8147    4.4136    4.1016    4.2765
    4.0962    3.7438    4.3909    4.2058    4.4644
    4.2742    3.9327    4.3790    3.8922    4.3864
    4.3494    3.6773    4.5247    3.7794    4.3898

Fit INST on the same stimulus sequences, 11 per group

stim = [S2(:).dist] ;   % 476x55
fdbk = [S2(:).fdbk] ;

Search parameters

These control the optimizer via the PARAMSEARCH Toolbox. FMINCON is used (instead of LSQNONLIN) because experimental fits of synthetic ARL2 profiles showed that it works somewhat better in our setting. Still, the cutoff parameters cannot be fitted with these methods, probably because the error derivative is not continuous w.r.t. the cutoffs. The plan is to explore the two cutoff parameters on a grid, while optimising w.r.t three others: mem_k, temper, and hist. The perceptual noise is fixed a priori to perc_k=.040. The history weight parameter is poorly determined by the data. Therefore, at first we fix it at H=.040 and do a grid search on the cutoff parameters. When we narrow down the relevant region of the parameter space, we will revisit H again.

fminconp = 1 ;
flags = [0 1 1 0 0 0 0 0] ;   % [perc_k mem_k temper hist cutoff xx xx c_ratio]
Sparams = protoanchor_search_params(flags,fminconp) ;
optns = optimset(Sparams.optns,'TolX',1e-2,'Display','off') ;
Sparams.optns = optns ;
Sparams.target_vals = grARL2
H = .040 ; Sparams.params.history = H ;   % for the time being
Sparams.params

Sparams = 
     model_name: 'protoanchor_ARL2'
         params: [1x1 struct]
      p2v_templ: {'VAL = PARAMS.mem_k * 10 ;'  'VAL = PARAMS.temper * 10 ;'}
      v2p_templ: {'PARAMS.mem_k = VAL / 10 ;'  'PARAMS.temper = VAL / 10 ;'}
         bounds: [2x2 double]
     ARL_params: [1x1 struct]
    funfun_name: 'fmincon'
          optns: [1x1 struct]
     lsfun_name: 'sumsqerr'
    target_vals: [9x5 double]
ans = 
        scale: 'LINEAR'
        N_cat: 7
      SM_conv: 1.0000e-03
       cat_sz: 0.0400
       perc_k: 0.0400
        mem_k: 0.0700
        avail: [7x1 double]
    instances: [7x3 double]
      cutoffs: [-2.4000 -0.8000 0.7200 2.1600]
       temper: 0.0400
      history: 0.0400
        alpha: 0.3000
        decay: 0.5000
          ITI: 4

"Chance event" sets

The model is fitted using the deterministic version PROTOANCHOR_CHEV that uses previously cached values of all randomly generated values. More than one such "chev" sets is used to evaluate its impact on the fit.

N_chevs = 2 ;

% Pack into "model data" structures suitable for passing as arguments.
Mdata = cell(1,N_chevs) ;
for k = 1:N_chevs
    Mdata{k}.stim = stim ;
    Mdata{k}.fdbk = fdbk ;
    Mdata{k}.group_idx = gr_idx ;
    %Mdata{k}.params = Sparams.params ;    % will be set by SUMSQERR
    Mdata{k}.chev = make_protoanchor_chev(size(stim,1),size(stim,2)) ;
    Mdata{k}.ARL_params = Sparams.ARL_params ;
end

Cutoff parameter grid

cutoff_levels = [.50 .55 .60 .65 .70 .80 .90 1.00]' ;  % c- in categ_sz units
N_cutoff_levels = length(cutoff_levels) ;
c_ratio_levels = [.5 .6 .7 .8 .9] ;          % c-/c+ ratio
N_c_ratio_levels = length(c_ratio_levels) ;

Main work

A grid search over 5x3x2 cells takes five hours on a fast machine. Therefore, the results will only be computed once and cached to disk.

fname = 'grid_search_INST_ARL2_part1.mat' ;
if (file_exists(fname))
    load(fname) ;
else
  grid_SSE = zeros([N_cutoff_levels,N_c_ratio_levels,N_chevs]) ;
  grid_optX = cell([N_cutoff_levels,N_c_ratio_levels,N_chevs]) ;
  grid_det = cell([N_cutoff_levels,N_c_ratio_levels,N_chevs]) ;
  tic
  for k1 = 1:N_cutoff_levels
      c = cutoff_levels(k1) ;
      fprintf('\n cutoff %4.2f: ',c) ;
      for k2 = 1:N_c_ratio_levels
          cr = c_ratio_levels(k2) ;
          cutoffs = c .* [-3 -1 cr 3*cr] ;
          Sparams.params.cutoffs = cutoffs ;
          for k = 1:N_chevs
              [par,SSE,optX,det] = paramsearch(Mdata{k},Sparams) ;  % <-- Sic!
              grid_SSE(k1,k2,k) = SSE ;
              grid_optX{k1,k2,k} = optX ;
              grid_det{k1,k2,k} = det ;
              fprintf('%5.3f ',SSE) ;
              save(fname,'grid_SSE','grid_optX','grid_det','Sparams',...
                   'cutoff_levels','c_ratio_levels') ;
          end
          fprintf(', ') ;
      end
  end
  fprintf('\n') ;
  toc
end % if file_exists ...

Plot the SSE

SSE = mean(grid_SSE,3)   % average across the chev's
[X,Y] = meshgrid(c_ratio_levels,cutoff_levels) ;
contourf(X,Y,SSE) ; colorbar ; title('SSE') ;
xlabel('c-/c+ ratio') ; ylabel('c- cutoff') ;

SSE =
    1.7888    1.5530    1.4217    1.6504    2.0318
    1.8492    1.5656    1.4132    1.5419    2.2204
    1.9386    1.5639    1.3889    1.5422    2.1584
    2.0268    1.6076    1.4028    1.7743    2.4702
    2.3314    1.4707    1.4331    1.8232    2.4112
    2.2935    1.4553    1.4324    1.9165    2.3614
    2.2761    1.5362    1.5492    1.8376    2.6434
    2.0776    1.9287    1.5657    2.1168    2.7561

Plot the memory noise parameter

mem_k = zeros(size(grid_optX)) ;
temper = zeros(size(grid_optX)) ;
for k1 = 1:N_cutoff_levels
    for k2 = 1:N_c_ratio_levels
        for k = 1:N_chevs
            mem_k(k1,k2,k) = grid_optX{k1,k2,k}(1) / 10 ;
            temper(k1,k2,k) = grid_optX{k1,k2,k}(2) / 10 ;
        end
    end
end

mem_k = mean(mem_k,3)
contourf(X,Y,mem_k) ; colorbar ; title('Memory noise') ;
xlabel('c-/c+ ratio') ; ylabel('c- cutoff') ;

mem_k =
    0.0200    0.0527    0.0542    0.0714    0.0696
    0.0200    0.0408    0.0464    0.0426    0.0664
    0.0200    0.0365    0.0489    0.0412    0.0643
    0.0200    0.0200    0.0345    0.0428    0.0250
    0.0200    0.0310    0.0340    0.0254    0.0326
    0.0273    0.0251    0.0346    0.0213    0.0368
    0.0200    0.0252    0.0201    0.0240    0.0202
    0.0200    0.0452    0.0209    0.0200    0.0200

Plot the temperature parameter

temper = mean(temper,3)
contourf(X,Y,temper) ; colorbar ; title('Softmax temperature') ;
xlabel('c-/c+ ratio') ; ylabel('c- cutoff') ;

temper =
    0.0800    0.0558    0.0266    0.0355    0.0329
    0.0800    0.0569    0.0246    0.0226    0.0368
    0.0800    0.0504    0.0317    0.0205    0.0384
    0.0744    0.0441    0.0319    0.0306    0.0416
    0.0567    0.0418    0.0462    0.0521    0.0429
    0.0240    0.0420    0.0280    0.0653    0.0406
    0.0200    0.0200    0.0281    0.0451    0.0456
    0.0200    0.0274    0.0281    0.0475    0.0438

Interim conclusion

The optimum seems to be located for: cutoff = .60 category_sz units cutoff_ratio = 0.7 Thus, cutoffs = .60 .* [-3 -1 .7 3*.7] = [-1.8 -.6 .42 1.26]. Memory noise seems to be about mem_k = .050 (or less), which is less than the default (.070) in PetrovAnderson05. The temperature parameter seems poorly determined, but temper = .040 seems a reasonable starting point for subsequent explorations. This is again less than the default (.050) in PetrovAnderson05.

optX = [.040 .050].*10 ;
Mparams = vector2params(optX,Sparams.params,Sparams.v2p_templ) ;
Mparams.cutoffs = .60 .* [-3 -1 .7 3*.7]

Mparams = 
        scale: 'LINEAR'
        N_cat: 7
      SM_conv: 1.0000e-03
       cat_sz: 0.0400
       perc_k: 0.0400
        mem_k: 0.0400
        avail: [7x1 double]
    instances: [7x3 double]
      cutoffs: [-1.8000 -0.6000 0.4200 1.2600]
       temper: 0.0500
      history: 0.0400
        alpha: 0.3000
        decay: 0.5000
          ITI: 4

Generate and ARL2 profiles with these parameters

for k = 1:N_chevs
    Mdata{k}.params = Mparams ;
    opt_grARL2(:,:,k) = protoanchor_ARL2(Mdata{k}) ;
end
opt_grARL2 = mean(opt_grARL2,3) ;
[grARL2 opt_grARL2]
opt_SSE = (grARL2-opt_grARL2) ;
opt_SSE = sum(opt_SSE(:).^2)

ans =
    4.3358    4.1354    4.2075    4.1789    4.2197    4.1094    4.0922    4.1151    4.1485    4.1193
    4.0525    4.0594    4.1666    4.6730    4.1990    4.1061    4.0228    4.2128    4.1561    4.2446
    4.1878    3.9444    4.2666    4.7687    4.0890    4.1031    3.9662    4.2516    4.2408    4.3463
    4.4236    4.4190    4.4370    3.9850    4.3336    4.1511    4.0063    4.3381    4.1448    4.3778
    4.3422    4.2662    4.4461    3.9343    4.3488    4.1914    4.0867    4.4002    4.0964    4.3611
    4.2404    3.8147    4.4136    4.1016    4.2765    4.1994    4.0561    4.3598    4.1028    4.3661
    4.0962    3.7438    4.3909    4.2058    4.4644    4.1582    4.0096    4.3427    4.1044    4.4449
    4.2742    3.9327    4.3790    3.8922    4.3864    4.1771    4.0186    4.3711    4.0743    4.4204
    4.3494    3.6773    4.5247    3.7794    4.3898    4.2143    4.0323    4.4161    4.0481    4.4310
opt_SSE =
    1.5061

Plot the "optimal" ARL2 profiles

ls = {'r-' 'b-' 'k<-' 'b--' 'k<--'} ;   % group line styles
lw = [1 2 2  2 2] ;                     % group line widths
els  = {'r:' 'b:' 'k<:' 'b:' 'k<:'} ;   % empirical-data line styles

x2 = [14 , 56:56:448]' ;
xtick = [28 , 140:112:476] ;
ytick = [3.6:.2:4.8] ;
ax = [0 476 3.6 4.8] ;

clf ; hold on ;
for k = 1:N_groups
    hh = plot(x2,opt_grARL2(:,k),ls{k}) ;
    set(hh,'LineWidth',lw(k),'MarkerSize',4) ;
end
for k = 1:N_groups
    plot(x2,grARL2(:,k),els{k}) ;
end
hold off ; box on ;
axis(ax) ; set(gca,'xtick',xtick,'ytick',ytick,'xgrid','on') ;
hh = xlabel('Trial') ; set(hh,'FontSize',14) ;
hh = ylabel('Average response level') ; set(hh,'FontSize',14) ;
hh = legend('U1','L1','H1','L2','H2',3) ; set(hh,'FontSize',14) ;
title(sprintf('INST fits, SSE = %.3f',opt_SSE)) ;

Define Assimilation := ARL(N)-ARL(P) and plot Assim profile

See .../CMUExper/fdbk1/data/empir_ARL2_profile.m

%Assim = grARL2(:,[3 find(ugr==6)]) - grARL2(:,[2 find(ugr==5)])
opt_Assim = opt_grARL2(:,[3 find(ugr==6)]) - opt_grARL2(:,[2 find(ugr==5)])

ax = [0 476 -0.85 1.15] ;
ytick1 = [-.8:.2:1] ;
idx1 = [1 2 3 6 7] ;    % feedback blocks for the feedback-first groups
idx2 = [1 4 5 8 9] ;    % feedback blocks for the no-feedback-first grps

clf ;
hh = plot(x2,opt_Assim(:,1),'r-',x2,opt_Assim(:,2),'b--',...
          x2(idx1),opt_Assim(idx1,1),'ro',x2(idx2),opt_Assim(idx2,2),'bs') ;
    set(hh,'LineWidth',2) ;
    axis(ax) ;
    hh=refline(0,0) ;
    set(hh,'Color','k','LineStyle','-','LineWidth',1) ;
box on ;
axis(ax) ; set(gca,'xtick',xtick,'ytick',ytick1,'xgrid','on') ;
hh = xlabel('Trial') ; set(hh,'FontSize',14) ;
hh = ylabel('ARL(H) - ARL(L)') ; set(hh,'FontSize',14) ;
hh = legend('Feedback first','No-feedback first',4) ; set(hh,'FontSize',14) ;

% The fitting effort continues in fit_INST_ARL2_part2.m

opt_Assim =
    0.0230   -0.0292
    0.1900    0.0885
    0.2854    0.1055
    0.3318    0.2329
    0.3135    0.2647
    0.3037    0.2633
    0.3331    0.3405
    0.3525    0.3461
    0.3838    0.3829

Clean up

clear k k1 k2 X Y ax hh ;

Published with MATLAB® 7.4