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))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