function betaEstimate=QuestBetaAnalysis(q,fid)
% betaEstimate=QuestBetaAnalysis(q,[fid]);
%
% Analyzes the quest function with beta as a free parameter. It prints (in
% the file or files pointed to by fid) the mean estimates of alpha (as
% logC) and beta. Gamma is left at whatever value the user fixed it at.
%
% Note that normalization of the pdf, by QuestRecompute, is disabled because it
% would need to be done across the whole q vector. Without normalization,
% the pdf tends to underflow at around 1000 trials. You will have some warning
% of this because the printout mentions any values of beta that were dropped
% because they had zero probability. Thus you should keep the number of trials
% under around 1000, to avoid the zero-probability warnings.
%
% See Quest.
% Denis Pelli 5/6/99
% 8/23/99 dgp streamlined the printout
% 8/24/99 dgp add sd's to printout
% 10/13/04 dgp added comment explaining 1/beta
if nargin<1|nargin>2
error('Usage: QuestBetaAnalysis(q,[fid])')
end
if nargin<2fid=1;endfprintf('Nowre-analyzingwithboththresholdandbetaasfreeparameters....\n');forf=fidfprintf(f,'logC±sdbeta±sdgamma\n');endfori=1:length(q(:))betaEstimate(i)=QuestBetaAnalysis1(q(i),fid);endreturnfunctionbetaEstimate=QuestBetaAnalysis1(q,fid)fori=1:16q2(i)=q;q2(i).beta=2^(i/4);q2(i).dim=250;q2(i).grain=0.02;endqq=QuestRecompute(q2);%omitbetasthathavezeroprobabilityfori=1:length(qq)p(i)=sum(qq(i).pdf);endifany(p==0)fprintf('Omittingbetavalues');fprintf('%.1f',qq(find(p==0)).beta);fprintf('becausetheyhavezeroprobability.\n');endclearq2q2=qq(find(p));t2=QuestMean(q2);%estimatethresholdforeachpossiblebetap2=QuestPdf(q2,t2);%getprobabilityofeachofthese(threshold,beta)combinationssd2=QuestSd(q2);%getsdofthresholdforeachpossiblebetabeta2=[q2.beta];%forf=fid%fprintf(f,'beta');fprintf(f,'%7.1f',q2(:).beta);fprintf(f,'\n');%fprintf(f,'t');fprintf(f,'%7.2f',t2);fprintf(f,'\n');%fprintf(f,'sd');fprintf(f,'%7.2f',sd2);fprintf(f,'\n');%fprintf(f,'logp');fprintf(f,'%7.2f',log10(p2));fprintf(f,'\n');%end[p,i]=max(p2);%takemode,i.e.themostprobable(threshold,beta)combinationt=t2(i);%thresholdatthatmodesd=QuestSd(q2(i));%sdofthresholdestimateatthebetaofthatmodep=sum(p2);betaMean=sum(p2.*beta2)/p;betaSd=sqrt(sum(p2.*beta2.^2)/p-(sum(p2.*beta2)/p).^2);%betahasaveryskeweddistribution,withalongtailouttoverylargevalueofbeta,whereas1/betais%moresymmetric,witharoughlynormaldistribution.Thusitisstatisticallymoreefficienttoestimatethe%parameteras1/average(1/beta)thanasaverage(beta)."iBeta"standsforinversebeta,1/beta.%Theprintouttakestheconservativeapproachofbasingthemeanon1/beta,butreportingthesdofbeta.iBetaMean=sum(p2./beta2)/p;iBetaSd=sqrt(sum(p2./beta2.^2)/p-(sum(p2./beta2)/p).^2);forf=fid%fprintf(f,'Threshold%4.2f±%.2f;Betamode%.1fmean%.1f±%.1fimean1/%.1f±%.1f;Gamma%.2f\n',t,sd,q2(i).beta,betaMean,betaSd,1/iBetaMean,iBetaSd,q.gamma);%fprintf(f,'%5.2f%4.1f%5.2f\n',t,1/iBetaMean,q.gamma);fprintf(f,'%5.2f%5.2f%4.1f%4.1f%6.3f\n',t,sd,1/iBetaMean,betaSd,q.gamma);endbetaEstimate=1/iBetaMean;

