Estimate multiple psychometric functions
Usage: [xx,ii,m,v]=psycest(-n,p,q,xp) % initialize n models
[xx,ii,m,v]=psycest(i,x,r) % supply a trial result to psycest
psycest(i) % plot pdf of model i
[p,q]=psychest(0) % output model parameters (or print them if no outputs)
Inputs:
-n minus the number of models
p(:,n) parameters for each model
1 thresh [0.75]
2 miss [0.04]
3 guess [0.1]
4 SNR min [-20]
5 SNR max [20]
6 Slope min [0]
7 slope max [0.5]
q(:) parameters common to all models (vector or struct)
1 q.nx number of SNR values in pdf [40]
2 q.ns number of slope values in pdf [21]
3 q.nh number of probe SNR values to evaluate [30]
4 q.cs weighting of slope relative to threshold in cost function [1]
5 q.dh minimum step size in dB for probe SNRs [0.2]
6 q.sl min slope at threshold [0.02]
7 q.kp number of std deviations of the pdfs to keep [4]
8 q.hg amount to grow expected gains in ni trials [1.3]
9 q.cf cost function: 1=variance, 2=entropy [2]
10 q.pm psychometric model: 1=logistic, 2=cumulative gaussian [1]
11 q.lg use log slope in pdf: 0=no, 1=yes [1]
xp{n}(:) list of available probe SNR values
i model probed
x probe SNR value used
r response: 0=wrong, 1=right.
Outputs:
xx recommended probe SNR
ii recommended model to probe next
m(2,n,3) estimated srt and slope of all models
m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
v(3,n) estimated covariance matrix entries:
[var(srt) cov(srt,slope) var(slope)]'
Algorithm parameters:
The algorithm estimates the psychometric function for a number of models simultaneously.
A psychometric function specifies the probability of correct recognition as a function of
SNR and is completely specified by two parameters: the SNR (in dB) and the slope (in 1/dB)
at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify
some of the details of the psychometric function and can be different for each model:
p(1,n) gives the target recognition probability
p(2,n) gives the probability of incorrect recognition at very good SNRs (the "miss" probability).
If this value is made too low, a single unwarrented recognition error will have a large
effect of the estimated parameters and greatly lengthen the time to convergence. The default
value is 4%.
p(3,n) gives the probability of correct recognition at very poor SNRs (the "guess" probabiity).
This should be set to 1/N where N is the number of possible responses to a stimulus.
p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will
be adjusted by the algorithm if necessary but wrong values will delay convergence.
p(6:7,n) These given the initial min and max slope (in probability per dB) at the target
recognition probability.
The remaining parameters are shared between all the models and control how the algorithm operates.
q(1:2) These parameters determine the sampling density of the joint pdf of SNR and Slope.
q(3) This parameter specifies how many potential values of probe SNR to evaluate in order to
determine which is likely to give the most improvement in the parameter estimates.
q(4) This parameter gives the relative weight of SNR and Slope in the cost function. Increasing
its value will improve the slope estimate (or log(slope) estimate) at the expense of the
SNR estimate. Actually its value is not all that critical.
q(5) At each iteration psycest evaluates several potential probe SNR levels to see which is
likely to give the most improvement to the parameter estimates. This parameter fixes the
minimum spacing between these potential probe values. It only has an effect when the variances
of the parameter estimates are very small.
q(6) To prevent the routine getting lost, this parameter specifies the smallest reasonable value
of the Slope parameter at the target recognition probability. Under normal circumstances, it
has no effect.
q(7) For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope.
The sampling grid is reclculated after each iteration to encompass this number of standard
deviations in each dimension.
q(8) At each iteration, psycest advises which model to probe next on the basis of which
gives the greatest expected reduction in cost function. To ensure that all models
are periodically sampled, this expected reduction of each unprobed model is multiplied
by this parameter at each iteration. Thus a model will eventually be probed even if
its expected cost factor improvement is small.
q(9) This determines whether the recommended probe SNR to use at the next iteration is chosen to
minimize the expected variance or the entropy of the estimated SNR and slope distributions.
My experience is that entropy (the default) gives faster convergence.
q(10) This selects whether the underlying model is a logistic function (1) or a cumulative
gaussian (2). The choice makes little difference unless the "miss" or "guess" probabilities
are very very small.
q(11) This selects whether the internal pdf samples "slope" or "log(slope)". It is recommended
that you stick to the default of log(slope) since this results in more symmetrical
distributions.0001 function [xx,ii,m,v]=psycest(iq,x,r,xp) 0002 % Estimate multiple psychometric functions 0003 % 0004 % Usage: [xx,ii,m,v]=psycest(-n,p,q,xp) % initialize n models 0005 % [xx,ii,m,v]=psycest(i,x,r) % supply a trial result to psycest 0006 % psycest(i) % plot pdf of model i 0007 % [p,q]=psychest(0) % output model parameters (or print them if no outputs) 0008 % 0009 % Inputs: 0010 % -n minus the number of models 0011 % p(:,n) parameters for each model 0012 % 1 thresh [0.75] 0013 % 2 miss [0.04] 0014 % 3 guess [0.1] 0015 % 4 SNR min [-20] 0016 % 5 SNR max [20] 0017 % 6 Slope min [0] 0018 % 7 slope max [0.5] 0019 % q(:) parameters common to all models (vector or struct) 0020 % 1 q.nx number of SNR values in pdf [40] 0021 % 2 q.ns number of slope values in pdf [21] 0022 % 3 q.nh number of probe SNR values to evaluate [30] 0023 % 4 q.cs weighting of slope relative to threshold in cost function [1] 0024 % 5 q.dh minimum step size in dB for probe SNRs [0.2] 0025 % 6 q.sl min slope at threshold [0.02] 0026 % 7 q.kp number of std deviations of the pdfs to keep [4] 0027 % 8 q.hg amount to grow expected gains in ni trials [1.3] 0028 % 9 q.cf cost function: 1=variance, 2=entropy [2] 0029 % 10 q.pm psychometric model: 1=logistic, 2=cumulative gaussian [1] 0030 % 11 q.lg use log slope in pdf: 0=no, 1=yes [1] 0031 % xp{n}(:) list of available probe SNR values 0032 % i model probed 0033 % x probe SNR value used 0034 % r response: 0=wrong, 1=right. 0035 % 0036 % Outputs: 0037 % xx recommended probe SNR 0038 % ii recommended model to probe next 0039 % m(2,n,3) estimated srt and slope of all models 0040 % m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates 0041 % v(3,n) estimated covariance matrix entries: 0042 % [var(srt) cov(srt,slope) var(slope)]' 0043 % 0044 % Algorithm parameters: 0045 % 0046 % The algorithm estimates the psychometric function for a number of models simultaneously. 0047 % A psychometric function specifies the probability of correct recognition as a function of 0048 % SNR and is completely specified by two parameters: the SNR (in dB) and the slope (in 1/dB) 0049 % at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify 0050 % some of the details of the psychometric function and can be different for each model: 0051 % p(1,n) gives the target recognition probability 0052 % p(2,n) gives the probability of incorrect recognition at very good SNRs (the "miss" probability). 0053 % If this value is made too low, a single unwarrented recognition error will have a large 0054 % effect of the estimated parameters and greatly lengthen the time to convergence. The default 0055 % value is 4%. 0056 % p(3,n) gives the probability of correct recognition at very poor SNRs (the "guess" probabiity). 0057 % This should be set to 1/N where N is the number of possible responses to a stimulus. 0058 % p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will 0059 % be adjusted by the algorithm if necessary but wrong values will delay convergence. 0060 % p(6:7,n) These given the initial min and max slope (in probability per dB) at the target 0061 % recognition probability. 0062 % The remaining parameters are shared between all the models and control how the algorithm operates. 0063 % q(1:2) These parameters determine the sampling density of the joint pdf of SNR and Slope. 0064 % q(3) This parameter specifies how many potential values of probe SNR to evaluate in order to 0065 % determine which is likely to give the most improvement in the parameter estimates. 0066 % q(4) This parameter gives the relative weight of SNR and Slope in the cost function. Increasing 0067 % its value will improve the slope estimate (or log(slope) estimate) at the expense of the 0068 % SNR estimate. Actually its value is not all that critical. 0069 % q(5) At each iteration psycest evaluates several potential probe SNR levels to see which is 0070 % likely to give the most improvement to the parameter estimates. This parameter fixes the 0071 % minimum spacing between these potential probe values. It only has an effect when the variances 0072 % of the parameter estimates are very small. 0073 % q(6) To prevent the routine getting lost, this parameter specifies the smallest reasonable value 0074 % of the Slope parameter at the target recognition probability. Under normal circumstances, it 0075 % has no effect. 0076 % q(7) For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope. 0077 % The sampling grid is reclculated after each iteration to encompass this number of standard 0078 % deviations in each dimension. 0079 % q(8) At each iteration, psycest advises which model to probe next on the basis of which 0080 % gives the greatest expected reduction in cost function. To ensure that all models 0081 % are periodically sampled, this expected reduction of each unprobed model is multiplied 0082 % by this parameter at each iteration. Thus a model will eventually be probed even if 0083 % its expected cost factor improvement is small. 0084 % q(9) This determines whether the recommended probe SNR to use at the next iteration is chosen to 0085 % minimize the expected variance or the entropy of the estimated SNR and slope distributions. 0086 % My experience is that entropy (the default) gives faster convergence. 0087 % q(10) This selects whether the underlying model is a logistic function (1) or a cumulative 0088 % gaussian (2). The choice makes little difference unless the "miss" or "guess" probabilities 0089 % are very very small. 0090 % q(11) This selects whether the internal pdf samples "slope" or "log(slope)". It is recommended 0091 % that you stick to the default of log(slope) since this results in more symmetrical 0092 % distributions. 0093 0094 % Copyright (C) Mike Brookes 2009-2010 0095 % Version: $Id: psycest.m 713 2011-10-16 14:45:43Z dmb $ 0096 % 0097 % VOICEBOX is a MATLAB toolbox for speech processing. 0098 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0099 % 0100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0101 % This program is free software; you can redistribute it and/or modify 0102 % it under the terms of the GNU General Public License as published by 0103 % the Free Software Foundation; either version 2 of the License, or 0104 % (at your option) any later version. 0105 % 0106 % This program is distributed in the hope that it will be useful, 0107 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0108 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0109 % GNU General Public License for more details. 0110 % 0111 % You can obtain a copy of the GNU General Public License from 0112 % http://www.gnu.org/copyleft/gpl.html or by writing to 0113 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0115 0116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0117 % Bugs/Suggestions: 0118 % (1) allow lookahead by 2 tests rather than 1 0119 % (2)use a structure for input parameters 0120 % (3)should only resample the pdfs when necessary 0121 % (4)add a forgetting factor for old measurements 0122 % (8) use quadratic interpolation to find cost function minimum 0123 % (10) optionally output the whole pdf + its axis values 0124 % (11) optionally output all model probe values 0125 % (13) Should perhaps estimate and return the mean and slope compensated for the guess rate 0126 % i.e. if you want p, you actually test at guess+ p*(1-guess-miss)/(1-miss) 0127 % (17) remember probe snrs, models and results and recalculate the entire 0128 % array when changing precision 0129 % (18) could force probe SNRs to be a multiple of minimum step size 0130 % (19) use a non-uniform prior e.e. 1/(1+x^2) 0131 % (20) possibly use different parametrization (e.g. 1/slope or a second SRT threshold) 0132 % (22) save inputs so that pdfs can be recalculated when rescaling 0133 % (24) Parametrize slope as x=(100s^2-1)/20s to make the distribution more uniform 0134 % inverse is s=(sqrt(x^2+1)-x)/10; alternatively use log(slope) 0135 % (25) optionally have no prior (i.e. maximum likelihood) 0136 % (26) Check why the estimated variances are too small 0137 % (27) Adapt range of potential probes if optimum is at the limit 0138 % (28) Resample the pdf on the basis of the marginal distributions 0139 % (29) Select the probe range on the basis of marginal distributions 0140 % (30) Resample the pdfs by a factor of 2 always 0141 % (31) use uneven samples in pdf concentrated in the middle 0142 % (32) Use a parametric mixture including one-sided constants for the pdf 0143 % (33) Selection of subset of prescribed probe SNRs is not optimum 0144 % (34) use quadratic interpolation to select the probe SNR unless using fixed values 0145 % (35) expand the pdf grid based on the effective number of samples (like particle filters) 0146 0147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0148 % persistent variables 0149 % sizes: nx=SNR values in pdf 0150 % ns=slope values in pdf (1/std dev of cumulative gaussian) 0151 % nh=potential test SNRs 0152 % ni=simultaneous models being estimated 0153 % wq(ns*nx,ni) = prob of model; each column is a vectorized ns*nx matrix 0154 % nr(1,10) = parameter values: [SNR-pdf slope-pdf SNR-probe ...] 0155 % pr(7,ni) = input model parameters 0156 % qr(4,ni) = derived model parameters 0157 % xq(nr(1),ni) = SNR values in pdf 0158 % sq(nr(2),ni) = slope values in pdf 0159 % mq(2,ni) = estimated srt and slope of all models 0160 % vq(3,ni) = estimated covariance matrix entries:[var(srt) cov(srt,slope) var(slope)]' 0161 % xn(1,ni) = next probe value to use 0162 % hn(1,ni) = expected decrease in cost function after next probe 0163 0164 persistent wq xq sq nr pr qr mq vq xn hn hfact xz 0165 0166 % algorithm parametes that could be programmable 0167 pfloor=0.001; % pdf floor relative to uniform distribution 0168 trynsig=2; % number of standard deviations to explore 0169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0170 if iq<0 % initialization 0171 ni=-iq; % number of models 0172 pr=repmat([0.75 0.04 0.1 -20 20 0 0.5]',1,ni); % default parameters 0173 if nargin>1 0174 if size(x,2)>1 0175 if size(x,2)>ni 0176 error('initialization parameter argument has too many columns'); 0177 end 0178 pr(1:size(x,1),1:size(x,2))=x; 0179 else 0180 pr(1:size(x,1),:)=repmat(x,1,ni); 0181 end 0182 end 0183 nr=[40 21 30 1 0.2 0.02 4 1.3 2 1 1]'; % default parameter values 0184 nrf={'nx','ns','nh','cs','dh','sl','kp','hg','cf','pm','lg'}; % parameter field names 0185 numnr=length(nr); 0186 if nargin>2 0187 if isstruct(r) 0188 fnn=fieldnames(r); 0189 for i=1:length(fnn) 0190 mk=strcmp(fnn{i},nrf); 0191 if any(mk) 0192 nr(mk)=r.(fnn{i}); 0193 end 0194 end 0195 else 0196 nr(1:min(numnr,numel(r)))=r(:); 0197 end 0198 nr(1:3)=round(nr(1:3)); 0199 end 0200 pr(6,:)=max(pr(6,:),nr(6)); % low limit of slope in prob/dB 0201 nsxq=nr(1)*nr(2); 0202 xq=(0:nr(1)-1)'*(pr(5,:)-pr(4,:))/(nr(1)-1)+repmat(pr(4,:),nr(1),1); 0203 if nr(11) % if log slope 0204 sq=(1-nr(2):0)'*(log(pr(7,:))-log(max(pr(6,:),nr(6))))/(nr(2)-1)+repmat(log(pr(7,:)),nr(2),1); 0205 else 0206 sq=(0:nr(2)-1)'*(pr(7,:)-pr(6,:))/(nr(2)-1)+repmat(pr(6,:),nr(2),1); 0207 end 0208 wq=repmat(1/nsxq,nsxq,ni); % initialize to uniform pdf 0209 qr=zeros(5,ni); 0210 qr(1,:)=1-pr(2,:)-pr(3,:); % prob range covered by cumulative gaussian 0211 qr(2,:)=(pr(1,:)-pr(3,:))./qr(1,:); % cumulative gaussian probability at threshold 0212 switch nr(10) 0213 case 1 0214 qr(3,:)=log(qr(2,:)./(1-qr(2,:))); % x position of target in std measure 0215 qr(4,:)=qr(1,:).*qr(2,:).*(1-qr(2,:)); % slope*"stddev" at threshold 0216 case 2 0217 qr(3,:)=norminv(qr(2,:)); % x position of target in std measure 0218 qr(4,:)=qr(1,:).*normpdf(qr(3,:)); % slope*stddev at threshold 0219 otherwise 0220 error('Unrecognised psychometric model selection'); 0221 end 0222 mq=repmat([mean(xq,1); mean(sq,1)],[1,1,3]); % initial means 0223 vq=[var(xq,1,1); zeros(1,ni); var(sq,1,1)]; % initial variances 0224 % hn=[1-nr(4) 0 nr(4)]*vq; % artificially high value of cost function ensures all models are probed early on 0225 hn=repmat(Inf,1,ni); % very high initial cost function 0226 hfact=nr(8)^(1/ni); % growth factor to ensure that no model is neglected for too long 0227 xn=mq(1,:); % start at mean value 0228 if nargin>=4 && ~isempty(xp) 0229 if iscell(xp) 0230 xz=xp; 0231 else 0232 xz=repmat(num2cell(xp(:)',2),ni,1); % else replicate for each model 0233 end 0234 for i=1:ni 0235 [dummy,j]=min(abs(xz{i}-mq(1,i))); % select the closest available probe to the mean 0236 xn(i)=xz{i}(j); 0237 end 0238 else 0239 xz=cell(ni,1); % make empty cells 0240 end 0241 elseif iq>0 && nargin==3 % update pdfs with a new probe result 0242 nxq=nr(1); 0243 nsq=nr(2); 0244 nxh=nr(3); 0245 nsxq=nxq*nsq; 0246 % thresh=pr(1,iq); % target threshold 0247 guess=pr(3,iq); % guess rate (1/choices) 0248 pscale=qr(1,iq); % prob range left after subtracting miss and guess probs 0249 xtstd=qr(3,iq); % x position of target in std measure 0250 sqi=sq(:,iq); % slope (or log slope) values in PDF 0251 if nr(11) % if log slope 0252 sqis=exp(sqi)/qr(4,iq); % inverse std dev of gaussian (proportional to slope) 0253 else 0254 sqis=sqi/qr(4,iq); % inverse std dev of gaussian (proportional to slope) 0255 end 0256 xqi=xq(:,iq); 0257 wqi=wq(:,iq); 0258 % 0259 % update probabilities with the previous test result 0260 % If r==1, we multiply by a horizontally reflected version of the psychometric function that equals 0261 % the threshold (e.g. 0.75) at the probe SNR, x, for all slopes. 0262 % If r==0, we multiply by 1- this reflected version which therefore equals (1-thresh) at x. 0263 % 0264 r0=r==0; 0265 switch nr(10) 0266 case 1 0267 wqi=wqi.*(r0+(1-2*r0)*(guess+pscale*((1+exp(reshape(sqis*xqi'-xtstd,nsxq,1)-repmat(sqis,nxq,1)*x)).^(-1)))); % P(l | r,x) 0268 case 2 0269 wqi=wqi.*(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,nxq,1)*x-reshape(sqis*xqi'-xtstd,nsxq,1)))); % P(l | r,x) 0270 otherwise 0271 error('Unrecognised psychometric model selection'); 0272 end 0273 wqi=wqi./sum(wqi,1); % normalize 0274 % wq(:,iq)=wqi; % save updated probabilities (not needed if we always interpolate them) 0275 0276 % **** this section copied to graph plotting **** 0277 % Calculate mean and covariance and entropy 0278 0279 wqsx=reshape(wqi,nsq,nxq); 0280 px=sum(wqsx,1); % p(x0) 0281 ps=sum(wqsx,2); % p(s0) 0282 xe=px*xqi; % E(x0) 0283 se=ps'*sqi; % E(s0) 0284 0285 [pxpk,xm]=max(px); % marginal mode in x 0286 if xm>1 && xm<nxq % use quadratic interpolation if possible 0287 [dum,xm2]=quadpeak(px(xm-1:xm+1)'); 0288 xm=xm+xm2-2; 0289 end 0290 xm=(2-xm)*xqi(1)+(xm-1)*xqi(2); % marginal mode(x0) 0291 0292 [pspk,sm]=max(ps); 0293 if sm>1 && sm<nsq % use quadratic interpolation if possible 0294 [dum,sm2]=quadpeak(ps(sm-1:sm+1)); 0295 sm=sm+sm2-2; 0296 end 0297 sm=(2-sm)*sqi(1)+(sm-1)*sqi(2); % marginal mode(s0) 0298 0299 [wqpk,j]=max(wqi); 0300 i=1+floor((j-1)/nsq); 0301 j=j-nsq*(i-1); 0302 if i>1 && i<nxq && j>1 && j<nsq % use quadratic interpolation if possible 0303 [dum,ji]=quadpeak(wqsx(j-1:j+1,i-1:i+1)); 0304 i=i+ji(2)-2; 0305 j=j+ji(1)-2; 0306 end 0307 xj=(2-i)*xqi(1)+(i-1)*xqi(2); % joint mode x 0308 sj=(2-j)*sqi(1)+(j-1)*sqi(2); % joint mode: s 0309 xv=px*(xqi.^2)-xe^2; % Var(x0) 0310 sv=ps'*(sqi.^2)-se^2; % Var(s0) 0311 sxv=wqi'*(repmat(sqi,nxq,1).*reshape(repmat(xqi',nsq,1),nsxq,1))-xe*se; % Cov(s0*x0) 0312 0313 % ************ end of copied section 0314 0315 mq(:,iq,1)=[xe; se]; % save means 0316 mq(:,iq,2)=[xj; sj]; % save means 0317 mq(:,iq,3)=[xm; sm]; % save means 0318 vq(:,iq)=[xv; sxv; sv]; % save covariance matrix 0319 xh=(px*log(px)')*(xqi(1)-xqi(2)); % h(x0) 0320 sh=(ps'*log(ps))*(sqi(1)-sqi(2)); % h(s0) 0321 0322 % now estimate the next probe SNR $$$$ 0323 0324 if ~numel(xz{iq}) % if no list of probe SNRs was specified 0325 ytry=exp(0.25i*pi*(0:7))'; % points around the circle 0326 ytry=[real(ytry) imag(ytry)]; 0327 [vtry,dtry]=eig([xv sxv; sxv sv]); % eigendecomposition of covariance matrix 0328 tryxs=repmat([xe,se],8,1)+trynsig*ytry*sqrt(dtry)*vtry'; 0329 pmin=0.05; % target probe success probability 0330 if nr(11) 0331 tryxs(:,2)=qr(4,iq)*exp(-tryxs(:,2)); % convert log(slope) to std dev 0332 else 0333 tryxs(:,2)=qr(4,iq)*tryxs(:,2).^(-1); % convert slope to std dev 0334 end 0335 switch nr(10) 0336 case 1 % logistic 0337 qmax=max(tryxs(:,1)+(log((1-pmin)/pmin)-xtstd)*tryxs(:,2)); 0338 qmin=min(tryxs(:,1)+(log(pmin/(1-pmin))-xtstd)*tryxs(:,2)); 0339 case 2 % cumulative gaussian 0340 qmax=max(tryxs(:,1)+(norminv(1-pmin)-xtstd)*tryxs(:,2)); 0341 qmin=min(tryxs(:,1)+(norminv(pmin)-xtstd)*tryxs(:,2)); 0342 end 0343 dxt=max(nr(5),(qmax-qmin)/nxh); % minimum step size of nr(5) [0.2 dB] 0344 xt=(qmin+qmax)/2+((1:nxh)-(1+nxh)/2)*dxt; 0345 else % if a specific list of probe SNRs exists 0346 xzi=xz{iq}; % xzi is the list of available probe SNRs 0347 if numel(xzi)<=nxh % use all available probe SNRs if there are not too many 0348 xt=xzi; 0349 else 0350 [xt,ixt]=min(abs(xzi-xe)); % find the closest one to xe ** not necessarily optimum *** 0351 ixt=max(1,min(1+numel(xzi)-nxh,ixt-floor((1+nxh)/2))); % arrange symmetrically around xt 0352 xt=xzi(ixt:min(ixt+nxh-1,numel(xzi))); 0353 end 0354 end 0355 nxhp=length(xt); % xt are the potential probe SNRs 0356 switch nr(10) 0357 case 1 0358 prt=guess+pscale*((1+exp(repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)-repmat(sqis,nxq,1)*xt)).^(-1)); % P(r | l,x) 0359 case 2 0360 prt=guess+pscale*normcdf(repmat(sqis,nxq,1)*xt-repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)); % P(r | l,x) 0361 end 0362 wqt=repmat(wqi,1,nxhp); 0363 pl1=prt.*wqt; % posterior prob given success = p(l | x,r=1) unnormalized 0364 pl0=(wqt-pl1); % posterior prob given failure = p(l | x,r=0) unnormalized 0365 prh=sum(pl1,1); % p(r | x)=Sum{P(r | l,x)*P(l)} [note wqt is normalized] (row vector) 0366 pl1=pl1./repmat(prh,nsxq,1); % normalize 0367 pl0=pl0./repmat(1-prh,nsxq,1); % posterior prob given failure = p(l | x,r=0) normalized 0368 px1=squeeze(sum(reshape(pl1,nsq,nxq,[]),1)); % p(x0 | x,r=1) 0369 px0=squeeze(sum(reshape(pl0,nsq,nxq,[]),1)); % p(x0 | x,r=0) 0370 ps1=squeeze(sum(reshape(pl1,nsq,nxq,[]),2)); % p(s0 | x,r=1) 0371 ps0=squeeze(sum(reshape(pl0,nsq,nxq,[]),2)); % p(s0 | x,r=0) 0372 xet1=xqi'*px1; % E(x0 | x,r=1) 0373 xvt1=(xqi.^2)'*px1-xet1.^2; % Var(x0 | x,r=1) 0374 xet0=xqi'*px0; % E(x0 | x,r=0) 0375 xvt0=(xqi.^2)'*px0-xet0.^2; % Var(x0 | x,r=0) 0376 xvt=xvt1.*prh+xvt0.*(1-prh); % E(Var(x0 | x )) 0377 set1=sqi'*ps1; % E(s0 | x,r=1) 0378 svt1=(sqi.^2)'*ps1-set1.^2; % Var(s0 | x,r=1) 0379 set0=sqi'*ps0; % E(s0 | x,r=0) 0380 svt0=(sqi.^2)'*ps0-set0.^2; % Var(s0 | x,r=0) 0381 svt=svt1.*prh+svt0.*(1-prh); % E(Var(s0 | x )) 0382 xht1=sum(log(px1).*px1,1); % -H(x0 | x,r=1) 0383 xht0=sum(log(px0).*px0,1); % -H(x0 | x,r=0) 0384 xht=(xht1.*prh+xht0.*(1-prh))*(xqi(1)-xqi(2)); % h(x0 | x) 0385 sht1=sum(log(ps1).*ps1,1); % -H(s0 | x,r=1) 0386 sht0=sum(log(ps0).*ps0,1); % -H(s0 | x,r=0) 0387 sht=(sht1.*prh+sht0.*(1-prh))*(sqi(1)-sqi(2)); % h(s0 | x) 0388 switch nr(9) 0389 case 1 0390 hx=(xvt + nr(4)*svt)/(1+nr(4)); % cost function for each possible test SNR 0391 [hxmin,ix]=min(hx); % find the minimum of cost function 0392 hn(iq)=(xv + nr(4)*sv)/(1+nr(4))-hxmin; % expected decrease in cost function 0393 case 2 0394 hx=(xht + nr(4)*sht)/(1+nr(4)); % cost function for each possible test SNR 0395 [hxmin,ix]=min(hx); % find the minimum of cost function 0396 hn(iq)=(xh + nr(4)*sh)/(1+nr(4))-hxmin; % expected decrease in cost function 0397 otherwise 0398 error('Unrecognised cost function option'); 0399 end 0400 xn(iq)=xt(ix); % next probe value for this model 0401 % fprintf('Probe range: %.3g %.3g; choose %.3g\n',xt(1),xt(end),xt(ix)); 0402 0403 0404 % rescale the pdfs if necessary 0405 0406 ssd=sqrt(sv); 0407 xsd=sqrt(xv); 0408 if nr(11) 0409 sq2=linspace(max(log(nr(6)),se-nr(7)*ssd),se+nr(7)*ssd,nsq)'; 0410 else 0411 sq2=linspace(max(nr(6),se-nr(7)*ssd),se+nr(7)*ssd,nsq)'; 0412 end 0413 xq2=linspace(xe-nr(7)*xsd,xe+nr(7)*xsd,nxq)'; % new x axis values 0414 % do linear interpolation in the x direction 0415 wqi=reshape(wqi,nsq,nxq); % turn into a matrix for easy interpolation 0416 xqf=(xq2-xq(1,iq))/(xq(2,iq)-xq(1,iq)); 0417 xqj=ceil(xqf); 0418 xqf=xqj-xqf; 0419 xqg=1-xqf; 0420 xqf((xqj<=0) | (xqj>nxq))=0; 0421 xqg((xqj<0) | (xqj>=nxq))=0; 0422 wq2=wqi(:,min(max(xqj,1),nxq)).*repmat(xqf',nsq,1)+wqi(:,min(max(xqj+1,1),nxq)).*repmat(xqg',nsq,1); 0423 % do linear interpolation in the s direction 0424 sqf=(sq2-sq(1,iq))/(sq(2,iq)-sq(1,iq)); 0425 sqj=ceil(sqf); 0426 sqf=sqj-sqf; 0427 sqg=1-sqf; 0428 sqf((sqj<=0) | (sqj>nsq))=0; 0429 sqg((sqj<0) | (sqj>=nsq))=0; 0430 wq2=wq2(min(max(sqj,1),nsq),:).*repmat(sqf,1,nxq)+wq2(min(max(sqj+1,1),nsq),:).*repmat(sqg,1,nxq); 0431 % now normalize and apply a floor 0432 wq2=wq2(:); % turn back into a vector 0433 wq2=wq2/sum(wq2,1); % normalize 0434 wq2=max(wq2,pfloor/nsxq); %impose a floor 0435 wq2=wq2/sum(wq2,1); % normalize again 0436 sq(:,iq)=sq2; 0437 xq(:,iq)=xq2; 0438 wq(:,iq)=wq2; 0439 elseif iq==0 % output model parameters 0440 xx=pr; 0441 ii=nr; 0442 if ~nargout 0443 pdesc={'Threshold','Miss prob','Guess Prob','Min SNR','Max SNR','Min Slope','Max Slope'}; 0444 fprintf('\n*********\nModel-specific Parameters\n'); 0445 for i=1:7 0446 fprintf('%4d) %s: ',i,pdesc{i}); 0447 if size(pr,2)>1 0448 fprintf('%.5g, ',pr(i,1:size(pr,2)-1)); 0449 end 0450 fprintf('%.5g\n',pr(i,size(pr,2))); 0451 end 0452 qdesc={'nx SNR values in PDF', ... 0453 'ns Slope values in PDF', ... 0454 'nh Probe SNR values to evaluate', ... 0455 'cs Weighting of slope relative to threshold in cost function', ... 0456 'dh Min step size in dB for probe SNRs', ... 0457 'sl Min slope at threshold', ... 0458 'kp Std deviations of the pdfs to keep', ... 0459 'hg Amount to grow expected gains in ni trials', ... 0460 'cf Cost function', ... 0461 'pm Psychometric model'}; 0462 qoptf=[9 10]; % fields with options 0463 qoptval={'variance','entropy'; ... 0464 'logistic','cumulative gaussian'}; 0465 fprintf('\nShared Parameters\n'); 0466 for i=1:10 0467 fprintf('%4d) %s: ',i,qdesc{i}); 0468 j=find(qoptf==i,1); 0469 if numel(j) 0470 fprintf('%d=%s\n',nr(i),qoptval{j,nr(i)}); 0471 else 0472 fprintf('%.5g\n',nr(i)); 0473 end 0474 end 0475 end 0476 end 0477 0478 % now select the appropriate model to probe next 0479 0480 if iq~=0 0481 [hnmin,ii]=max(hn); % chose model with the biggest expected decrease 0482 hn=hn*hfact; % increase values to ensure they all get a chance 0483 xx=xn(ii); 0484 m=mq; 0485 v=vq; 0486 if nr(11) 0487 m(2,:,:)=exp(m(2,:,:)); % convert to real slope 0488 v(2,:)=v(2,:).*m(2,:,1); % correct the covariance 0489 v(3,:)=v(3,:).*m(2,:,1).^2; % and the slope variance 0490 end 0491 end 0492 0493 if ~nargout && iq>0 0494 sqi=sq(:,iq); 0495 nsq=length(sqi); 0496 xqi=xq(:,iq); 0497 nxq=length(xqi); 0498 wqi=wq(:,iq); 0499 nsxq=nxq*nsq; 0500 0501 % **** this section copied from main processing **** 0502 % Calculate mean and covariance and entropy 0503 0504 wqsx=reshape(wqi,nsq,nxq); 0505 px=sum(wqsx,1); % p(x0) 0506 ps=sum(wqsx,2); % p(s0) 0507 xe=px*xqi; % E(x0) 0508 se=ps'*sqi; % E(s0) 0509 [pxpk,xm]=max(px); 0510 if xm>1 && xm<nxq % use quadratic interpolation if possible 0511 [dum,xm2]=quadpeak(px(xm-1:xm+1)'); 0512 xm=xm+xm2-2; 0513 end 0514 xm=(2-xm)*xqi(1)+(xm-1)*xqi(2); % marginal mode(x0) 0515 [pspk,sm]=max(ps); 0516 if sm>1 && sm<nsq % use quadratic interpolation if possible 0517 [dum,sm2]=quadpeak(ps(sm-1:sm+1)); 0518 sm=sm+sm2-2; 0519 end 0520 sm=(2-sm)*sqi(1)+(sm-1)*sqi(2); % marginal mode(s0) 0521 [wqpk,j]=max(wqi); 0522 i=1+floor((j-1)/nsq); 0523 j=j-nsq*(i-1); 0524 if i>1 && i<nxq && j>1 && j<nsq % use quadratic interpolation if possible 0525 [dum,ji]=quadpeak(wqsx(j-1:j+1,i-1:i+1)); 0526 i=i+ji(2)-2; 0527 j=j+ji(1)-2; 0528 end 0529 xj=(2-i)*xqi(1)+(i-1)*xqi(2); % joint mode x 0530 sj=(2-j)*sqi(1)+(j-1)*sqi(2); % joint mode: s 0531 xv=px*(xqi.^2)-xe^2; % Var(x0) 0532 sv=ps'*(sqi.^2)-se^2; % Var(s0) 0533 sxv=wqi'*(repmat(sqi,nxq,1).*reshape(repmat(xqi',nsq,1),nsxq,1))-xe*se; % Cov(s0*x0) 0534 0535 % ************ end of copied section 0536 0537 % draw final 2-D distribution 0538 nxq=nr(1); 0539 nsq=nr(2); 0540 axm=[3 -2; 2 -1]; 0541 xqi2=[axm*xqi(1:2);xqi]; 0542 sqi2=[axm*sqi(1:2);sqi]; 0543 imagesc(xqi2,sqi2,[zeros(1,2) px/pxpk; zeros(1,nxq+2);ps/pspk zeros(nsq,1) wqsx/wqpk]); 0544 hold on 0545 plot(xe,se,'wo',xm,sm,'w^',xj,sj,'w*'); 0546 plot(xqi2(2),se,'wo',xqi2(2),sm,'w^',xqi2(2),sj,'w*'); 0547 plot(xe,sqi2(2),'wo',xm,sqi2(2),'w^',xj,sqi2(2),'w*'); 0548 t=linspace(0,2*pi,200); 0549 xcir=cos(t); 0550 scir=sin(t); 0551 vcir=sqrt((sv*xcir.^2+xv*scir.^2-2*sxv*xcir.*scir)/(xv*sv-sxv^2)); 0552 plot(xe+xcir./vcir,se+scir./vcir,'w-'); 0553 hold off 0554 axis 'xy'; 0555 % colorbar; 0556 % cblabel('Relative Probability'); 0557 xlabel(sprintf('SNR @ %d%%SRT (dB)',round(pr(1)*100))); 0558 if nr(11) 0559 ylabel('Log psychometric Slope at threshold (prob/dB)'); 0560 else 0561 ylabel('Psychometric Slope at threshold (prob/dB)'); 0562 end 0563 title('Joint pdf: o mean, * mode, \Delta marginal mode'); 0564 end