ZEST_marvit.m

function thresh=ZEST_marvit(response,init)
% thresh=ZEST(response,init)
% 
% Replicates ZEST analyes in Marvit et al. (2003), JASA. Based largely on
% code received by Petr from Marvit.
%
% 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 (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 init, 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.
% 
% For details about the parameters, see Marvit, et al. (2003), JASA
% 113(6):3348-3361.
%
% 15 May 2014 BH 

% CHeck Default values. Currently for yes/no
persistent T numbersteps A B C maxrange minrange fa miss beta eta q meanpdf convert

PLOT_pdf = 0;
fid = fopen('T.txt','w');

% Initialize
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
	end
	if isfield(init,'zestB') B=init.zestB;, else, 
		B=3.0;  % Falling slope
	end
	if isfield(init,'zestC') C=init.zestC;, else, 
		C=1.5;  % Rising slope
	end	
	if isfield(init,'zestmaxrange') maxrange=init.zestmaxrange;, else, 
		maxrange=-2.5;
	end
	if isfield(init,'zestminrange') minrange=init.zestminrange;, 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.50;   %gamma in the text, false alarm rate (guess rate for 2AFC)
	end
	if isfield(init,'zestmiss') miss=init.zestmiss;, else, 
		miss=0.01; %delta in the text, miss rate (1/2 inattention rate for 2AFC)
	end
	if isfield(init,'zestbeta') beta=init.zestbeta;, else, 
		beta=.6;    %beta in the text, slope of response function
	end
	if isfield(init,'zesteta') eta=init.zesteta;, else, 
		eta=0;     %eta in the text, "sweat factor" or response criterion parameter
	end
	
	% Starting params
    if isfield(init,'zestinit_diffLvl')
        delta_L = init.zestinit_diffLvl;
    else
        delta_L = 3; % start w/ 3 dB
    end
    
    % CONVERT BACK TO LOG
    init_t = log10(delta_L); %delta_L; 

%     % But convert it to log(dB)
%     init_t=log10(20*log10((1+m)/(1-m))); % Use initial "guess" as midpoint of PDF
% 	
	% Finally, flag to know how to return the threshold
	% POssibilities include 'delta_L' (level change in dB), 'log_delta_L' 
	if isfield(init,'zestconvert') convert=init.zestconvert, else,
		convert='delta_L';
	end
		
	% Create a discrete vector of the extent/range of possible DLs
	T=linspace(minrange,maxrange,2000); % Linear DL's can vary from .01 dB to 20 dB in .01 dB increments
	
    %fprintf(fid,'%f ',T);
    %fclose(fid);
    numbersteps=length(T);
	
	% Calculate the initial PDF
	q=A./(B*exp(-C*(T-init_t)) + C*exp(B*(T-init_t)));
	meanpdf=init_t;
else % if nargin>1
	% If we just have a response, calculate the next thresh. The prior
	% thresh estimate (stimulus value) 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==0
		p=1-p;
	end

	% Compute the next q (the next pdf)
	q=p.*q;
	meanpdf=sum(T.*q)/sum(q); % Calculate new midpoint
    v = sum(((T-meanpdf).^2).*q)/sum(q);
    s = sqrt(v);
    s_dB = 10^s;
	%thresh=meanpdf;
end %if nargin>1

% FInal conversions, if any.
if strcmp(convert,'delta_L')
    % de-logify
	thresh=10^meanpdf;
else
	thresh=meanpdf;
end

if PLOT_pdf
    if ~isnan(response)
        hold on
    end
    figure(1), clf
    h = plot(T,q);
    old_x_scale = get(gca,'xtick');
    new_x_scale = 10.^old_x_scale;
    set(gca,'xticklabel',new_x_scale)
    line([meanpdf meanpdf],get(gca,'YLim'))
end

end