function q=QuestCreate(tGuess,tGuessSd,pThreshold,beta,delta,gamma,grain,range)
% q=QuestCreate(tGuess,tGuessSd,pThreshold,beta,delta,gamma,[grain],[range])
%
% Create a struct q with all the information necessary to measure
% threshold. Threshold "t" is measured on an abstract "intensity"
% scale, which usually corresponds to log10 contrast.
%
% QuestCreate saves in struct q the parameters for a Weibull psychometric function:
% p2=delta*gamma+(1-delta)*(1-(1-gamma)*exp(-10.^(beta*(x2+xThreshold))));
% where x represents log10 contrast relative to threshold. The Weibull
% function itself appears only in QuestRecompute, which uses the
% specified parameter values in q to compute a psychometric function
% and store it in q. All the other Quest functions simply use the
% psychometric function stored in "q". QuestRecompute is called solely
% by QuestCreate and QuestBetaAnalysis (and possibly by a few user
% programs). Thus, if you prefer to use a different kind of
% psychometric function, called Foo, you need only create your own
% QuestCreateFoo, QuestRecomputeFoo, and (if you need it)
% QuestBetaAnalysisFoo, based on QuestCreate, QuestRecompute, and
% QuestBetaAnalysis, and you can use them with the rest of the Quest
% package unchanged. You would only be changing a few lines of code,
% so it would quite easy to do.
%
% Several users of Quest have asked questions on the Psychtoolbox forum
% about how to restrict themselves to a practical testing range. That is
% not what tGuessSd and "range" are for; they should be large, e.g. I
% typically set tGuessSd=3 and range=5 when intensity represents log
% contrast. If necessary, you should restrict the range yourself, outside
% of Quest. Here, in QuestCreate, you tell Quest about your prior beliefs,
% and you should try to be open-minded, giving Quest a generously large
% range to consider as possible values of threshold. For each trial you
% will later ask Quest to suggest a test intensity. It is important to
% realize that what Quest returns is just what you asked for, a
% suggestion. You should then test at whatever intensity you like, taking
% into account both the suggestion and any practical constraints (e.g. a
% maximum and minimum contrast that you can achieve, and quantization of
% contrast). After running the trial you should call QuestUpdate with the
% contrast that you actually used and the observer's response to add your
% new datum to the database. Don't restrict "tGuessSd" or "range" by the
% limitations of what you can display. Keep open the possibility that
% threshold may lie outside the range of contrasts that you can produce,
% and let Quest consider all possibilities.
%
% There is one exception to the above advice of always being generous with
% tGuessSd. Occasionally we find that we have a working Quest-based
% program that measures threshold, and we discover that we need to measure
% the proportion correct at a particular intensity. Instead of writing a
% new program, or modifying the old one, it is often more convenient to
% instead reduce tGuessSd to practically zero, e.g. a value like 0.001,
% which has the effect of restricting all threshold estimates to be
% practically identical to tGuess, making it easy to run any number of
% trials at that intensity. Of course, in this case, the final threshold
% estimate from Quest should be ignored, since it is merely parroting back
% to you the assertion that threshold is equal to the initial guess
% "tGuess". What's of interest is the final proportion correct; at the
% end, call QuestTrials or add an FPRINTF statement to report it.
%
% tGuess is your prior threshold estimate.
% tGuessSd is the standard deviation you assign to that guess. Be generous.
% pThreshold is your threshold criterion expressed as probability of
% response==1. An intensity offset is introduced into the psychometric
% function so that threshold (i.e. the midpoint of the table) yields
% pThreshold.
% beta, delta, and gamma are the parameters of a Weibull psychometric function.
% beta controls the steepness of the psychometric function. Typically 3.5.
% delta is the fraction of trials on which the observer presses blindly.
% Typically 0.01.
% gamma is the fraction of trials that will generate response 1 when
% intensity==-inf.
% grain is the quantization (step size) of the internal table. E.g. 0.01.
% range is the intensity difference between the largest and smallest
% intensity that the internal table can store. E.g. 5. This interval will
% be centered on the initial guess tGuess, i.e.
% tGuess+(-range/2:grain:range/2). "range" is used only momentarily here,
% to determine "dim", which is retained in the quest struct. "dim" is the
% number of distinct intensities that the internal table can store, e.g.
% 500. QUEST assumes that intensities outside of this interval have zero
% prior probability, i.e. they are impossible values for threshold. The
% cost of making "range" too big is some extra storage and computation,
% which are usually negligible. The cost of making "range" too small is
% that you prejudicially exclude what are actually possible values for
% threshold. Getting out-of-range warnings from QuestUpdate is one
% possible indication that your stated range is too small.
%
% See Quest.
% 6/8/96 dgp Wrote it.
% 6/11/96 dgp Optimized the order of stuffing for faster unstuffing.
% 11/10/96 dhb Added warning about correctness after DGP told me.
% 3/1/97 dgp Fixed error in sign of xThreshold in formula for p2.
% 3/1/97 dgp Updated to use Matlab 5 structs.
% 3/3/97 dhb Added missing semicolon to first struct eval.
% 3/5/97 dgp Fixed sd: use exp instead of 10^.
% 3/5/97 dgp Added some explanation of the psychometric function.
% 6/24/97 dgp For simulations, now allow specification of grain and dim.
% 9/30/98 dgp Added "dim" fix from Richard Murray.
% 4/12/99 dgp dropped support for Matlab 4.
% 5/6/99 dgp Simplified "dim" calculation; just round up to even integer.
% 8/15/99 dgp Explain how to use other kind of psychometric function.
% 2/10/02 dgp Document grain and range.
% 9/11/04 dgp Explain why supplied "range" should err on the high side.
% 10/13/04 dgp Explain why tGuesSd and range should be large, generous.
% 10/13/04 dgp Set q.normalizePdf to 1, to avoid underflow errors that otherwise accur after around 1000 trials.
%
% Copyright (c) 1996-2004 Denis Pelli
if nargin<6|nargin>8
error('Usage: q=QuestCreate(tGuess,tGuessSd,pThreshold,beta,delta,gamma,[grain],[range])')
end
if nargin<7grain=0.01;endifnargin<8dim=500;elseifrange<=0error('"range"mustbegreaterthanzero.')enddim=range/grain;dim=2*ceil(dim/2);%rounduptoanevenintegerendq.updatePdf=1;%boolean:0forno,1foryesq.warnPdf=1;%booleanq.normalizePdf=1;%boolean.ThisaddsafewmspercalltoQuestUpdate,butotherwisethepdfwillunderflowafterabout1000trials.q.tGuess=tGuess;q.tGuessSd=tGuessSd;q.pThreshold=pThreshold;q.beta=beta;q.delta=delta;q.gamma=gamma;q.grain=grain;q.dim=dim;q=QuestRecompute(q);%THISCODEWASINTHEOLDVERSION.I'VEPASTED"q."INTOTHEOBVIOUSPLACES.%THISISRETAINEDSOLELYTOHELPDEBUGANYBUGSINTHENEWCODE.%%prepareallthearrays%q.i=-dim/2:dim/2;%q.x=i*grain;%q.pdf=exp(-0.5*(q.x/tGuessSd).^2);%q.pdf=q.pdf/sum(q.pdf);%normalizethepdf%i2=-dim:dim;%q.x2=i2*q.grain;%q.p2=delta*gamma+(1-delta)*(1-(1-gamma)*exp(-10.^(beta*q.x2)));%index=find(diff(q.p2));%subsetthatisstrictlymonotonic%q.xThreshold=interp1(q.p2(index),q.x2(index),q.pThreshold);%q.p2=delta*gamma+(1-delta)*(1-(1-gamma)*exp(-10.^(beta*(q.x2+q.xThreshold))));%q.s2=fliplr([1-q.p2;q.p2]);%%%BestquantileOrderdependsonlyonminandmaxofpsychometricfunction.%%For2-intervalforcedchoice,ifpL=0.5andpH=1thenbestquantileOrder=0.60%%Wewritex*log(x+eps)inplaceofx*log(x)togetzeroinsteadofNANwhenxiszero.%pL=q.p2(1);%pH=q.p2(end);%pE=pH*log(pH+eps)-pL*log(pL+eps)+(1-pH+eps)*log(1-pH+eps)-(1-pL+eps)*log(1-pL+eps);%pE=1/(1+exp(pE/(pL-pH)));%q.quantileOrder=(pE-pL)/(pH-pL);

