Home > psychophysics > zest > zest.m

zest

PURPOSE ^

thresh=zest(response,init)

SYNOPSIS ^

function [thresh,pdfinfo]=zest(response,params)

DESCRIPTION ^

 thresh=zest(response,init)

 This routine generalizes the ZEST procedure. The first invocation should
 have a second parameter "params" (actually a structure) that specifies the
 various parameters to initialize the ZEST routine. In particular, params
 should have the following fields:
   zestA   PDF scale factor
   zestB   falling PDF slope
   zestC   rising PDF slope
   zestmaxrange    the maximum value for possible estimates (usually in
                   log units)
   zestminrange    the minimum value for possible estimates (usually in
                   log units)
   zestfa  False alarm rate estimate (0<= fa <=1)
   zestmiss    Miss rate estimate (0<= miss <=)
   zestbeta    Response function slope
   zesteta     "sweat factor" or response criterion parameter
 If there is a passed structure params, then the value of "response" (though
 necessary) is ignored.

 Otherwise, just response (=0 or 1) is passed to ZEST. Based on teh
 current set of parameters, the routine returns the most probable mean
 vlaue of the PDF that can either be used for the next trial OR as a final
 threshold estimate.

 To be done: return the variance of the current PDF so one could stop at a
 particular confidence interval rather than after so many trials.

 Note that this code is currently written with default values for a
 modulation detection task. 
 
 For details about the parameters, see Marvit, et al. (2003), JASA
 113(6):3348-3361.

 In the main code, we specify modulation depth as percentage, but it is
 better expressed and manipulated in a dB scale. So: Convert modulation
 depth, m, to delta L in dB 
   DLdB=20*log10((1+m)./(1-m))

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [thresh,pdfinfo]=zest(response,params)
0002 % thresh=zest(response,init)
0003 %
0004 % This routine generalizes the ZEST procedure. The first invocation should
0005 % have a second parameter "params" (actually a structure) that specifies the
0006 % various parameters to initialize the ZEST routine. In particular, params
0007 % should have the following fields:
0008 %   zestA   PDF scale factor
0009 %   zestB   falling PDF slope
0010 %   zestC   rising PDF slope
0011 %   zestmaxrange    the maximum value for possible estimates (usually in
0012 %                   log units)
0013 %   zestminrange    the minimum value for possible estimates (usually in
0014 %                   log units)
0015 %   zestfa  False alarm rate estimate (0<= fa <=1)
0016 %   zestmiss    Miss rate estimate (0<= miss <=)
0017 %   zestbeta    Response function slope
0018 %   zesteta     "sweat factor" or response criterion parameter
0019 % If there is a passed structure params, then the value of "response" (though
0020 % necessary) is ignored.
0021 %
0022 % Otherwise, just response (=0 or 1) is passed to ZEST. Based on teh
0023 % current set of parameters, the routine returns the most probable mean
0024 % vlaue of the PDF that can either be used for the next trial OR as a final
0025 % threshold estimate.
0026 %
0027 % To be done: return the variance of the current PDF so one could stop at a
0028 % particular confidence interval rather than after so many trials.
0029 %
0030 % Note that this code is currently written with default values for a
0031 % modulation detection task.
0032 %
0033 % For details about the parameters, see Marvit, et al. (2003), JASA
0034 % 113(6):3348-3361.
0035 %
0036 % In the main code, we specify modulation depth as percentage, but it is
0037 % better expressed and manipulated in a dB scale. So: Convert modulation
0038 % depth, m, to delta L in dB
0039 %   DLdB=20*log10((1+m)./(1-m))
0040 
0041 % 01/27/06 Petr Janata
0042 % - isolated the initialization routine a bit more
0043 % - persistent variable that did not need to be persistent were removed
0044 % - changed beta to b to avoid overriding beta()
0045 
0046 % CHeck Default values. Currently for 2AFC
0047 persistent is_inited
0048 persistent T fa miss b eta q meanpdf convert use_log
0049 
0050 % Initialize
0051 if isempty(is_inited)
0052   fprintf('Initializing PDF ...\n');
0053   
0054   % Parameters for Initial PDF for ZEST
0055   if isfield(params,'zestA') A=params.zestA;, else, 
0056     A=1;    % Scale factor to start with unity pdf area
0057   end
0058   if isfield(params,'zestB') B=params.zestB;, else, 
0059     B=3.0;  % Falling slope
0060   end
0061   if isfield(params,'zestC') C=params.zestC;, else, 
0062     C=1.5;  % Rising slope
0063   end    
0064   if isfield(params,'zestmaxrange') maxrange=params.zestmaxrange;, else, 
0065     maxrange=-2.5;
0066   end
0067   if isfield(params,'zestminrange') minrange=params.zestminrange;, else, 
0068     minrange=2.5;
0069   end
0070   
0071   % Parameters for Response function (Weibull function) (P(yes) given stimulus)
0072   if isfield(params,'zestfa') fa=params.zestfa;, else, 
0073     fa=0.50;   %gamma in the text, false alarm rate (guess rate for 2AFC)
0074   end
0075   if isfield(params,'zestmiss') miss=params.zestmiss;, else, 
0076     miss=0.01; %delta in the text, miss rate (1/2 inattention rate for 2AFC)
0077   end
0078   if isfield(params,'zestbeta') b=params.zestbeta;, else, 
0079     b=.6;    %beta in the text, slope of response function
0080   end
0081   if isfield(params,'zesteta') eta=params.zesteta;, else, 
0082     eta=0;     %eta in the text, "sweat factor" or response criterion parameter
0083   end
0084   
0085   % If a vector of DLs has been provided, use it
0086   if isfield(params,'T')
0087     T = params.T;
0088   else
0089     % Create a discrete vector of the extent/range of possible DLs
0090     T=linspace(minrange,maxrange,1000); % Linear DL's can vary from .01 dB to 20 dB in .01 dB increments
0091   end
0092   
0093   % Starting params
0094   if isfield(params,'zestinit') 
0095     m=params.zestinit; 
0096   else 
0097     m=.5; % Initially, start with 50% AMdepth
0098   end
0099   
0100   % Are we doing everything on a logscale
0101   if isfield(params,'logscale')
0102     use_log = params.logscale;
0103   else
0104     use_log = 1; 
0105   end
0106   
0107   % But convert it to log(dB)
0108   if use_log
0109     % Use initial "guess" as midpoint of PDF
0110     init_t=log10(20*log10((1+m)/(1-m))); 
0111   else                     
0112     init_t = m;
0113   end
0114   
0115   % Finally, flag to know how to return the threshold
0116   % POssibilities include 'm' (AMdepth->(0,100)%), 'mdB', 'logmdB'
0117   
0118   % 01/27/06 PJ changed (incorrect) reference to params.convert to params.zestconvert
0119   if isfield(params,'zestconvert') convert=params.zestconvert; else,   
0120     convert='m';
0121   end
0122   
0123   % Calculate the initial PDF
0124   q=A./(B*exp(-C*(T-init_t)) + C*exp(B*(T-init_t)));
0125   
0126   % Normalize area under q to be 1
0127   q = q/sum(q);
0128   
0129   meanpdf=init_t;
0130 
0131   thresh = convert_thresh(meanpdf, convert, use_log);
0132 
0133   pdfinfo.q = q;
0134   pdfinfo.T = T;
0135   pdfinfo.meanpdf = meanpdf;
0136   
0137   is_inited = 1;
0138   return
0139 end % if isempty(is_inited)
0140 
0141 % If we just have a response, calculate the next thresh. The prior
0142 % thresh estimate (stimulus value) was meanpdf.
0143 
0144 %   Psychometric function (Weibull):p is  model prob of resp given log_lev_diff if true log_DL is T
0145 p=1-miss-((1-fa-miss)*exp(-10.^(b*(meanpdf-T+eta))));
0146 if response==0
0147   p=1-p;
0148 end
0149 
0150 % Compute the next q (the next pdf)
0151 q=p.*q;
0152 
0153 % Normalize q
0154 q = q/sum(q);
0155 
0156 meanpdf=sum(T.*q)/sum(q); % Calculate new midpoint
0157 
0158 % Make sure that the new midpoint falls on a location on our scale
0159 if isempty(find(T == meanpdf))
0160   tidx = max(find(T < meanpdf));
0161   if isempty(tidx)
0162     tidx = 1;
0163   end
0164   meanpdf = T(tidx);
0165 end
0166 
0167 thresh = convert_thresh(meanpdf, convert, use_log);
0168 
0169 pdfinfo.q = q;
0170 pdfinfo.T = T;
0171 pdfinfo.meanpdf = meanpdf;
0172 
0173 return
0174   
0175 function thresh=convert_thresh(meanpdf, convert_type, use_log)
0176   switch convert_type
0177     case 'm'
0178       if use_log
0179     % De-logify
0180     antimean=10^meanpdf;
0181     % Now, de-dBify
0182     thresh=((10^(antimean/20))-1)/((10^(antimean/20))+1);
0183       else
0184     thresh = meanpdf;
0185       end
0186     case 'mdB'
0187       if use_log
0188     % Just de-logify
0189     thresh=10^meanpdf;
0190       end
0191     otherwise
0192       thresh=meanpdf;
0193   end
0194

Generated on Sun 24-Mar-2019 04:00:51 by m2html © 2003