zest.m

function [thresh, sd_pdf] = zest(response, init)
%
% Implements ZEST threshold procedures with special considerations for BAT
% test. Parts of code originally written by Peter Marvit (author of JASA paper
% on ZEST cited below). Optimized and elaborated by Brian Hurley.
%
% This routine generalizes the ZEST procedure. The first invocation should
% have a second parameter "init" (actually a structure) that specifies the
% various parameters to initialize the ZEST routine. In particular, init
% should have the following fields:
%   zestA           PDF scale factor
%   zestB           falling PDF slope
%   zestC           rising PDF slope
%   zestmaxrange    the maximum value for possible estimates (in log units)
%   zestminrange    the minimum value for possible estimates (in log units)
%   zestfa          False alarm rate estimate (0<= fa <=1)
%   zestmiss        Miss rate estimate (0<= miss <=1)
%   zestbeta        Response function slope
%   zesteta         "sweat factor" or response criterion parameter
%
% If the init structure is passed, then the value of "response" (though
% necessary) is ignored.
%
% Once initialized, just response (=0 or 1) is passed to ZEST. Based on the
% current set of parameters, the routine returns the most probable mean
% vlaue of the probability density function ('thresh') that can either be 
% used for the next trial OR as a final threshold estimate. Common criteria 
% used for terminating ZEST are (1) some fixed number of trials reached, and 
% (2) standard deviation of PDF (returned here as 'sd_pdf') falls to a 
% critical value. Criterion 2 is best for minimizing number of trials.
%
% The function returns 'thresh,' which is the trial-to-trial estimate of
% threshold.
%
% For details about the parameters, see Marvit, et al. (2003), JASA
% 113(6):3348-3361.
%
% ----------------------
% written May 2014 by Brian K. Hurley. Revised Sep 2016 for use with BAT
% test.

%
% Default values provided for parameters not provided
persistent T numbersteps A B C maxrange minrange fa miss beta eta q meanpdf convert

%% INITIALiZE THRESHOLD PARAMS
if nargin > 1
  
  % Parameters for Initial PDF for ZEST
  if isfield(init, 'zestA')
    A = init.zestA;
  else
    A = 1;    % Scale factor to start with unity PDF area; A in Marvit et al.
  end
  if isfield(init, 'zestB')
    B = init.zestB;
  else
    B = 2.5;  % Falling slope of PDF; B in Marvit et al.
  end
  if isfield(init, 'zestC')
    C = init.zestC;
  else
    C = 2.5;  % Rising slope of PDF; C in Marvit et al.
  end
  if isfield(init, 'zestmaxrange')
    maxrange = init.zestmaxrange;
  else
    maxrange = 2.5; % log max value of possible estimates
  end
  if isfield(init, 'zestminrange')
    minrange = init.zestminrange; % log min value of possible estimates (pre-log min value must be > 0; e.g., .01 rather than 0)
  else
    minrange = -2.5;
  end
  
  % Parameters for Response function (Weibull function) (P(yes) given stimulus)
  if isfield(init, 'zestfa')
    fa = init.zestfa;
  else
    fa = 0.10;   % false-alarm rate; gamma in Marvit et. al
  end
  if isfield(init, 'zestmiss')
    miss = init.zestmiss;
  else
    miss = 0.02; % false-negative (miss) rate; delta in Marvit et al.
  end
  if isfield(init, 'zestbeta')
    beta = init.zestbeta;
  else
    beta = 6;    % slope of weibull response function; beta in Marvit et al.
  end
  if isfield(init, 'zesteta')
    eta = init.zesteta;
  else
    eta = 0;     % "sweat factor" or response criterion; eta in Marvit et al.
  end
  
  % Starting params
  if isfield(init, 'zest_init_dl') % initial difference limen
    delta_L = init.zest_init_dl;
  else
    delta_L = 15; % initial difference limen and starting deviant magnitude; delta L in Marvit et al.
  end
  
  % convert to log
  init_t = log10(delta_L);
  
  % Finally, flag the output values that should be converted back to
  % original units.
  % Possibilities include 'delta_L' (level change threshold),'sd_pdf' (SD
  % of threshold PDF)
  if isfield(init, 'zestconvert')
    convert = init.zestconvert;
  else
    convert = {'delta_L', 'sd_pdf'};
  end
  
  % Create a discrete vector of the extent/range of possible DLs
  num_steps = round((10^maxrange)/.01); %
  % Linear DL's can vary from minrange to maxrange in .01 unit increments
  T = linspace(minrange, maxrange, num_steps);
  
  numbersteps = length(T);
  
  % Calculate the initial PDF
  q = A ./ (B*exp(-C*(T-init_t)) + C*exp(B*(T-init_t)));
  meanpdf = init_t;
  
  % De-log initial meanpdf if desired
  if any(strcmp(convert,'delta_L'))
    thresh=10^meanpdf;
  end
  
else % if nargin == 1
  % If we just have a response, calculate the next thresh. The prior
  % thresh estimate (stimulus difference level) was meanpdf.
  
  % Psychometric function (Weibull)
  % p is  model prob of resp given log_lev_diff if true log_DL is T
  p = 1 - miss - ((1-fa-miss) * exp(-10.^(beta * (meanpdf-T+eta))));
  % if response was No
  if response==0
    p=1-p;
  end
  
  % Compute the next q (next PDF)
  q = p .* q;
  
  % Calculate new PDF mean (threshold)
  meanpdf = sum(T.*q) / sum(q);
  
  % Calculate variance and standard deviation of pdf
  log_v = sum(((T-meanpdf).^2).*q) / sum(q);
  log_s = sqrt(log_v);
  
  % Final conversions, if any. Transforms speficied vars from log back to
  % original units.
  if any(strcmp(convert,'delta_L'))
    % de-logify
    thresh=10^meanpdf;
  else
    thresh=meanpdf;
  end
  if any(strcmp(convert, 'sd_pdf'))
    sd_pdf = 10^log_s;
  else
    sd_pdf = log_s;
  end
  
end %if nargin>1

end