function QuestBetaAnalysis(q,fid)
% QuestBetaAnalysis(q,[fid])
% Analyzes the quest function with beta as a free parameter. It returns the
% mean estimates of alpha (as logC) and beta. Gamma is left at whatever value
% the user fixed it at.
% 
% Denis Pelli 5/6/99
% 8/23/99 dgp streamlined the printout
% 8/24/99 dgp add sd's to printout
if nargin<1 | nargin>2
	error('Usage: QuestBetaAnalysis(q,[fid])')
end
if nargin<2
	fid=1;
end
fprintf('Now re-analyzing with beta as a free parameter. . . .\n');
for f=fid
%	fprintf(f,'logC	beta	gamma\n');
	fprintf(f,'logC 	 ±sd 	 beta	 ±sd	 gamma\n');
end
for i=1:length(q(:))
	QuestBetaAnalysis1(q(i),fid);
end
return

function 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
q2=QuestRecompute(q2);
t2=QuestMean(q2);
p2=QuestPdf(q2,t2);
sd2=QuestSd(q2);
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);
t=t2(i);
sd=QuestSd(q2(i));
p=sum(p2);
betaMean=sum(p2.*beta2)/p;
betaSd=sqrt(sum(p2.*beta2.^2)/p-(sum(p2.*beta2)/p).^2);
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