function p=QuestPdf(q,t)
% p=QuestPdf(q,t)
%
% The (possibly unnormalized) probability density of candidate threshold "t".
% q and t may be vectors of the same size, in which case the returned p is a vector of that size.
%
% See Quest.
% Denis Pelli
% 5/6/99 dgp wrote it
% Copyright (c) 1996-1999 Denis Pelli
if nargin~=2
error('Usage: p=QuestPdf(q,t)')
end
if size(q)~=size(t)
error('both arguments must have the same dimensions')
end
if length(q)>1
p=zeros(size(q));
for i=1:length(q(:))
p(i)=QuestPdf(q(i),t(i));
end
return
end
i=round((t-q.tGuess)/q.grain)+1+q.dim/2;
i=min(length(q.pdf),max(1,i));
p=q.pdf(i);

function t=QuestQuantile(q,quantileOrder)
% intensity=QuestQuantile(q,[quantileOrder])
%
% Gets a quantile of the pdf in the struct q. You may specify the desired
% quantileOrder, e.g. 0.5 for median, or, making two calls, 0.05 and 0.95
% for a 90% confidence interval. If the "quantileOrder" argument is not
% supplied, then it's taken from the "q" struct. QuestCreate uses
% QuestRecompute to compute the optimal quantileOrder and saves that in the
% "q" struct; this quantileOrder yields a quantile that is the most
% informative intensity for the next trial.
%
% This is based on work presented at a conference, but otherwise unpublished:
% Pelli, D. G. (1987). The ideal psychometric procedure. Investigative
% Ophthalmology & Visual Science, 28(Suppl), 366.
%
% See Quest.
% Denis Pelli, 6/9/96
% 6/17/96 dgp, worked around "nonmonotonic" (i.e. not strictly monotonic)
% interp1 error.
% 3/1/97 dgp updated to use Matlab 5 structs.
% 4/12/99 dgp removed support for Matlab 4.
%
% Copyright (c) 1996-1999 Denis Pelli
if nargin>2
error('Usage: intensity=QuestQuantile(q,[quantileOrder])')
end
if length(q)>1
if nargin>1
error('Can''t accept quantileOrder for q vector. Set each q.quantileOrder instead.')
end
t=zeros(size(q));
for i=1:length(q(:))
t(i)=QuestQuantile(q(i));
end
return
end
if nargin<2quantileOrder=q.quantileOrder;endp=cumsum(q.pdf);if~isfinite(p(end))error('pdfisnotfinite')endifp(end)==0error('pdfisallzero')endindex=find(diff([-1p])>0);
if length(index)<2error(sprintf('pdfhasonly%gnonzeropoint(s)',length(index)))endt=q.tGuess+interp1(p(index),q.x(index),quantileOrder*p(end));%40ms

