/usr/share/psychtoolbox-3/Quest/QuestBetaAnalysis.m is in psychtoolbox-3-common 3.0.11.20140816.dfsg1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 | 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
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;
|