listening test for DFG (rbo) deleted

applications/ListeningTests/DFG_Ex1_WallVariation/QuestBetaAnalysis.m

-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<2
	fid=1;endfprintf('Now re-analyzing with both threshold and beta as free parameters. ...\n');for f=fid
	fprintf(f,'logC 	 ±sd 	 beta	 ±sd	 gamma\n');endfor i=1:length(q(:))
	betaEstimate(i)=QuestBetaAnalysis1(q(i),fid);endreturn
function betaEstimate=QuestBetaAnalysis1(q,fid)for i=1:16
	q2(i)=q;
	q2(i).beta=2^(i/4);
	q2(i).dim=250;
	q2(i).grain=0.02;endqq=QuestRecompute(q2);
% omit betas that have zero probabilityfor i=1:length(qq)
	p(i)=sum(qq(i).pdf);endif any(p==0)
	fprintf('Omitting beta values ');
	fprintf('%.1f ',qq(find(p==0)).beta);
	fprintf('because they have zero probability.\n');endclear q2q2=qq(find(p));
t2=QuestMean(q2); % estimate threshold for each possible betap2=QuestPdf(q2,t2); % get probability of each of these (threshold,beta) combinationssd2=QuestSd(q2); % get sd of threshold for each possible betabeta2=[q2.beta];% for f=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,'log p');fprintf(f,'	%7.2f',log10(p2));fprintf(f,'\n');% end[p,i]=max(p2); % take mode, i.e. the most probable (threshold,beta) combinationt=t2(i); % threshold at that modesd=QuestSd(q2(i)); % sd of threshold estimate at the beta of that modep=sum(p2);betaMean=sum(p2.*beta2)/p;betaSd=sqrt(sum(p2.*beta2.^2)/p-(sum(p2.*beta2)/p).^2);% beta has a very skewed distribution, with a long tail out to very large value of beta, whereas 1/beta is % more symmetric, with a roughly normal distribution. Thus it is statistically more efficient to estimate the% parameter as 1/average(1/beta) than as average(beta). "iBeta" stands for inverse beta, 1/beta.% The printout takes the conservative approach of basing the mean on 1/beta, but reporting the sd of beta.iBetaMean=sum(p2./beta2)/p;iBetaSd=sqrt(sum(p2./beta2.^2)/p-(sum(p2./beta2)/p).^2);for f=fid
	%	fprintf(f,'Threshold %4.2f ± %.2f; Beta mode %.1f mean %.1f ± %.1f imean 1/%.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;
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestCreate.m

-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<7
	grain=0.01;endif nargin<8
	dim=500;else
	if range<=0
		error('"range" must be greater than zero.')
	end
	dim=range/grain;
	dim=2*ceil(dim/2);	% round up to an even integerendq.updatePdf=1; % boolean: 0 for no, 1 for yesq.warnPdf=1; % booleanq.normalizePdf=1; % boolean. This adds a few ms per call to QuestUpdate, but otherwise the pdf will underflow after about 1000 trials.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);
% THIS CODE WAS IN THE OLD VERSION. I'VE PASTED "q." INTO THE OBVIOUS PLACES.% THIS IS RETAINED SOLELY TO HELP DEBUG ANY BUGS IN THE NEW CODE.% % prepare all the arrays% 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);			% normalize the pdf % 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)); 		% subset that is strictly monotonic% 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]);% % % Best quantileOrder depends only on min and max of psychometric function.% % For 2-interval forced choice, if pL=0.5 and pH=1 then best quantileOrder=0.60% % We write x*log(x+eps) in place of x*log(x) to get zero instead of NAN when x is zero.% 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);
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestMean.m

-function t=QuestMean(q)
% t=QuestMean(q)
%
% Get the mean threshold estimate.
% 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-2002 Denis Pelli

if nargin~=1
	error('Usage: t=QuestMean(q)')
end
if length(q)>1
	t=zeros(size(q));
	for i=1:length(q(:))
		t(i)=QuestMean(q(i));
	end
	return
end

t=q.tGuess+sum(q.pdf.*q.x)/sum(q.pdf);	% mean of our pdf
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestMode.m

-function [t,p]=QuestMode(q)
% [t,p]=QuestMode(q)
%
% "t" is the mode threshold estimate
% "p" is the value of the (unnormalized) pdf at t.
% 
% 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-2004 Denis Pelli

if nargin~=1
	error('Usage: t=QuestMode(q)')
end
if length(q)>1
	t=zeros(size(q));
	for i=1:length(q(:))
		t(i)=QuestMode(q(i));
	end
	return
end

[p,iMode]=max(q.pdf);
t=q.x(iMode)+q.tGuess;
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestP.m

-function p=QuestP(q,x)
% p=QuestP(q,x)
%
% The probability of a correct (or yes) response at intensity x, assuming
% threshold is at x=0.
%
% See Quest.

% 7/25/04   awi     cosmetic (help text layout).

% Copyright (c) 1996-2004 Denis Pelli

if x<q.x2(1)
	p=q.p2(1)elseif x>q.x2(end)
	p=q.p2(end)
else
	p=interp1(q.x2,q.p2,x);
end
if ~isfinite(p)
	q
	error(sprintf('psychometric function %g at %.2g',p,x))
end
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestPdf.m

-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);
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestQuantile.m

