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;
end
fprintf('Now re-analyzing with both threshold and beta as free parameters. ...\n');
for f=fid
	fprintf(f,'logC 	 +-sd 	 beta	 +-sd	 gamma\n');
end
for i=1:length(q(:))
	betaEstimate(i)=QuestBetaAnalysis1(q(i),fid);
end
return

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;
end
qq=QuestRecompute(q2);

% omit betas that have zero probability
for i=1:length(qq)
	p(i)=sum(qq(i).pdf);
end
if any(p==0)
	fprintf('Omitting beta values ');
	fprintf('%.1f ',qq(find(p==0)).beta);
	fprintf('because they have zero probability.\n');
end
clear q2
q2=qq(find(p));

t2=QuestMean(q2); % estimate threshold for each possible beta
p2=QuestPdf(q2,t2); % get probability of each of these (threshold,beta) combinations
sd2=QuestSd(q2); % get sd of threshold for each possible beta
beta2=[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) combination
t=t2(i); % threshold at that mode
%t=QuestMean(q2(i)); % mean threshold estimate
sd=QuestSd(q2(i)); % sd of threshold estimate at the beta of that mode
p=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);
end
betaEstimate=1/iBetaMean;