% QuestDemo.m
%
% By commenting and uncommenting five lines below, you can use
% this file to implement three QUEST-related procedures for measuring
% threshold.
%
% QuestMode: In the original algorithm of Watson & Pelli (1983)
% each trial and the final estimate are at the MODE of the posterior pdf.
%
% QuestMean: In the improved algorithm of King-Smith et al. (1994).
% each trial and the final estimate are at the MEAN of the posterior pdf.
%
% QuestQuantile & QuestMean: In the ideal algorithm of Pelli (1987)
% each trial is at the best QUANTILE, and the final estimate is at 
% the MEAN of the posterior pdf.
%
% King-Smith, P. E., Grigsby, S. S., Vingrys, A. J., Benes, S. C., and Supowit, A.
% (1994) Efficient and unbiased modifications of the QUEST threshold method: theory, 
% simulations, experimental evaluation and practical implementation. 
% Vision Res, 34 (7), 885-912.
%
% Pelli, D. G. (1987) The ideal psychometric procedure. Investigative Ophthalmology 
% & Visual Science, 28 (Suppl), 366.
%
% Watson, A. B. and Pelli, D. G. (1983) QUEST: a Bayesian adaptive psychometric 
% method. Percept Psychophys, 33 (2), 113-20.

% Copyright (c) 1996-9 Denis G. Pelli
%
% 3/3/97  dhb  Cosmetic editing.

% NOTE: GetSecs.mex is part of the Psychophysics Toolbox and only works on
% Macs.  If you are running QuestDemo on some other machine, use CPUTIME 
% instead of GetSecs.

fprintf('The intensity scale is abstract, but usually we think of it as representing log contrast.\n');

% We'll need this for the simulation.
tActual=input('Specify true threshold of simulated observer: ');

% Get psychometric function parameters.
tGuess=input('Estimate threshold: ');
tGuessSd=2.0; % sd of Gaussian before clipping to specified range
pThreshold=0.82;
beta=3.5;delta=0.01;gamma=0.5;
q=QuestCreate(tGuess,tGuessSd,pThreshold,beta,delta,gamma);

% fprintf('Your initial guess was %g ± %g\n',tGuess,tGuessSd);
% fprintf('Quest''s initial threshold estimate is %g ± %g\n',QuestMean(q),QuestSd(q));

% Simulate a series of trials.
trialsDesired=100;
wrongRight=str2mat('wrong','right');
timeZero=GetSecs;
for k=1:trialsDesired
	% Get recommended level.  Choose your favorite algorithm.
	tTest=QuestQuantile(q);
	%tTest=QuestMean(q);
	%tTest=QuestMode(q);
	
	tTest=tTest+Sample([-0.1,0,0.1]);
	
	% Simulate a trial
	timeSplit=GetSecs; % omit simulation and printing from reported time/trial.
 	response=QuestSimulate(q,tTest,tActual);
 	fprintf('Trial %3d at %4.1f is %s\n',k,tTest,wrongRight(response+1,:));
	timeZero=timeZero+GetSecs-timeSplit;
	
	% Update the pdf
	q=QuestUpdate(q,tTest,response);
end

% Print results of timing.
fprintf('%.0f ms/trial\n',1000*(GetSecs-timeZero)/trialsDesired);

% Get final estimate.
t=QuestMean(q);
sd=QuestSd(q);
fprintf('Mean threshold estimate is %4.2f ± %.2f\n',t,sd);
%t=QuestMode(q);
%fprintf('Mode threshold estimate is %4.2f\n',t);

%fprintf('\nQuest beta analysis. Beta controls the steepness of the Weibull function.\n');
QuestBetaAnalysis(q); % optional
fprintf('Actual parameters of simulated observer:\n');
fprintf('logC	beta	gamma\n');
fprintf('%5.2f	%4.1f	%5.2f\n',tActual,q.beta,q.gamma);
