# v_psycest

## PURPOSE

V_PSYCEST estimate multiple psychometric functions

## SYNOPSIS

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

## DESCRIPTION

    V_PSYCEST estimate multiple psychometric functions

     Usage: [xx,ii,m,v]=v_psycest(-n,p,q,xp,lf)    % initialize n models
            [xx,ii,m,v]=v_psycest(i,x,r)           % supply a trial result to v_psycest
            [xx,ii,m,v]=v_psycest(i,x,r,o)         % supply a trial result to v_psycest and plot
                        v_psycest(i,o)             % plot pdf of model i with plot options o
            [xx,ii,m,v,mr,vr]=v_psycest(i,o)       % Also determine robust outputs which exclude outliers [takes longer]
                  [p,q,msr]=psycest(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=v_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 (exclusing outliers)
                       m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
              vr(3,n)  robust estimated covariance matrix entries  (exclusing outliers):
                       [var(srt) cov(srt,slope) var(slope)]'
              msr(:,7)  List probe snrs and results: [model probe-snr result xe se xv sv]

     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 v_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, v_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 v_entropy of the estimated SNR and slope distributions.
              My experience is that v_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:

- v_axisenlarge V_AXISENLARGE - enlarge the axes of a figure (f,h)
- v_quadpeak V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array
- v_texthvc V_TEXTHVC - write text on graph with specified alignment and colour

This function is called by:

- v_psycdigit V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli

## SOURCE CODE

    0001 function [xx,ii,m,v,mr,vr]=v_psycest(iq,x,r,xp,lf)
    0002 %V_PSYCEST estimate multiple psychometric functions
    0003 %
    0004 % Usage: [xx,ii,m,v]=v_psycest(-n,p,q,xp,lf)    % initialize n models
    0005 %        [xx,ii,m,v]=v_psycest(i,x,r)           % supply a trial result to v_psycest
    0006 %        [xx,ii,m,v]=v_psycest(i,x,r,o)         % supply a trial result to v_psycest and plot
    0007 %                    v_psycest(i,o)             % plot pdf of model i with plot options o
    0008 %        [xx,ii,m,v,mr,vr]=v_psycest(i,o)       % Also determine robust outputs which exclude outliers [takes longer]
    0009 %              [p,q,msr]=psycest(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=v_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 (exclusing outliers)
    0061 %                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
    0062 %          vr(3,n)  robust estimated covariance matrix entries  (exclusing outliers):
    0063 %                   [var(srt) cov(srt,slope) var(slope)]'
    0064 %          msr(:,7)  List probe snrs and results: [model probe-snr result xe se xv sv]
    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 v_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, v_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 v_entropy of the estimated SNR and slope distributions.
    0108 %          My experience is that v_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: v_psycest.m 10865 2018-09-21 17:22:45Z 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 Lesser General Public License as published by
    0137 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0144 %
    0145 %   You can obtain a copy of the GNU Lesser General Public License from
    0146 %   https://www.gnu.org/licenses/ .
    0147 %
    0148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0149
    0150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0151 % Bugs/Suggestions:
    0152 % (2)use a structure for input parameters
    0153 % (3) v_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 v_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 v_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]=v_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 v_entropy h(x0)
    0550     sh=(ps'*log(ps))*(sqi(1)-sqi(2));              % differential v_entropy h(s0)
    0551     if nargin==3
    0552         switch nr(9) % cost function: 1=variance, 2=v_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]=v_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=v_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=v_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 v_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                     %     v_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                         v_axisenlarge([-1 -1.03]);
    0846                         xlabel('Trial Number');
    0847                         switch nr(9)  % cost function: 1=variance, 2=v_entropy
    0848                             case 1                                              % minimize variance
    0849                                 ylabel('Weighted Variance: SRT and Total');
    0850                             case 2                                              % minimize v_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                             v_texthvc(nresi,xlim(2,iq),'max o','rtr');
    0879                         end
    0880                         if nr(19)>guess && xlim(1,iq)>gcaylim(1) % plot lower outlier bound
    0881                             plot([1;nresi],[1;1]*xlim(1,iq),'-r');
    0882                             v_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                         v_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                         v_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