QuestBetaAnalysis.m 3.21 KB
Newer Older
1
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;