-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<2
	quantileOrder=q.quantileOrder;endp=cumsum(q.pdf);if ~isfinite(p(end))
	error('pdf is not finite')endif p(end)==0
	error('pdf is all zero')endindex=find(diff([-1 p])>0);
if length(index)<2
	error(sprintf('pdf has only %g nonzero point(s)',length(index)))endt=q.tGuess+interp1(p(index),q.x(index),quantileOrder*p(end)); % 40 ms
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestRecompute.m

-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.pThreshold
	error(sprintf('psychometric function range [%.2f %.2f] omits %.2f threshold',q.p2(1),q.p2(end),q.pThreshold))endif any(~isfinite(q.p2))
	error('psychometric function p2 is not finite')endindex=find(diff(q.p2)); 		% subset that is strictly monotonicif length(index)<2
	error(sprintf('psychometric function has only %g strictly monotonic point(s)',length(index)))endq.xThreshold=interp1(q.p2(index),q.x2(index),q.pThreshold);if ~isfinite(q.xThreshold)
	q
	error(sprintf('psychometric function has no %.2f threshold',q.pThreshold))endq.p2=q.delta*q.gamma+(1-q.delta)*(1-(1-q.gamma)*exp(-10.^(q.beta*(q.x2+q.xThreshold))));if any(~isfinite(q.p2))
	q
	error('psychometric function p2 is not finite')endq.s2=fliplr([1-q.p2;q.p2]);if ~isfield(q,'intensity') | ~isfield(q,'response')
	q.intensity=[];
	q.response=[];endif any(~isfinite(q.s2(:)))
	error('psychometric function s2 is not finite')end
% Best quantileOrder depends only on min and max of psychometric function.% For 2-interval forced choice, if pL=0.5 and pH=1 then best quantileOrder=0.60% We write x*log(x+eps) in place of x*log(x) to get zero instead of NaN when x is zero.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);
if any(~isfinite(q.pdf))
	error('prior pdf is not finite')end% recompute the pdf from the historical record of trialsfor k=1:length(q.intensity)
	inten=max(-1e10,min(1e10,q.intensity(k))); % make intensity finite
	ii=size(q.pdf,2)+q.i-round((inten-q.tGuess)/q.grain);
	if ii(1)<1
		ii=ii+1-ii(1);
	end
	if ii(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
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestSd.m

-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);
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestSimulate.m

-function response=QuestSimulate(q,tTest,tActual)
% response=QuestSimulate(q,intensity,tActual)
%
% Simulate the response of an observer with threshold tActual.
% 
% See Quest.

% Denis Pelli, 6/8/96
% 3/1/97 dgp restrict intensity parameter to range of x2.
% 3/1/97 dgp updated to use Matlab 5 structs.
% 4/12/99 dgp dropped support for Matlab 4.

% Copyright (c) 1996-2004 Denis Pelli

if nargin~=3
	error('Usage: response=QuestSimulate(q,tTest,tActual)')
end
if length(q)>1
	error('can''t deal with q being a vector')
end
t=min(max(tTest-tActual,q.x2(1)),q.x2(end));
response=interp1(q.x2,q.p2,t) > rand(1);
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/QuestUpdate.m

-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)<1
			ii=ii+1-ii(1);
		else
			ii=ii+size(q.s2,2)-ii(end);
		end
	end
	q.pdf=q.pdf.*q.s2(response+1,ii); % 4 ms
	if q.normalizePdf
		q.pdf=q.pdf/sum(q.pdf);		% keep the pdf normalized	% 3 ms
	endend
% keep a historical record of the trialsq.intensity=[q.intensity,intensity];q.response=[q.response,response];
\ No newline at end of file

applications/ListeningTests/DFG_Ex1_WallVariation/test_rbo_questMean.m

-pThreshold  = 0.75;
-beta        = 3.5;
-delta       = 0.01;
-gamma       = 0.5;
-
-tGuess      = 4;
-tGuessSd    = 3;
-
-q           = QuestCreate(tGuess,tGuessSd,pThreshold,beta,delta,gamma);
-q.normalizePdf = 1; 
-
-trialsDesired=20;
-
-wrongRight={'wrong','right'};
-res         = ones(trialsDesired,1);
-
-for k=1:trialsDesired
-	% Get recommended level.  Choose your favorite algorithm.
-%	tTest=QuestQuantile(q);	% Recommended by Pelli (1987), and still our favorite.
- 	tTest=QuestMean(q);		% Recommended by King-Smith et al. (1994)
-% 	tTest=QuestMode(q);		% Recommended by Watson & Pelli (1983)
-	
-%   HERE IS MY CODE
-    %response = round(rand(1,1));
-
-    response  = res(k);
-    
-    q=QuestUpdate(q,tTest,response); % Add the new datum (actual test intensity and observer response) to the database.
-    disp(num2str(tTest))
-end
-
-% Ask Quest for the final estimate of threshold.
-t=QuestMean(q);		% Recommended by Pelli (1989) and King-Smith et al. (1994). Still our favorite.
-sd=QuestSd(q);
-fprintf('Final threshold estimate (meansd) is %.2f  %.2f\n',t,sd);
-
-% Optionally, reanalyze the data with beta as a free parameter.
-QuestBetaAnalysis(q); % optional
-fprintf('Actual parameters of simulated observer:\n');
-fprintf('logC	beta	gamma\n');
-