function q=QuestRecompute(q)
% q=QuestRecompute(q)
%
% Call this immediately after changing a parameter of the psychometric
% function. QuestRecompute uses the specified parameters in "q" to
% recompute the psychometric function. It then uses the newly computed
% psychometric function and the history in q.intensity and q.response
% to recompute the pdf. (QuestRecompute does nothing if q.updatePdf is
% false.)
%
% QuestCreate saves in struct q the parameters for a Weibull psychometric function:
% p2=delta*gamma+(1-delta)*(1-(1-gamma)*exp(-10.^(beta*(x2+xThreshold))));
% where x represents log10 contrast relative to threshold. The Weibull
% function itself appears only in QuestRecompute, which uses the
% specified parameter values in q to compute a psychometric function
% and store it in q. All the other Quest functions simply use the
% psychometric function stored in "q". QuestRecompute is called solely
% by QuestCreate and QuestBetaAnalysis (and possibly by a few user
% programs). Thus, if you prefer to use a different kind of
% psychometric function, called Foo, you need only create your own
% QuestCreateFoo, QuestRecomputeFoo, and (if you need it)
% QuestBetaAnalysisFoo, based on QuestCreate, QuestRecompute, and
% QuestBetaAnalysis, and you can use them with the rest of the Quest
% package unchanged. You would only be changing a few lines of code,
% so it would quite easy to do.
%
% "dim" is the number of distinct intensities that the internal tables in q can store,
% e.g. 500. This vector, of length "dim", with increment size "grain",
% will be centered on the initial guess tGuess, i.e.
% tGuess+[-range/2:grain:range/2]. QUEST assumes that intensities outside
% of this interval have zero prior probability, i.e. they are impossible
% values for threshold. The cost of making "dim" too big is some extra
% storage and computation, which are usually negligible. The cost of
% making "dim" too small is that you prejudicially exclude what are
% actually possible values for threshold. Getting out-of-range warnings
% from QuestUpdate is one possible indication that your stated range is
% too small.
%
% See QuestCreate, QuestUpdate, QuestQuantile, QuestMean, QuestMode,
% QuestSd, and QuestSimulate.
% 4/29/99 dgp Wrote it.
% 8/15/99 dgp Explain how to use other kind of psychometric function.
% 9/11/04 dgp Explain why supplied "dim" should err on the high side.
% Copyright (c) 1996-2004 Denis Pelli
if nargin~=1
error('Usage: q=QuestRecompute(q)')
end
if length(q)>1
for i=1:length(q(:))
q(i).normalizePdf=0; % any norming must be done across the whole set of pdfs, because it's actually one big multi-dimensional pdf.
q(i)=QuestRecompute(q(i));
end
return
end
if ~q.updatePdf
return
end
if q.gamma>q.pThreshold
warning(sprintf('reducing gamma from %.2f to 0.5',q.gamma))
q.gamma=0.5;
end
% prepare all the arrays
q.i=-q.dim/2:q.dim/2;
q.x=q.i*q.grain;
q.pdf=exp(-0.5*(q.x/q.tGuessSd).^2);
q.pdf=q.pdf/sum(q.pdf);
i2=-q.dim:q.dim;
q.x2=i2*q.grain;
q.p2=q.delta*q.gamma+(1-q.delta)*(1-(1-q.gamma)*exp(-10.^(q.beta*q.x2)));
if q.p2(1)>=q.pThreshold | q.p2(end)<=q.pThresholderror(sprintf('psychometricfunctionrange[%.2f%.2f]omits%.2fthreshold',q.p2(1),q.p2(end),q.pThreshold))endifany(~isfinite(q.p2))error('psychometricfunctionp2isnotfinite')endindex=find(diff(q.p2));%subsetthatisstrictlymonotoniciflength(index)<2error(sprintf('psychometricfunctionhasonly%gstrictlymonotonicpoint(s)',length(index)))endq.xThreshold=interp1(q.p2(index),q.x2(index),q.pThreshold);if~isfinite(q.xThreshold)qerror(sprintf('psychometricfunctionhasno%.2fthreshold',q.pThreshold))endq.p2=q.delta*q.gamma+(1-q.delta)*(1-(1-q.gamma)*exp(-10.^(q.beta*(q.x2+q.xThreshold))));ifany(~isfinite(q.p2))qerror('psychometricfunctionp2isnotfinite')endq.s2=fliplr([1-q.p2;q.p2]);if~isfield(q,'intensity')|~isfield(q,'response')q.intensity=[];q.response=[];endifany(~isfinite(q.s2(:)))error('psychometricfunctions2isnotfinite')end%BestquantileOrderdependsonlyonminandmaxofpsychometricfunction.%For2-intervalforcedchoice,ifpL=0.5andpH=1thenbestquantileOrder=0.60%Wewritex*log(x+eps)inplaceofx*log(x)togetzeroinsteadofNaNwhenxiszero.pL=q.p2(1);pH=q.p2(size(q.p2,2));pE=pH*log(pH+eps)-pL*log(pL+eps)+(1-pH+eps)*log(1-pH+eps)-(1-pL+eps)*log(1-pL+eps);pE=1/(1+exp(pE/(pL-pH)));q.quantileOrder=(pE-pL)/(pH-pL);ifany(~isfinite(q.pdf))error('priorpdfisnotfinite')end%recomputethepdffromthehistoricalrecordoftrialsfork=1:length(q.intensity)inten=max(-1e10,min(1e10,q.intensity(k)));%makeintensityfiniteii=size(q.pdf,2)+q.i-round((inten-q.tGuess)/q.grain);ifii(1)<1ii=ii+1-ii(1);endifii(end)>size(q.s2,2)
ii=ii+size(q.s2,2)-ii(end);
end
q.pdf=q.pdf.*q.s2(q.response(k)+1,ii); % 4 ms
if q.normalizePdf & mod(k,100)==0
q.pdf=q.pdf/sum(q.pdf); % avoid underflow; keep the pdf normalized % 3 ms
end
end
if q.normalizePdf
q.pdf=q.pdf/sum(q.pdf); % keep the pdf normalized % 3 ms
end
if any(~isfinite(q.pdf))
error('pdf is not finite')
end

