Home > voicebox > psycest.m

psycest

PURPOSE ^

Estimate multiple psychometric functions

SYNOPSIS ^

function [xx,ii,m,v,mr,vr]=psycest(iq,x,r,xp,lf)

DESCRIPTION ^

 Estimate multiple psychometric functions

 Usage: [xx,ii,m,v]=psycest(-n,p,q,xp,lf) % initialize n models
        [xx,ii,m,v]=psycest(i,x,r)     % supply a trial result to psycest
        [xx,ii,m,v]=psycest(i,x,r,o)   % supply a trial result to psycest and plot
                    psycest(i,o)       % plot pdf of model i with plot options o
        [xx,ii,m,v,mr,vr]=psycest(i,o) % Determine robust outputs in addition expcluding outliers [takes longer]
              [p,q,msr]=psychest(0)    % output model parameters (or print them if no outputs)

 Inputs (see usage examples for argument positions):
         -n        minus the number of models
          p(:,n)   parameters for each model
                      1  thresh [0.75]
                      2  miss [0.04]
                      3  guess [0.1]
                      4  SNR min [-20]
                      5  SNR max [20]
                      6  Slope min (but over-ridden if <q.sl) [0]
                      7  slope max [0.5]
          q(:)     parameters common to all models (vector or struct)
                      1  q.nx  number of SNR values in pdf [40]
                      2  q.ns  number of slope values in pdf [21]
                      3  q.nh  number of probe SNR values to evaluate [30]
                      4  q.cs  weighting of slope relative to SRT in cost function [1]
                      5  q.dh  minimum step size in dB for probe SNRs [0.2]
                      6  q.sl  min slope at threshold (must be >0) [0.005]
                      7  q.kp  number of std deviations of the pdfs to keep [4]
                      8  q.hg  amount to grow expected gains in ni trials [1.3]
                      9  q.cf  cost function: 1=variance, 2=entropy [2]
                     10  q.pm  psychometric model: 1=logistic, 2=cumulative gaussian [1]
                     11  q.lg  use log slope in pdf: 0=no, 1=yes [1]
                     12  q.pp  Number of prior standard deviations in initial semi-range [1]
                     13  q.pf  Probability floor (integrated over entire grid) [0.0001]
                     14  q.ts  Number of std devs to explore [2]
                     15  q.dp  Maximum probe SNR shift (dB) [10]
                     16  q.it  Grid interpolation threshold [0.5]
                     17  q.at  Axis change threshold [0.1]
                     18  q.la  Look 2-ahead when choosing probe [1]
                     19  q.op  Outlier probability [0.01]
                     20  q.rx  Minimum range factor per iteration [0.5]
          xp{n}(:) list of available probe SNR values
          lf       log file ID
          i        model probed
          x        probe SNR value used
          r        response: 0=wrong, 1=right.
          o        plot option string:
                   'p' Plot pdf
                   'h' Plot probe history
                   'x' Plot expected cost function for probe SNRs
                   'c' Plot cost function evolution

 Outputs:
          xx       recommended probe SNR
          ii       recommended model to probe next
          m(2,n,3) estimated srt and slope of all models
                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
          v(3,n)   estimated covariance matrix entries:
                   [var(srt) cov(srt,slope) var(slope)]'
          mr(2,n,3)robust estimated srt and slope of all models
                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
          vr(3,n)  robust estimated covariance matrix entries:
                   [var(srt) cov(srt,slope) var(slope)]'
          msr(:,3)  List probe snrs and results: [model probe-snr result]

 Algorithm parameters:

 The algorithm estimates the psychometric function for a number of models simultaneously.
 A psychometric function specifies the probability of correct recognition as a function of
 SNR and is completely specified by two free parameters: the SNR (in dB) and the slope (in 1/dB)
 at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify
 some of the details of the psychometric function and can be different for each model:
   p(1,n) gives the target recognition probability
   p(2,n) gives the probability of incorrect recognition at very good SNRs (the "miss" probability).
          If this value is made too low, a single unwarrented recognition error will have a large
          effect of the estimated parameters and greatly lengthen the time to convergence. The default
          value is 4%.
   p(3,n) gives the probability of correct recognition at very poor SNRs (the "guess" probabiity).
          This should be set to 1/N where N is the number of possible responses to a stimulus.
 p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will
           be adjusted by the algorithm if necessary but wrong values will delay convergence.
 p(6:7,n) These given the initial min and max slope (in probability per dB) at the target
          recognition probability.
 The remaining parameters are shared between all the models and control how the algorithm operates.
   q(1:2) These parameters determine the sampling size of the joint pdf in SNR and Slope respectively.
   q(3)   This parameter specifies how many potential values of probe SNR to evaluate in order to
          determine which is likely to give the most improvement in the parameter estimates.
   q(4)   This parameter gives the relative weight of SNR and Slope in the cost function. Increasing
          its value will improve the slope estimate (or log(slope) estimate) at the expense of the
          SNR estimate. Actually its value is not all that critical.
   q(5)   At each iteration psycest evaluates several potential probe SNR levels to see which is
          likely to give the most improvement to the parameter estimates. This parameter fixes the
          minimum spacing between these potential probe values. It only has an effect when the variances
          of the parameter estimates are very small.
   q(6)   To prevent the routine getting lost, this parameter specifies the smallest reasonable value
          of the Slope parameter at the target recognition probability. Under normal circumstances, it
          has no effect.
   q(7)   For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope.
          The sampling grid is reclculated after each iteration to encompass this number of standard
          deviations in each dimension.
   q(8)   At each iteration, psycest advises which model to probe next on the basis of which
          gives the greatest expected reduction in cost function. To ensure that all models
          are periodically sampled, this expected reduction of each unprobed model is multiplied
          by this parameter at each iteration. Thus a model will eventually be probed even if
          its expected cost factor improvement is small.
   q(9)   This determines whether the recommended probe SNR to use at the next iteration is chosen to
          minimize the expected variance or the entropy of the estimated SNR and slope distributions.
          My experience is that entropy (the default) gives faster convergence.
  q(10)   This selects whether the underlying model is a logistic function (1) or a cumulative
          gaussian (2). The choice makes little difference unless the "miss" or "guess" probabilities
          are very very small.
  q(11)   This selects whether the internal pdf samples "slope" (0) or "log(slope)" (1). It is recommended
          that you stick to the default of log(slope) since this results in more symmetrical distributions. [1]
  q(12)   This sets the number of std deviations of the prior corresponding to the limits of the
          initial distribution set in p(4:7,n). A small value corresponds to a very weak prior. [1]
  q(13)   Probability floor (integrated over entire grid) [0.0001]
  q(14)   Number of std devs to explore when evaluating possible probe SNRs [2]
  q(15)   Maximum probe SNR shift beyond the range of previous probes (dB) [10]
  q(16)   Grid interpolation threshold. The grid is only interpolated if
          it has shrunk by less that this factor [0.5]
  q(17)   If the desired grid limits change by less than this fraction of
          the grid length, they will be left unaltered. [0.1]
  q(18)   Look 2-ahead when choosing probe [1]
  q(19)   Outlier probability; probe results with lower probability than
          this will be ignored [0.01]
  q(20)   The ranges of SRT and log slope will not be reduced per iteration below this factor [0.5]

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [xx,ii,m,v,mr,vr]=psycest(iq,x,r,xp,lf)
0002 % Estimate multiple psychometric functions
0003 %
0004 % Usage: [xx,ii,m,v]=psycest(-n,p,q,xp,lf) % initialize n models
0005 %        [xx,ii,m,v]=psycest(i,x,r)     % supply a trial result to psycest
0006 %        [xx,ii,m,v]=psycest(i,x,r,o)   % supply a trial result to psycest and plot
0007 %                    psycest(i,o)       % plot pdf of model i with plot options o
0008 %        [xx,ii,m,v,mr,vr]=psycest(i,o) % Determine robust outputs in addition expcluding outliers [takes longer]
0009 %              [p,q,msr]=psychest(0)    % output model parameters (or print them if no outputs)
0010 %
0011 % Inputs (see usage examples for argument positions):
0012 %         -n        minus the number of models
0013 %          p(:,n)   parameters for each model
0014 %                      1  thresh [0.75]
0015 %                      2  miss [0.04]
0016 %                      3  guess [0.1]
0017 %                      4  SNR min [-20]
0018 %                      5  SNR max [20]
0019 %                      6  Slope min (but over-ridden if <q.sl) [0]
0020 %                      7  slope max [0.5]
0021 %          q(:)     parameters common to all models (vector or struct)
0022 %                      1  q.nx  number of SNR values in pdf [40]
0023 %                      2  q.ns  number of slope values in pdf [21]
0024 %                      3  q.nh  number of probe SNR values to evaluate [30]
0025 %                      4  q.cs  weighting of slope relative to SRT in cost function [1]
0026 %                      5  q.dh  minimum step size in dB for probe SNRs [0.2]
0027 %                      6  q.sl  min slope at threshold (must be >0) [0.005]
0028 %                      7  q.kp  number of std deviations of the pdfs to keep [4]
0029 %                      8  q.hg  amount to grow expected gains in ni trials [1.3]
0030 %                      9  q.cf  cost function: 1=variance, 2=entropy [2]
0031 %                     10  q.pm  psychometric model: 1=logistic, 2=cumulative gaussian [1]
0032 %                     11  q.lg  use log slope in pdf: 0=no, 1=yes [1]
0033 %                     12  q.pp  Number of prior standard deviations in initial semi-range [1]
0034 %                     13  q.pf  Probability floor (integrated over entire grid) [0.0001]
0035 %                     14  q.ts  Number of std devs to explore [2]
0036 %                     15  q.dp  Maximum probe SNR shift (dB) [10]
0037 %                     16  q.it  Grid interpolation threshold [0.5]
0038 %                     17  q.at  Axis change threshold [0.1]
0039 %                     18  q.la  Look 2-ahead when choosing probe [1]
0040 %                     19  q.op  Outlier probability [0.01]
0041 %                     20  q.rx  Minimum range factor per iteration [0.5]
0042 %          xp{n}(:) list of available probe SNR values
0043 %          lf       log file ID
0044 %          i        model probed
0045 %          x        probe SNR value used
0046 %          r        response: 0=wrong, 1=right.
0047 %          o        plot option string:
0048 %                   'p' Plot pdf
0049 %                   'h' Plot probe history
0050 %                   'x' Plot expected cost function for probe SNRs
0051 %                   'c' Plot cost function evolution
0052 %
0053 % Outputs:
0054 %          xx       recommended probe SNR
0055 %          ii       recommended model to probe next
0056 %          m(2,n,3) estimated srt and slope of all models
0057 %                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
0058 %          v(3,n)   estimated covariance matrix entries:
0059 %                   [var(srt) cov(srt,slope) var(slope)]'
0060 %          mr(2,n,3)robust estimated srt and slope of all models
0061 %                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
0062 %          vr(3,n)  robust estimated covariance matrix entries:
0063 %                   [var(srt) cov(srt,slope) var(slope)]'
0064 %          msr(:,3)  List probe snrs and results: [model probe-snr result]
0065 %
0066 % Algorithm parameters:
0067 %
0068 % The algorithm estimates the psychometric function for a number of models simultaneously.
0069 % A psychometric function specifies the probability of correct recognition as a function of
0070 % SNR and is completely specified by two free parameters: the SNR (in dB) and the slope (in 1/dB)
0071 % at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify
0072 % some of the details of the psychometric function and can be different for each model:
0073 %   p(1,n) gives the target recognition probability
0074 %   p(2,n) gives the probability of incorrect recognition at very good SNRs (the "miss" probability).
0075 %          If this value is made too low, a single unwarrented recognition error will have a large
0076 %          effect of the estimated parameters and greatly lengthen the time to convergence. The default
0077 %          value is 4%.
0078 %   p(3,n) gives the probability of correct recognition at very poor SNRs (the "guess" probabiity).
0079 %          This should be set to 1/N where N is the number of possible responses to a stimulus.
0080 % p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will
0081 %           be adjusted by the algorithm if necessary but wrong values will delay convergence.
0082 % p(6:7,n) These given the initial min and max slope (in probability per dB) at the target
0083 %          recognition probability.
0084 % The remaining parameters are shared between all the models and control how the algorithm operates.
0085 %   q(1:2) These parameters determine the sampling size of the joint pdf in SNR and Slope respectively.
0086 %   q(3)   This parameter specifies how many potential values of probe SNR to evaluate in order to
0087 %          determine which is likely to give the most improvement in the parameter estimates.
0088 %   q(4)   This parameter gives the relative weight of SNR and Slope in the cost function. Increasing
0089 %          its value will improve the slope estimate (or log(slope) estimate) at the expense of the
0090 %          SNR estimate. Actually its value is not all that critical.
0091 %   q(5)   At each iteration psycest evaluates several potential probe SNR levels to see which is
0092 %          likely to give the most improvement to the parameter estimates. This parameter fixes the
0093 %          minimum spacing between these potential probe values. It only has an effect when the variances
0094 %          of the parameter estimates are very small.
0095 %   q(6)   To prevent the routine getting lost, this parameter specifies the smallest reasonable value
0096 %          of the Slope parameter at the target recognition probability. Under normal circumstances, it
0097 %          has no effect.
0098 %   q(7)   For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope.
0099 %          The sampling grid is reclculated after each iteration to encompass this number of standard
0100 %          deviations in each dimension.
0101 %   q(8)   At each iteration, psycest advises which model to probe next on the basis of which
0102 %          gives the greatest expected reduction in cost function. To ensure that all models
0103 %          are periodically sampled, this expected reduction of each unprobed model is multiplied
0104 %          by this parameter at each iteration. Thus a model will eventually be probed even if
0105 %          its expected cost factor improvement is small.
0106 %   q(9)   This determines whether the recommended probe SNR to use at the next iteration is chosen to
0107 %          minimize the expected variance or the entropy of the estimated SNR and slope distributions.
0108 %          My experience is that entropy (the default) gives faster convergence.
0109 %  q(10)   This selects whether the underlying model is a logistic function (1) or a cumulative
0110 %          gaussian (2). The choice makes little difference unless the "miss" or "guess" probabilities
0111 %          are very very small.
0112 %  q(11)   This selects whether the internal pdf samples "slope" (0) or "log(slope)" (1). It is recommended
0113 %          that you stick to the default of log(slope) since this results in more symmetrical distributions. [1]
0114 %  q(12)   This sets the number of std deviations of the prior corresponding to the limits of the
0115 %          initial distribution set in p(4:7,n). A small value corresponds to a very weak prior. [1]
0116 %  q(13)   Probability floor (integrated over entire grid) [0.0001]
0117 %  q(14)   Number of std devs to explore when evaluating possible probe SNRs [2]
0118 %  q(15)   Maximum probe SNR shift beyond the range of previous probes (dB) [10]
0119 %  q(16)   Grid interpolation threshold. The grid is only interpolated if
0120 %          it has shrunk by less that this factor [0.5]
0121 %  q(17)   If the desired grid limits change by less than this fraction of
0122 %          the grid length, they will be left unaltered. [0.1]
0123 %  q(18)   Look 2-ahead when choosing probe [1]
0124 %  q(19)   Outlier probability; probe results with lower probability than
0125 %          this will be ignored [0.01]
0126 %  q(20)   The ranges of SRT and log slope will not be reduced per iteration below this factor [0.5]
0127 
0128 %      Copyright (C) Mike Brookes and Clement Doire 2009-2016
0129 %      Version: $Id: psycest.m 9215 2016-12-17 22:37:49Z dmb $
0130 %
0131 %   VOICEBOX is a MATLAB toolbox for speech processing.
0132 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0133 %
0134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0135 %   This program is free software; you can redistribute it and/or modify
0136 %   it under the terms of the GNU General Public License as published by
0137 %   the Free Software Foundation; either version 2 of the License, or
0138 %   (at your option) any later version.
0139 %
0140 %   This program is distributed in the hope that it will be useful,
0141 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0142 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0143 %   GNU General Public License for more details.
0144 %
0145 %   You can obtain a copy of the GNU General Public License from
0146 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0147 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0149 
0150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0151 % Bugs/Suggestions:
0152 % (2)use a structure for input parameters
0153 % (3) entropy seems to change signifiacntly when the grid is rescaled
0154 % (4)add a forgetting factor for old measurements
0155 % (8) use quadratic interpolation to find cost function minimum
0156 % (10) optionally output the whole pdf + its axis values
0157 % (11) optionally output all model probe values
0158 % (13) Should perhaps estimate and return the mean and slope compensated for the guess rate
0159 %      i.e. if you want p, you actually test at guess+ p*(1-guess-miss)/(1-miss)
0160 % (17) remember probe snrs, models and results and recalculate the entire
0161 %      array when changing precision
0162 % (18) could force probe SNRs to be a multiple of minimum step size
0163 % (19) use a non-uniform prior e.e. 1/(1+x^2)
0164 % (20) possibly use different parametrization (e.g. 1/slope or a second SRT threshold)
0165 % (22) save inputs so that pdfs can be recalculated when rescaling
0166 % (24) Parametrize slope as x=(100s^2-1)/20s to make the distribution more uniform
0167 %      inverse is s=(sqrt(x^2+1)-x)/10; alternatively use log(slope)
0168 % (25) optionally have no prior (i.e. maximum likelihood)
0169 % (26) Check why the estimated variances are too small
0170 % (27) Adapt range of potential probes if optimum is at the limit
0171 % (28) Resample the pdf on the basis of the marginal distributions
0172 % (29) Select the probe range on the basis of marginal distributions
0173 % (30) Resample the pdfs by a factor of 2 always
0174 % (31) use uneven samples in pdf concentrated in the middle
0175 % (32) Use a parametric mixture including one-sided constants for the pdf
0176 % (33) Selection of subset of prescribed probe SNRs is not optimum
0177 % (34) use quadratic interpolation to select the probe SNR unless using fixed values
0178 % (35) expand the pdf grid based on the effective number of samples (like particle filters)
0179 
0180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0181 % persistent variables
0182 % sizes: nx=SNR values in pdf
0183 %        ns=slope values in pdf (1/std dev of cumulative gaussian)
0184 %        nh=potential test SNRs
0185 %        ni=simultaneous models being estimated
0186 % wq(ns*nx,ni)  = log prob of model; each column is a vectorized ns*nx matrix
0187 % nr(1,10)      = parameter values: [SNR-pdf slope-pdf SNR-probe ...]
0188 % pr(7,ni)      = input model parameters
0189 % qr(4,ni)      = derived model parameters
0190 % xq(nr(1),ni)  = SNR values in pdf
0191 % sq(nr(2),ni)  = slope values in pdf
0192 % mq(2,ni,3)    = estimated srt and slope of all models
0193 % vq(3,ni)      = estimated covariance matrix entries:[var(srt) cov(srt,slope) var(slope)]'
0194 % xn(1,ni)      = next probe value to use
0195 % hn(1,ni)      = expected decrease in cost function after next probe
0196 % hfact
0197 % xz
0198 % res(nres,7)   = results [model-number, probe-SNR, result]
0199 % nres          = total number of probe results
0200 % nresq(1,ni)   = number of probe results for each model
0201 % xmm(2,ni)     = [min(probe-SNR); max(probe-SNR)]
0202 % mq0(2,ni)      = prior means [x-mean; s-mean]
0203 % pq0(2,ni)      = scale factors for prior distribution
0204 % wfl           = floor relative to peak of log-pdf array
0205 % sqmin         = minimum value of slope or log-slope
0206 % LOG           = file ID of logging file
0207 % mqr(2,ni,3)   = robust mean estimates (excluding outliers)
0208 % vqr(3,ni)     = robust variance estimates (excluding outliers)
0209 % nresr         = nres value at which robust estimates were last calculated
0210 
0211 persistent wq xq sq nr pr qr mq vq xn hn hfact xz res nres nresq xmm mq0 pq0 wfl sqmin LOG mqr vqr nresr xlim
0212 
0213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0214 %              Initialization (iq<0)            %
0215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0216 if iq<0  % initialization
0217     if nargin>4 && ~isempty(lf)
0218         LOG=lf;
0219         fprintf(LOG,'******************************\npsycest Initialization: %s\n******************************\n',datestr(now));
0220     else
0221         LOG=0;
0222     end
0223     ni=-iq;                 % number of models
0224     nres=0;                 % total number of probe-results so far
0225     nresq=zeros(1,ni);      % number of probe-results for each model
0226     res=zeros(1,7);         % array for saving results [expands as needed]
0227     pr=repmat([0.75 0.04 0.1 -20 20 0 0.5]',1,ni);          % default parameters
0228     if nargin>1
0229         if size(x,2)>1
0230             if size(x,2)>ni
0231                 error('initialization parameter argument has too many columns');
0232             end
0233             pr(1:size(x,1),1:size(x,2))=x;
0234         else
0235             pr(1:size(x,1),:)=repmat(x,1,ni);
0236         end
0237     end
0238     nr=[40 21 30 1 0.2 0.02 4 1.3 2 1 1 1 0.0001 2 10 0.5 0.1 1 0.01 0.5]';                      % default parameter values
0239     nrf={'nx','ns','nh','cs','dh','sl','kp','hg','cf','pm','lg','pp', 'pf', 'ts', 'dp', 'it','at','la','op','rx'};     % parameter field names
0240     numnr=length(nr);
0241     if nargin>2 && ~isempty(r) % replace defaults with any parameter values specified in the call
0242         if isstruct(r)
0243             fnn=fieldnames(r);
0244             for i=1:length(fnn)
0245                 mk=strcmp(fnn{i},nrf);
0246                 if any(mk)
0247                     nr(mk)=r.(fnn{i});
0248                 end
0249             end
0250         else
0251             nr(1:min(numnr,numel(r)))=r(:); % if parameters are specified as a vector, copy them across
0252         end
0253         nr(1:3)=round(nr(1:3));             % first three parameters must be integers
0254     end
0255     pr(6,:)=max(pr(6,:),nr(6));                % low limit of slope in prob/dB
0256     nxq=nr(1);
0257     nsq=nr(2);
0258     nsxq=nxq*nsq;
0259     xq=(0:nxq-1)'*(pr(5,:)-pr(4,:))/(nxq-1)+repmat(pr(4,:),nxq,1);
0260     if nr(11)  % if log slope
0261         sqmin=log(nr(6)); % low limit for axis
0262         sq=(1-nsq:0)'*(log(pr(7,:))-log(max(pr(6,:),nr(6))))/(nsq-1)+repmat(log(pr(7,:)),nsq,1);
0263     else
0264         sqmin=nr(6); % low limit for axis
0265         sq=(1-nsq:0)'*(pr(7,:)-max(pr(6,:),nr(6)))/(nsq-1)+repmat(pr(7,:),nsq,1);
0266     end
0267     mq0=[mean(xq,1); mean(sq,1)];     % find mean of prior
0268     pq0=-2*nr(12)^2*[xq(end,:)-xq(1,:); sq(end,:)-sq(1,:)].^(-2); % scale factor for prior distribution
0269     wq=repmat(pq0(2,1)*(sq(:,1)-mq0(2,1)).^2,1,nxq)+repmat(pq0(1,1)*(xq(:,1)-mq0(1,1))'.^2,nsq,1); % log probs of prior (same for all models)
0270     wfl=log(nr(13)/nsxq); % floor relative to peak of wq array
0271     wq=repmat(wq(:),1,ni);              % initialize to +-q.pp std deviations of gaussian prior
0272     qr=zeros(5,ni);
0273     qr(1,:)=1-pr(2,:)-pr(3,:);             % prob range covered by cumulative gaussian
0274     if any(qr(1,:)<=0)
0275         [dum,i]=min(qr(1,:));
0276         error('Model %d: guess prob (%.2f) + miss prob (%.2f) is >=1',i,pr(3,i),pr(2,:));
0277     end
0278     qr(2,:)=(pr(1,:)-pr(3,:))./qr(1,:);                 % cumulative gaussian probability at threshold
0279     if any(abs(qr(2,:)-0.5)>=0.5)
0280         [dum,i]=max(qr(2,:)-0.5);
0281         error('Model %d: target SRT threshold (%.2f) must lie in the range guess prob (%.2f) to (1-miss prob) (%.2f)',i,pr(1,i),pr(3,i),1-pr(2,:));
0282     end
0283     switch nr(10)                                       % switch on model type
0284         case 1                                          % logistic model
0285             qr(3,:)=log(qr(2,:)./(1-qr(2,:)));            % x position of target in std measure
0286             qr(4,:)=qr(1,:).*qr(2,:).*(1-qr(2,:));        % slope*"stddev" at threshold
0287         case 2                                          % cumulative gaussian model
0288             qr(3,:)=norminv(qr(2,:));                   % x position of target in std measure
0289             qr(4,:)=qr(1,:).*normpdf(qr(3,:));           % slope*stddev at threshold
0290         otherwise
0291             error('Unrecognised psychometric model selection');
0292     end
0293     if nr(9)< 1 || nr(9)>2
0294         error('Unrecognised cost function option');
0295     end
0296     mq=repmat(mq0,[1,1,3]);        % initial means, joint mode and marginal mode all are equal
0297     vq=[var(xq,1,1); zeros(1,ni); var(sq,1,1)];         % initial variances (actually ignores prior probs)
0298     mqr=mq;         % robust means and modes
0299     vqr=vq;         % robust variances
0300     nresr=0;        % robust calculation time
0301     xlim=[-Inf(1,ni); Inf(1,ni)]; % SNR limits for outliers
0302     %     hn=[1-nr(4) 0 nr(4)]*vq;  % artificially high value of cost function ensures all models are probed early on
0303     hn=Inf(1,ni);           % very high initial cost function
0304     hfact=nr(8)^(1/ni);     % growth factor to ensure that no model is neglected for too long
0305     xn=mq0(1,:);             % start at mean value
0306     xmm=[xn; xn];  % min/max of probe values for each model
0307     if nargin>=4 && ~isempty(xp)
0308         if iscell(xp)
0309             xz=xp;
0310         else
0311             xz=repmat(num2cell(xp(:)',2),ni,1);        % else replicate for each model
0312         end
0313         for i=1:ni
0314             [dum,j]=min(abs(xz{i}-mq0(1,i)));      % select the closest available probe to the mean
0315             xn(i)=xz{i}(j);
0316         end
0317     else
0318         xz=cell(ni,1);          % make empty cells
0319     end
0320 elseif iq>0
0321     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0322     %     iq>0 means update or plot model iq           %
0323     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0324     %        first check the plotting options          %
0325     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0326     if nargin<3
0327         if nargin==2
0328             po=x; % plot options are 2nd argument
0329         elseif ~nargout
0330             po='p';
0331         else
0332             po='';
0333         end
0334     else
0335         if nargin>3
0336             po=xp; % plot options are 4th argument
0337         else
0338             po='';
0339         end
0340     end
0341     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0342     %     now get parameters of current model (iq)     %
0343     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0344     nxq=nr(1);                  % number of SNR values
0345     nsq=nr(2);                  % number of slope values
0346     nxh=nr(3);                  % number of probe-SNR values
0347     nsxq=nxq*nsq;               % size of log pdf array
0348     thresh=pr(1,iq);            % target probability at threshold
0349     guess=pr(3,iq);             % guess rate (1/choices)
0350     pscale=qr(1,iq);            % prob range left after subtracting miss and guess probs
0351     xtstd=qr(3,iq);             % x position of target in std measure
0352     xqi=xq(:,iq);               % SRT values of the pdf array
0353     sqi=sq(:,iq);               % slope (or log slope) values in PDF
0354     wqi=wq(:,iq);                  % log pdf array
0355     mqi= mq(:,iq,1);            % [xe; se] means
0356     vqi=vq(:,iq);               % [xv; sxv; sv]  covariance matrix
0357     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0358     %     rescale the pdfs if necessary         %
0359     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0360     ssd=sqrt(vqi(3));           % std deviation of slope
0361     xsd=sqrt(vqi(1));           % std deviation of SRT
0362     % determine range of new grid
0363     xqrange=xqi(end)-xqi(1);    % range of old SRT grid
0364     sqrange=sqi(end)-sqi(1);    % range of old slope grid
0365     if  nresq(iq)<2              % keep old limits if nresq(iq)<2
0366         xq2lim=xqi([1 end]);
0367         sq2lim=sqi([1 end]);
0368     else
0369         xqsemirange=max(nr(7)*xsd,0.5*nr(20)*xqrange);
0370         xq2lim=[mqi(1)-xqsemirange, mqi(1)+xqsemirange];
0371         sqsemirange=max(nr(7)*ssd,0.5*nr(20)*sqrange);
0372         sq2lim=[max(sqmin,mqi(2)-sqsemirange),mqi(2)+sqsemirange];
0373         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0374         % To avoid unnecessary changes, keep old limits if the change  %
0375         % is less than nr(17) times the previous range [nr(17)=0.1]    %
0376         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0377         if abs(xq2lim(1)-xqi(1))<nr(17)*xqrange
0378             xq2lim(1)=xqi(1);
0379         end
0380         if abs(xq2lim(2)-xqi(end))<nr(17)*xqrange
0381             xq2lim(2)=xqi(end);
0382         end
0383         
0384         if abs(sq2lim(1)-sqi(1))<nr(17)*sqrange
0385             sq2lim(1)=sqi(1);
0386         end
0387         if abs(sq2lim(2)-sqi(end))<nr(17)*sqrange
0388             sq2lim(2)=sqi(end);
0389         end
0390     end
0391     xq2=linspace(xq2lim(1),xq2lim(2),nxq)';   % new x axis values
0392     sq2=linspace(sq2lim(1),sq2lim(2),nsq)';   % new s axis values
0393     wqup=2;                                   % update flag
0394     if xq2(1)<xqi(1) || xq2(end)>xqi(end) || sq2(1)<sqi(1) || sq2(end)>sqi(end)
0395         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0396         % if extrapolating, recalculate log-pdfs from saved data  %
0397         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0398         if LOG, fprintf(LOG,'N=%d:%d, extrapolate: x(%.2g,%.2g)->(%.2g,%.2g) and s(%.2g,%.2g)->(%.2g,%.2g)\n',iq,nresq(iq),xqi([1 end]),xq2([1 end]),sqi([1 end]),sq2([1 end])); end
0399         wq2=repmat(pq0(2,iq)*(sq2-mq0(2,iq)).^2,1,nxq)+repmat(pq0(1,iq)*(xq2-mq0(1,iq))'.^2,nsq,1); % log probs of prior for model iq
0400         ires=find(res(1:nres,1)==iq); % find the results for this model
0401         if nr(11)                   % if log slope
0402             sqis=exp(sq2)/qr(4,iq);    % inverse std dev of gaussian (proportional to slope)
0403         else
0404             sqis=sq2/qr(4,iq);         % inverse std dev of gaussian (proportional to slope)
0405         end
0406         switch nr(10)   % switch on model type
0407             case 1
0408                 for i=1:length(ires)
0409                     j=ires(i);
0410                     r0=res(j,3)==0;
0411                     wq2=wq2+log(r0+(1-2*r0)*(guess+pscale*((1+exp(sqis*xq2'-xtstd-repmat(sqis,1,nxq)*res(j,2))).^(-1)))); %  P(l | r,x)
0412                 end
0413             case 2
0414                 for i=1:length(ires)
0415                     j=ires(i);
0416                     r0=res(j,3)==0;
0417                     wq2=wq2+log(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,1,nxq)*res(j,2)-sqis*xq2'+xtstd))); %  P(l | r,x)
0418                 end
0419         end
0420     else                                                    % possibly do quadratic interpolation
0421         wq2=max(reshape(wqi,nsq,nxq)-max(wqi(:)),wfl);        % turn into a normalized, clipped matrix for easy interpolation
0422         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0423         % use quadratic interpolation in SRT axis %
0424         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0425         if (xq2(end)-xq2(1))/(xqrange)>nr(16)          % if range has shrunk by < nr(16), leave the SRT axis unchanged
0426             xq2=xqi;                                  % copy the old SRT axis
0427             wqup=1;                                   % update flag
0428         else                                                % else do quadratic interpolation
0429             if LOG, fprintf(LOG,'N=%d:%d, interpolate x: x(%.2g,%.2g)->(%.2g,%.2g)\n',iq,nresq(iq),xqi([1 end]),xq2([1 end])); end
0430             xqf=(xq2-xqi(1))/(xqi(2)-xqi(1));
0431             xqj=ceil(xqf);
0432             xqf=xqj-xqf;
0433             xqg=1-xqf;
0434             xqh=0.25*xqf.*xqg; % Quadratic coefficient
0435             xqf((xqj<=0) | (xqj>nxq))=0;
0436             xqg((xqj<0) | (xqj>=nxq))=0;
0437             xqh((xqj<1) | (xqj>=nxq-1))=0;
0438             wq2=wq2(:,min(max(xqj,1),nxq)).*repmat(xqf'+xqh',nsq,1)+wq2(:,min(max(xqj+1,1),nxq)).*repmat(xqg'+xqh',nsq,1) - ...
0439                 (wq2(:,min(max(xqj-1,1),nxq))+wq2(:,min(max(xqj+2,1),nxq))).*repmat(xqh',nsq,1);
0440         end
0441         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0442         % use quadratic interpolation in slope axis %
0443         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0444         if (sq2(end)-sq2(1))/(sqrange)>nr(16)       % if range has shrunk by < nr(16), leave the slope axis unchanged
0445             sq2=sqi;                                % copy the old slope axis
0446             wqup=wqup-1;                            % update flag
0447         else
0448             if LOG, fprintf(LOG,'N=%d:%d, interpolate s: s(%.2g,%.2g)->(%.2g,%.2g)\n',iq,nresq(iq),sqi([1 end]),sq2([1 end])); end
0449             sqf=(sq2-sqi(1))/(sqi(2)-sqi(1));
0450             sqj=ceil(sqf);
0451             sqf=sqj-sqf;
0452             sqg=1-sqf;
0453             sqh=0.25*sqf.*sqg; % Quadratic coefficient
0454             sqf((sqj<=0) | (sqj>nsq))=0;
0455             sqg((sqj<0) | (sqj>=nsq))=0;
0456             sqh((sqj<1) | (sqj>=nsq-1))=0;
0457             wq2=wq2(min(max(sqj,1),nsq),:).*repmat(sqf+sqh,1,nxq)+wq2(min(max(sqj+1,1),nsq),:).*repmat(sqg+sqh,1,nxq) - ...
0458                 (wq2(min(max(sqj-1,1),nsq),:)+wq2(min(max(sqj+2,1),nsq),:)).*repmat(sqh,1,nxq);
0459         end
0460     end
0461     if wqup>0                               % check if we have updated either axis
0462         wq2=max(wq2(:)-max(wq2(:)),wfl);    % turn back into a normalized, clipped vector
0463         sqi=sq2;                            % update slope (or log slope) values in PDF
0464         xqi=xq2;                            % update SRT values of the pdf array
0465         wqi=wq2;                            % log pdf array
0466         sq(:,iq)=sqi;                       % save new s axis (log slope or slope)
0467         xq(:,iq)=xqi;                       % save new x axis (SRT)
0468         wq(:,iq)=wqi;                       % save log-pdf
0469     end
0470     if nr(11)                   % if log slope
0471         sqis=exp(sqi)/qr(4,iq);    % inverse std dev of gaussian (proportional to slope)
0472     else
0473         sqis=sqi/qr(4,iq);         % inverse std dev of gaussian (proportional to slope)
0474     end
0475     if nargin>=3
0476         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0477         % update pdfs with a new probe result  %
0478         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0479         if nres>=size(res,1)        % ensure there is space for a new result
0480             res=[res; zeros(nres,7)]; % double the size of the results array
0481         end
0482         nres=nres+1;                % increment the total number of results
0483         nresq(iq)=nresq(iq)+1;        % increment the number of results for this model
0484         res(nres,1:3)=[iq,x,r];       % save the latest result
0485         xmm(:,iq)=[min(x,xmm(1,iq)); max(x,xmm(2,iq))]; % update min/max probe values
0486         hn=hn*hfact;         % increase entropy gains to ensure all models get a chance
0487         
0488         %     disp(res(1:nres,:));  % ********** diagnostic
0489         
0490         %
0491         % update log probabilities with the previous test result
0492         % If r==1, we multiply (add in log) by a horizontally reflected version of the psychometric function that equals
0493         % the threshold (e.g. 0.5 or 0.75) at the probe SNR, x, for all slopes.
0494         % If r==0, we multiply (add in log) by 1- this reflected version which therefore equals (1-thresh) at x.
0495         %
0496         r0=r==0;
0497         switch nr(10)   % switch on model type
0498             case 1
0499                 wqi=wqi+log(r0+(1-2*r0)*(guess+pscale*((1+exp(reshape(sqis*xqi'-xtstd,nsxq,1)-repmat(sqis,nxq,1)*x)).^(-1)))); %  P(l | r,x)
0500             case 2
0501                 wqi=wqi+log(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,nxq,1)*x-reshape(sqis*xqi'-xtstd,nsxq,1)))); %  P(l | r,x)
0502             otherwise
0503                 error('Unrecognised psychometric model selection');
0504         end
0505         wq(:,iq)=wqi;             % save updated probabilities
0506     end
0507     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0508     % Calculate mean and covariance and entropy  %
0509     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0510     ewqi=exp(wqi-max(wqi));             % unnormalized probability vector
0511     ewqi=ewqi/sum(ewqi);
0512     wqsx=reshape(ewqi,nsq,nxq);         % normalized probabilities
0513     px=sum(wqsx,1);                             %  p(x0)
0514     ps=sum(wqsx,2);                             %  p(s0)
0515     xe=px*xqi;                                  % E(x0)
0516     se=ps'*sqi;                                 % E(s0)
0517     
0518     [pxpk,xm]=max(px);                          % marginal mode in x
0519     if xm>1 && xm<nxq                           % use quadratic interpolation in log prob if possible
0520         [dum,xm2]=quadpeak(log(px(xm-1:xm+1))');
0521         xm=xm+xm2-2;
0522     end
0523     xm=(2-xm)*xqi(1)+(xm-1)*xqi(2);             % marginal mode(x0)
0524     
0525     [pspk,sm]=max(ps);
0526     if sm>1 && sm<nsq                           % use quadratic interpolation in log prob if possible
0527         [dum,sm2]=quadpeak(log(ps(sm-1:sm+1)));
0528         sm=sm+sm2-2;
0529     end
0530     sm=(2-sm)*sqi(1)+(sm-1)*sqi(2);             % marginal mode(s0)
0531     
0532     [wqpk,j]=max(ewqi);
0533     i=1+floor((j-1)/nsq);
0534     j=j-nsq*(i-1);
0535     if i>1 && i<nxq && j>1 && j<nsq             % use quadratic interpolation in log prob if possible
0536         [dum,ji]=quadpeak(wqi(repmat((j-1:j+1)',1,3)+repmat(nsq*(i-2:i),3,1)));
0537         i=i+ji(2)-2;
0538         j=j+ji(1)-2;
0539     end
0540     xj=(2-i)*xqi(1)+(i-1)*xqi(2);               % joint mode  x
0541     sj=(2-j)*sqi(1)+(j-1)*sqi(2);               % joint mode: s
0542     xv=px*(xqi.^2)-xe^2;                        % Var(x0)
0543     sv=ps'*(sqi.^2)-se^2;                       % Var(s0)
0544     sxv=ewqi'*(repmat(sqi-se,nxq,1).*reshape(repmat(xqi'-xe,nsq,1),nsxq,1)); % Cov(s0*x0)
0545     mq(:,iq,1)=[xe; se];                        % save means
0546     mq(:,iq,2)=[xj; sj];                        % save joint mode
0547     mq(:,iq,3)=[xm; sm];                        % save marginal modes
0548     vq(:,iq)=[xv; sxv; sv];                     % save covariance matrix
0549     xh=(px*log(px)')*(xqi(1)-xqi(2));           % differential entropy h(x0)
0550     sh=(ps'*log(ps))*(sqi(1)-sqi(2));              % differential entropy h(s0)
0551     if nargin==3
0552         switch nr(9) % cost function: 1=variance, 2=entropy
0553             case 1
0554                 res(nres,4:7)=[xe se xv sv];               % save info for plotting history
0555             case 2
0556                 res(nres,4:7)=[xe se xh sh];                                           % find the minimum of cost function
0557         end
0558     end
0559     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0560     % now calculate robust mean and covariance for all models   %
0561     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0562     if nargout>4 || ~isempty(po)       % need to update robust statistics if >4 outputs or if plotting
0563         if  nr(19)==0               % if no outliers allowed in any case
0564             mr=mq;
0565             vr=vq;
0566         elseif nresr==nres          % if robust statistics are all up-to-date
0567             mr=mqr;
0568             vr=vqr;
0569         else                        % update robust statistics for all models
0570             for jq=1:length(nresq)             % loop for each model
0571                 if any(res(nresr+1:nres,1)==jq)   % if we have had any recent probe results for this model
0572                     guessj=pr(3,jq);                % guess rate (1/choices)
0573                     pscalej=qr(1,jq);            % prob range left after subtracting miss and guess probs
0574                     xtstdj=qr(3,jq);               % x position of target in std measure
0575                     xqj=xq(:,jq);               % SRT values of the pdf array
0576                     sqj=sq(:,jq);               % slope (or log slope) values in PDF
0577                     % first determine the SNR range to include
0578                     xej=mq(1,jq,1);  % mean SRT
0579                     sej=mq(2,jq,1);   % mean slope or log slope
0580                     resj=res(res(:,1)==jq,:);  % results for this model
0581                     if nr(11)                   % if log slope
0582                         sqisj=exp(sej)/qr(4,jq);    % inverse std dev of gaussian (proportional to slope) at mean slope
0583                     else
0584                         sqisj=sej/qr(4,jq);         % inverse std dev of gaussian (proportional to slope) at mean slope
0585                     end
0586                     % find upper bound for SNR with r=0
0587                     switch nr(10)   % switch on model type
0588                         case 1
0589                             xlim(2,jq)=xej-(xtstdj+log(nr(19)/(pscalej-nr(19))))/sqisj; % outlier if r=0 and x > xlim
0590                         case 2
0591                             xlim(2,jq)=xej+(norminv((pscalej-nr(19))/pscalej)-xtstdj)/sqisj; % outlier if r=1 and x < xlim
0592                     end
0593                     mou=resj(:,2)>=xlim(2,jq) & resj(:,3)==0;       % find high outliers
0594                     if nr(19)>guessj  % find lower bound for SNR with r=1
0595                         switch nr(10)   % switch on model type
0596                             case 1
0597                                 xlim(1,jq)=xej-(xtstdj+log(pscalej/(nr(19)-guessj)-1))/sqisj; % outlier if r=1 and x < xlim
0598                             case 2
0599                                 xlim(1,jq)=xej+(norminv((nr(19)-guessj)/pscalej)-xtstdj)/sqisj; % outlier if r=1 and x < xlim
0600                         end
0601                         mou=mou | resj(:,2)<=xlim(1,jq) & resj(:,3)==1;       % add in low outliers
0602                     end
0603                     if any(mou) % if any results must be excluded
0604                         % now recalculate pdf excluding bad points
0605                         if nr(11)                   % if log slope
0606                             sqisj=exp(sqj)/qr(4,jq);    % inverse std dev of gaussian (proportional to slope) for all grid points
0607                         else
0608                             sqisj=sqj/qr(4,jq);         % inverse std dev of gaussian (proportional to slope) for all grid points
0609                         end
0610                         wqj=repmat(pq0(2,jq)*(sqj-mq0(2,jq)).^2,1,nxq)+repmat(pq0(1,jq)*(xqj-mq0(1,jq))'.^2,nsq,1); % log probs of prior
0611                         switch nr(10)               % switch on model type
0612                             case 1                  % nr(10)=1: logistic model
0613                                 for j=find(mou==0)'
0614                                     r0=resj(j,3)==0;
0615                                     wqj=wqj+log(r0+(1-2*r0)*(guessj+pscalej*((1+exp(sqisj*xqj'-xtstdj-repmat(sqisj,1,nxq)*resj(j,2))).^(-1)))); %  P(l | r,x)
0616                                 end
0617                             case 2                  % nr(10)=2: cumulative gaussian model
0618                                 for j=find(mou==0)'
0619                                     r0=resj(j,3)==0;
0620                                     wqj=wqj+log(r0+(1-2*r0)*(guessj+pscalej*normcdf(repmat(sqisj,1,nxq)*resj(j,2)-sqisj*xqj'+xtstdj))); %  P(l | r,x)
0621                                 end
0622                         end
0623                         % now calculate robust means and modes
0624                         wqj=wqj(:);                         % turn 2D array into a column vector
0625                         ewqj=exp(wqj-max(wqj));             % unnormalized probability vector
0626                         ewqj=ewqj/sum(ewqj);
0627                         wqsxr=reshape(ewqj,nsq,nxq);         % normalized probabilities
0628                         pxr=sum(wqsxr,1);                             %  p(x0)
0629                         psr=sum(wqsxr,2);                             %  p(s0)
0630                         xer=pxr*xqj;                                  % E(x0)
0631                         ser=psr'*sqj;                                 % E(s0)
0632                         [pxpk,xmr]=max(pxr);                          % marginal mode in x
0633                         if xmr>1 && xmr<nxq                           % use quadratic interpolation in log prob if possible
0634                             [dum,xm2]=quadpeak(log(pxr(xmr-1:xmr+1))');
0635                             xmr=xmr+xm2-2;
0636                         end
0637                         xmr=(2-xmr)*xqj(1)+(xmr-1)*xqj(2);             % marginal mode(x0) (requires the xqi to be uniformly spaced)
0638                         [pspk,smr]=max(psr);
0639                         if smr>1 && smr<nsq                           % use quadratic interpolation in log prob if possible
0640                             [dum,sm2]=quadpeak(log(psr(smr-1:smr+1)));
0641                             smr=smr+sm2-2;
0642                         end
0643                         smr=(2-smr)*sqj(1)+(smr-1)*sqj(2);            % marginal mode(s0)
0644                         [wqpk,j]=max(ewqj);                         % find max of 1-D array
0645                         i=1+floor((j-1)/nsq);                       % convert to 2-D index (i,j)
0646                         j=j-nsq*(i-1);
0647                         if i>1 && i<nxq && j>1 && j<nsq             % use quadratic interpolation in log prob if possible
0648                             [dum,ji]=quadpeak(wqj(repmat((j-1:j+1)',1,3)+repmat(nsq*(i-2:i),3,1)));
0649                             i=i+ji(2)-2;
0650                             j=j+ji(1)-2;
0651                         end
0652                         xjr=(2-i)*xqj(1)+(i-1)*xqj(2);               % joint mode  x
0653                         sjr=(2-j)*sqj(1)+(j-1)*sqj(2);               % joint mode: s
0654                         xvr=pxr*(xqj.^2)-xer^2;                       % Var(x0)
0655                         svr=psr'*(sqj.^2)-ser^2;                      % Var(s0)
0656                         sxvr=ewqj'*(repmat(sqj-ser,nxq,1).*reshape(repmat(xqj'-xer,nsq,1),nsxq,1)); % Cov(s0*x0)
0657                         mqr(:,jq,:)=reshape([xer xjr xmr; ser sjr smr],[2 1 3]); % save means, joint modes and marginal modes
0658                         vqr(:,jq)=[xvr; sxvr; svr];                    % save covariance matrix
0659                     else   % if there are no outliers, just copy non-robust results
0660                         mqr(:,jq,:)= mq(:,jq,:);                    % use non-robust means
0661                         vqr(:,jq)= vq(:,jq);                          % and non-robust variances
0662                     end
0663                 end
0664             end
0665             nresr=nres;                 % all robust means/modes are now up-to-date
0666             mr=mqr;                     % send to output arguments
0667             vr=vqr;
0668         end
0669         if nr(11)    % if using log-slope
0670             mr(2,:,:)=exp(mr(2,:,:));  % convert to real slope
0671             vr(2,:)=vr(2,:).*mr(2,:,1); % correct the covariance
0672             vr(3,:)=vr(3,:).*mr(2,:,1).^2; % and the slope variance
0673         end
0674     end
0675     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0676     % now estimate the next probe SNR            %
0677     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0678     if ~numel(xz{iq})                               % if no list of probe SNRs was specified
0679         ytry=exp(0.25i*pi*(0:7))';                  % points around the circle
0680         ytry=[real(ytry) imag(ytry)];
0681         [vtry,dtry]=eig([xv sxv; sxv sv]);          % eigendecomposition of covariance matrix
0682         tryxs=repmat([xe,se],8,1)+nr(14)*ytry*sqrt(dtry)*vtry';
0683         pmin=0.05;                                  % target probe success probability
0684         if nr(11)
0685             tryxs(:,2)=qr(4,iq)*exp(-tryxs(:,2));     % convert log(slope) to std dev
0686         else
0687             tryxs(:,2)=qr(4,iq)*tryxs(:,2).^(-1);     % convert slope to std dev
0688         end
0689         dp=nr(15);  % maximum shift of probe value outside previous range
0690         switch nr(10)   % switch by model type
0691             case 1 % logistic
0692                 qmax=min(xmm(2,iq)+dp,max(tryxs(:,1)+(log((1-pmin)/pmin)-xtstd)*tryxs(:,2)));
0693                 qmin=max(xmm(1,iq)-dp,min(tryxs(:,1)+(log(pmin/(1-pmin))-xtstd)*tryxs(:,2)));
0694             case 2 % cumulative gaussian
0695                 qmax=min(xmm(2,iq)+dp, max(tryxs(:,1)+(norminv(1-pmin)-xtstd)*tryxs(:,2)));
0696                 qmin=max(xmm(1,iq)-dp,min(tryxs(:,1)+(norminv(pmin)-xtstd)*tryxs(:,2)));
0697         end
0698         dxt=max(nr(5),(qmax-qmin)/nxh);             % minimum step size of nr(5) [0.2 dB]
0699         xt=(qmin+qmax)/2+((1:nxh)-(1+nxh)/2)*dxt;
0700     else                                            % if a specific list of probe SNRs exists
0701         xzi=xz{iq};                                 % xzi is the list of available probe SNRs
0702         if numel(xzi)<=nxh                          % use all available probe SNRs if there are not too many
0703             xt=xzi;
0704         else
0705             [xt,ixt]=min(abs(xzi-xe));                                  % find the closest one to xe ** not necessarily optimum ***
0706             ixt=max(1,min(1+numel(xzi)-nxh,ixt-floor((1+nxh)/2)));      % arrange symmetrically around xt
0707             xt=xzi(ixt:min(ixt+nxh-1,numel(xzi)));
0708         end
0709     end
0710     nxhp=length(xt);  % xt are the potential probe SNRs
0711     % Now find the probe value that minimizes the expected value of the cost function
0712     % In the following: l = parameters of psychometric function, x = the probe SNR and r = the probe result
0713     switch nr(10)
0714         case 1
0715             prt=guess+pscale*((1+exp(repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)-repmat(sqis,nxq,1)*xt)).^(-1)); %  P(r=1 | l,x)
0716         case 2
0717             prt=guess+pscale*normcdf(repmat(sqis,nxq,1)*xt-repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)); %  P(r=1 | l,x)
0718     end
0719     wqt=repmat(ewqi,1,nxhp);
0720     hminr=zeros(2,1);                   % space for minimum expected cost function for each r0
0721     hminj=zeros(1,nxhp);   % space for minimum expected cost function for each x0
0722     if nr(18) % if doing look 2-ahead
0723         pl1=prt.*wqt;                       % posterior prob given success = p(l | x0,r0=1) unnormalized
0724         pl0=wqt-pl1;                        % posterior prob given failure = p(l | x0,r0=0) unnormalized
0725         pr0=sum(pl1,1);                     % p(r0=1 | x0)=Sum{P(r0=1 | l,x0)*P(l)} [note each column of wqt is normalized] (row vector)
0726         plxr=reshape([pl0./repmat(1-pr0,nsxq,1) pl1./repmat(pr0,nsxq,1)],nsxq,nxhp,2); % posterior prob p(l | x0,r0) column-normalized
0727         nx0=nxhp; % outer loop for each x0
0728         nr0=2; % inner loop for each r0
0729     else  % if only doing look 1-ahead
0730         nx0=1;  % only execute outer loop once
0731         nr0=1;  % only execute inner loop once
0732         pr0=0; % dummy value (used at end of outer loop)
0733     end
0734     hx2=zeros(nx0,nxhp,nr0);        % space for square array of expected cost functions (for possible plotting only)
0735     for j=1:nx0                     % loop for each possible probe SNR, x0, (or only once is look-1-ahead)
0736         for jr=1:nr0                 % loop for each possible test result, r0=jr-1 (or only once is look-1-ahead)
0737             if nr(18) % if doing look 2-ahead
0738                 wqt=repmat(plxr(:,j,jr),1,nxhp);         % posterior prob p(l | x0=xt(j),r0=jr-1) column-normalized
0739             end
0740             pl1=prt.*wqt;                     % posterior prob given success = p(l | x,r=1) unnormalized
0741             pl0=wqt-pl1;                      % posterior prob given failure = p(l | x,r=0) unnormalized
0742             prh=sum(pl1,1);                     % p(r | x)=Sum{P(r | l,x)*P(l)} [note each column of wqtjr is normalized] (row vector)
0743             pl1=pl1./repmat(prh,nsxq,1);        % posterior prob given success = p(l | x,r=1) normalized
0744             pl0=pl0./repmat(1-prh,nsxq,1);      % posterior prob given failure = p(l | x,r=0) normalized
0745             
0746             px1=squeeze(sum(reshape(pl1,nsq,nxq,[]),1));    % p(x0 | x,r=1)
0747             px0=squeeze(sum(reshape(pl0,nsq,nxq,[]),1));    % p(x0 | x,r=0)
0748             ps1=squeeze(sum(reshape(pl1,nsq,nxq,[]),2));    % p(s0 | x,r=1)
0749             ps0=squeeze(sum(reshape(pl0,nsq,nxq,[]),2));    % p(s0 | x,r=0)
0750             xet1=xqi'*px1;                                  % E(x0 | x,r=1)
0751             xvt1=(xqi.^2)'*px1-xet1.^2;                     % Var(x0 | x,r=1)
0752             xet0=xqi'*px0;                                  % E(x0 | x,r=0)
0753             xvt0=(xqi.^2)'*px0-xet0.^2;                     % Var(x0 | x,r=0)
0754             xvt=xvt1.*prh+xvt0.*(1-prh);                    % E(Var(x0 | x ))
0755             set1=sqi'*ps1;                                  % E(s0 | x,r=1)
0756             svt1=(sqi.^2)'*ps1-set1.^2;                     % Var(s0 | x,r=1)
0757             set0=sqi'*ps0;                                  % E(s0 | x,r=0)
0758             svt0=(sqi.^2)'*ps0-set0.^2;                     % Var(s0 | x,r=0)
0759             svt=svt1.*prh+svt0.*(1-prh);                    % E(Var(s0 | x ))
0760             xht1=sum(log(px1).*px1,1);                      % -H(x0 | x,r=1)
0761             xht0=sum(log(px0).*px0,1);                      % -H(x0 | x,r=0)
0762             xht=(xht1.*prh+xht0.*(1-prh))*(xqi(1)-xqi(2));    % h(x0 | x)
0763             sht1=sum(log(ps1).*ps1,1);                      % -H(s0 | x,r=1)
0764             sht0=sum(log(ps0).*ps0,1);                      % -H(s0 | x,r=0)
0765             sht=(sht1.*prh+sht0.*(1-prh))*(sqi(1)-sqi(2));    % h(s0 | x)
0766             switch nr(9) % cost function: 1=variance, 2=entropy
0767                 case 1
0768                     hx=(xvt + nr(4)*svt)/(1+nr(4));               % expected cost function for each possible test SNR
0769                 case 2
0770                     hx=(xht + nr(4)*sht)/(1+nr(4));                % expected cost function for each possible test SNR                                          % find the minimum of cost function
0771             end
0772             hx2(j,:,jr)=hx; % save expected cost function for possible plotting
0773             [hminr(jr),ix]=min(hx);       % find the minimum of cost function for each value of r (0 or 1)
0774             %         fprintf('Probe range: %.3g %.3g; choose %.3g\n',xt(1),xt(end),xt(ix));
0775         end
0776         hminj(j)=[1-pr0(j) pr0(j)]*hminr; % expected cost function for each x0 assuming x1 is chosen optimally
0777     end
0778     if nr(18) % if doing look 2-ahead
0779         [hminr(1),ix]=min(hminj);                         % find the minimum of cost function
0780     end
0781     % now ix indexes the optimal probe choice and hminr(1) is the resultant expected cost function
0782     xn(iq)=xt(ix);  % save the next probe snr for this model
0783     % calculate the expected decrease in cost function (to decide which model to probe next)
0784     switch nr(9)  % cost function: 1=variance, 2=entropy
0785         case 1                                              % minimize variance
0786             hn(iq)=(xv + nr(4)*sv)/(1+nr(4))-hminr(1);        % save the expected decrease in cost function for this model
0787         case 2                                              % minimize entropy
0788             hn(iq)=(xh + nr(4)*sh)/(1+nr(4))-hminr(1);        % save the expected decrease in cost function for this model
0789     end
0790     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0791     %  now do plotting                                 %
0792     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0793     if ~isempty(po)
0794         h=gcf;
0795         if isreal(h)
0796             fig0=round(h); % get figure number
0797         else
0798             fig0=get(h,'number');  % in new versions of matlab use this method
0799         end
0800         axm=[3 -2; 2 -1]; % matrix to extend an axis beackwards by two steps
0801         for i=1:length(po)
0802             poi=po(i);
0803             switch poi
0804                 case 'p'    % plot pdf
0805                     nxq=nr(1);
0806                     nsq=nr(2);
0807                     xqj=[axm*xqi(1:2);xqi];
0808                     sqj=[axm*sqi(1:2);sqi];
0809                     figure(fig0);
0810                     imagesc(xqj,sqj,[zeros(1,2) px/pxpk; zeros(1,nxq+2);ps/pspk zeros(nsq,1) wqsx/wqpk]);
0811                     hold on
0812                     plot(xe,se,'ko',xm,sm,'k^',xj,sj,'k*');
0813                     plot(xqj(2),se,'wo',xqj(2),sm,'w^',xqj(2),sj,'w*');
0814                     plot(xe,sqj(2),'wo',xm,sqj(2),'w^',xj,sqj(2),'w*');
0815                     if any(mqr(:,iq,:)~=mq(:,iq,:))
0816                         plot(mqr(1,iq,1),mqr(2,iq,1),'go',mqr(1,iq,3),mqr(2,iq,3),'g^',xj,sj,'g*');
0817                         plot(xqj(end),mqr(2,iq,1),'go',xqj(end),mqr(2,iq,3),'g^',xqj(end),mqr(2,iq,2),'g*');
0818                         plot(mqr(1,iq,1),sqj(end),'go',mqr(1,iq,3),sqj(end),'g^',mqr(1,iq,2),sqj(end),'g*');
0819                         title('Joint pdf: o mean, * mode, \Delta marg mode (green=robust)');
0820                     else
0821                         title('Joint pdf: o mean, * mode, \Delta marginal mode');
0822                     end
0823                     t=linspace(0,2*pi,200);
0824                     xcir=cos(t);
0825                     scir=sin(t);
0826                     vcir=sqrt((sv*xcir.^2+xv*scir.^2-2*sxv*xcir.*scir)/(xv*sv-sxv^2));
0827                     plot(xe+xcir./vcir,se+scir./vcir,'w-');
0828                     hold off
0829                     axis 'xy';
0830                     %     colorbar;
0831                     %     cblabel('Relative Probability');
0832                     xlabel(sprintf('SNR @ %d%%SRT (dB)',round(pr(1)*100)));
0833                     if nr(11)
0834                         ylabel('Log psychometric Slope at threshold (prob/dB)');
0835                     else
0836                         ylabel('Psychometric Slope at threshold (prob/dB)');
0837                     end
0838                     fig0=fig0+1;
0839                 case 'c' % plot cost evolution
0840                     resi=res(res(:,1)==iq,:);
0841                     nresi=size(resi,1);
0842                     if nresi>0
0843                         figure(fig0);
0844                         plot(1:nresi,resi(:,6)/(1+nr(4)),'--r',1:nresi,(resi(:,6) + nr(4)*resi(:,6))/(1+nr(4)),'-b');
0845                         axisenlarge([-1 -1.03]);
0846                         xlabel('Trial Number');
0847                         switch nr(9)  % cost function: 1=variance, 2=entropy
0848                             case 1                                              % minimize variance
0849                                 ylabel('Weighted Variance: SRT and Total');
0850                             case 2                                              % minimize entropy
0851                                 ylabel('Weighted Entropy: SRT and Total');
0852                         end
0853                         title ('Cost function evolution');
0854                         fig0=fig0+1;
0855                     end
0856                 case 'h' % plot history
0857                     resi=res(res(:,1)==iq,:);
0858                     nresi=size(resi,1);
0859                     if nresi>0
0860                         if nr(11)
0861                             resi(:,5)=exp(resi(:,5)); % convert to real slope
0862                             ylim=xe+([-2 0 thresh 1 3]-thresh)/exp(se);
0863                         else
0864                             ylim=xe+([-2 0 thresh 1 3]-thresh)/se;
0865                         end
0866                         figure(fig0);
0867                         plot(1:nresi,resi(:,[4 4 4])+resi(:,5).^(-1)*([0 thresh 1]-thresh),'-b');
0868                         hold on;
0869                         % could plot outliers in a different colour by checking mou(res(:,1)==iq)
0870                         msky=resi(:,3)==1;
0871                         plot(find(msky),resi(msky,2),'r+');
0872                         msky=resi(:,3)==0;
0873                         plot(find(msky),resi(msky,2),'bo');
0874                         gcaylim=get(gca,'ylim');
0875                         poutb=0;
0876                         if nr(19)>0 && xlim(2,iq)<gcaylim(2) % plot upper outlier bound
0877                             plot([1;nresi],[1;1]*xlim(2,iq),'-r');
0878                             texthvc(nresi,xlim(2,iq),'max o','rtr');
0879                         end
0880                         if nr(19)>guess && xlim(1,iq)>gcaylim(1) % plot upper outlier bound
0881                             plot([1;nresi],[1;1]*xlim(1,iq),'-r');
0882                             texthvc(nresi,xlim(1,iq),'min+','rbr');
0883                         end
0884                         hold off;
0885                         xlabel('Trial Number');
0886                         ylabel('Probe SNR (dB) [o,x=result]');
0887                         axisenlarge(-1.02);
0888                         title('SRT mean + slope intercepts');
0889                         fig0=fig0+1;
0890                     end
0891                 case 'x' % plot probe choice
0892                     figure(fig0);
0893                     if nr(18) % if doing look 2-ahead
0894                         xtaxm=xt(1:2)*axm';
0895                         hx2e=hx2(:,:,1).*repmat(1-pr0',1,nxhp)+hx2(:,:,2).*repmat(pr0',1,nxhp);
0896                         [dum,ix2m]=min(hx2,[],2); % ix2 indices to minimize hx2
0897                         % ==========  note that the following line assumes that the xt()
0898                         % are uniformly spaced (not necessarily true)
0899                         imagesc([xtaxm xt],xt,[hminj' repmat(max(hx2e(:)),nxhp,1) hx2e]);
0900                         axis xy;
0901                         colorbar;
0902                         hold on
0903                         [hnmin,ii]=max(hn);         % find the next model to probe
0904                         plot(xtaxm(1),xn(iq),'>w');
0905                         if iq==ii
0906                             plot(xt(1:2)*axm(1,:)',xn(iq),'*w'); % circle the next model to be probed
0907                         end
0908                         plot([xtaxm(1) xt(end)],xn(iq)*[1 1],':w');
0909                         plot(xt(ix2m(:,1,1)),xt,'ow');
0910                         plot(xt(ix2m(:,1,2)),xt,'+w');
0911                         hold off
0912                         xlabel('SNR probe 2 (dB) [+,o=best for R1=1,0]');
0913                         ylabel(sprintf('SNR probe 1 (dB) [> = %.2f]',xn(iq)));
0914                         title('Expected Cost after 2 probes');
0915                     else % if doing look 1-ahead only
0916                         plot(xt,hx2(1,:,1),'b');
0917                         axisenlarge([-1 -1.03]);
0918                         caylim=get(gca,'ylim');
0919                         hold on
0920                         plot(xn([iq iq]),caylim,':k');
0921                         xlabel(sprintf('Next probe SNR (dB) [min @ %.2f]',xn(iq)));
0922                         ylabel('Cost')
0923                         title('Expected Cost after next probe');
0924                     end
0925                     fig0=fig0+1;
0926             end
0927         end
0928     end
0929 end
0930 if iq~=0
0931     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0932     %  now select the appropriate model to probe next  %
0933     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0934     [hnmin,ii]=max(hn);         % chose model with the biggest expected decrease
0935     xx=xn(ii);
0936     m=mq;
0937     v=vq;
0938     if nr(11)    % if using log-slope
0939         m(2,:,:)=exp(m(2,:,:));  % convert to real slope
0940         v(2,:)=v(2,:).*m(2,:,1); % correct the covariance
0941         v(3,:)=v(3,:).*m(2,:,1).^2; % and the slope variance
0942     end
0943 else
0944     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0945     %     iq=0 means output  model parameters and probe results   %
0946     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0947     xx=pr;
0948     ii=nr;
0949     m=res(1:nres,:);
0950     if ~nargout                     % print them if no output arguments to call
0951         fid=LOG+(LOG==0);           % use standard out if LOG==0
0952         pdesc={'Threshold','Miss prob','Guess Prob','Min SNR','Max SNR','Min Slope','Max Slope'};
0953         fprintf(LOG,'\n*********\nModel-specific Parameters\n');
0954         for i=1:7
0955             fprintf(LOG,'%4d) %s: ',i,pdesc{i});
0956             if size(pr,2)>1
0957                 fprintf(LOG,'%.5g, ',pr(i,1:size(pr,2)-1));
0958             end
0959             fprintf(LOG,'%.5g\n',pr(i,size(pr,2)));
0960         end
0961         qdesc={'nx  SNR values in PDF', ...
0962             'ns  Slope values in PDF', ...
0963             'nh  Probe SNR values to evaluate', ...
0964             'cs  Weighting of slope relative to SRT in cost function', ...
0965             'dh  Min step size in dB for probe SNRs', ...
0966             'sl  Min slope at threshold', ...
0967             'kp  Std deviations of the pdfs to keep', ...
0968             'hg  Amount to grow expected gains in ni trials', ...
0969             'cf  Cost function', ...
0970             'pm  Psychometric model', ...
0971             'lg  Use log slope as parameter', ...
0972             'pm  Prior std devs in semi-width', ...
0973             'pf  Integrated probability floor', ...
0974             'ts  Number of std devs to explore', ...
0975             'dp  Max shift in probe value', ...
0976             'it  Grid interpolation threshold', ...
0977             'at  Grid boundary tolerance', ...
0978             'la  Look 2-ahead when choosing probe', ...
0979             'op  Oulier probability threshold', ...
0980             'rx  min grid shrink factor'};
0981         qoptf=[9 10]; % fields with options
0982         qoptval={'variance','entropy'; ...
0983             'logistic','cumulative gaussian'};
0984         fprintf(LOG,'\nShared Parameters\n');
0985         for i=1:length(nr)
0986             fprintf(LOG,'%4d) %s: ',i,qdesc{i});
0987             j=find(qoptf==i,1);
0988             if numel(j)
0989                 fprintf(LOG,'%d=%s\n',nr(i),qoptval{j,nr(i)});
0990             else
0991                 fprintf(LOG,'%.5g\n',nr(i));
0992             end
0993         end
0994         fprintf(LOG,'\n');
0995     end
0996 end

Generated on Tue 19-Sep-2017 12:07:31 by m2html © 2003