function sd=QuestSd(q)
% sd=QuestSd(q)
%
% Get the sd of the threshold distribution.
% If q is a vector, then the returned t is a vector of the same size.
%
% See Quest.
% Denis Pelli, 6/8/96
% 3/1/97 dgp updated to use Matlab 5 structs.
% 4/12/99 dgp dropped support for Matlab 4.
% Copyright (c) 1996-1999 Denis Pelli
if nargin~=1
error('Usage: sd=QuestSd(q)')
end
if length(q)>1
sd=zeros(size(q));
for i=1:length(q(:))
sd(i)=QuestSd(q(i));
end
return
end
p=sum(q.pdf);
sd=sqrt(sum(q.pdf.*q.x.^2)/p-(sum(q.pdf.*q.x)/p).^2);

function q=QuestUpdate(q,intensity,response)
% q=QuestUpdate(q,intensity,response)
%
% Update the struct q to reflect the results of this trial. The historical
% records q.intensity and q.response are always updated, but q.pdf is only
% updated if q.updatePdf is true. You can always call QuestRecompute to
% recreate q.pdf from scratch from the historical record.
%
% See Quest.
% Denis Pelli, 6/11/96
% 2/28/97 dgp Updated for Matlab 5: call round.
% 4/12/99 dgp Dropped support for Matlab 4.
% 4/30/99 dgp Give explanatory error message when intensity is out of bounds.
% 9/11/04 dgp For compatibility with Matlab 6.5, comment out the testing
% of WARNING level.
%
% Copyright (c) 1996-2004 Denis Pelli
if nargin~=3
error('Usage: q=QuestUpdate(q,intensity,response)')
end
if length(q)>1
error('can''t deal with q being a vector')
end
if response<0|response>=size(q.s2,1)
error(sprintf('response %g out of range 0 to %d',response,size(q.s2,1)-1))
end
if q.updatePdf
inten=max(-1e10,min(1e10,intensity)); % make intensity finite
ii=size(q.pdf,2)+q.i-round((inten-q.tGuess)/q.grain);
if ii(1)<1|ii(end)>size(q.s2,2)
if q.warnPdf
low=(1-size(q.pdf,2)-q.i(1))*q.grain+q.tGuess;
high=(size(q.s2,2)-size(q.pdf,2)-q.i(end))*q.grain+q.tGuess;
oldWarning=warning;
warning('on'); % no backtrace
warning(sprintf('QuestUpdate: intensity %.2f out of range %.2f to %.2f. Pdf will be inexact. Suggest that you increase "range" in call to QuestCreate.',intensity,low,high));
warning(oldWarning);
end
if ii(1)<1ii=ii+1-ii(1);elseii=ii+size(q.s2,2)-ii(end);endendq.pdf=q.pdf.*q.s2(response+1,ii);%4msifq.normalizePdfq.pdf=q.pdf/sum(q.pdf);%keepthepdfnormalized%3msendend%keepahistoricalrecordofthetrialsq.intensity=[q.intensity,intensity];q.response=[q.response,response];