V_PSYCEST estimate multiple psychometric functions
Usage: [xx,ii,m,v]=v_psycest(-n,p,q,xp,lf) % initialize n models
[xx,ii,m,v]=v_psycest(i,x,r) % supply a trial result to v_psycest
[xx,ii,m,v]=v_psycest(i,x,r,o) % supply a trial result to v_psycest and plot
v_psycest(i,o) % plot pdf of model i with plot options o
[xx,ii,m,v,mr,vr]=v_psycest(i,o) % Also determine robust outputs which exclude outliers [takes longer]
[p,q,msr]=psycest(0) % output model parameters (or print them if no outputs)
Inputs (see usage examples for argument positions):
-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 (but over-ridden if <q.sl) [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 SRT in cost function [1]
5 q.dh minimum step size in dB for probe SNRs [0.2]
6 q.sl min slope at threshold (must be >0) [0.005]
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=v_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]
12 q.pp Number of prior standard deviations in initial semi-range [1]
13 q.pf Probability floor (integrated over entire grid) [0.0001]
14 q.ts Number of std devs to explore [2]
15 q.dp Maximum probe SNR shift (dB) [10]
16 q.it Grid interpolation threshold [0.5]
17 q.at Axis change threshold [0.1]
18 q.la Look 2-ahead when choosing probe [1]
19 q.op Outlier probability [0.01]
20 q.rx Minimum range factor per iteration [0.5]
xp{n}(:) list of available probe SNR values
lf log file ID
i model probed
x probe SNR value used
r response: 0=wrong, 1=right.
o plot option string:
'p' Plot pdf
'h' Plot probe history
'x' Plot expected cost function for probe SNRs
'c' Plot cost function evolution
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)]'
mr(2,n,3)robust estimated srt and slope of all models (exclusing outliers)
m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
vr(3,n) robust estimated covariance matrix entries (exclusing outliers):
[var(srt) cov(srt,slope) var(slope)]'
msr(:,7) List probe snrs and results: [model probe-snr result xe se xv sv]
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 free 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 size of the joint pdf in SNR and Slope respectively.
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 v_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, v_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 v_entropy of the estimated SNR and slope distributions.
My experience is that v_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" (0) or "log(slope)" (1). It is recommended
that you stick to the default of log(slope) since this results in more symmetrical distributions. [1]
q(12) This sets the number of std deviations of the prior corresponding to the limits of the
initial distribution set in p(4:7,n). A small value corresponds to a very weak prior. [1]
q(13) Probability floor (integrated over entire grid) [0.0001]
q(14) Number of std devs to explore when evaluating possible probe SNRs [2]
q(15) Maximum probe SNR shift beyond the range of previous probes (dB) [10]
q(16) Grid interpolation threshold. The grid is only interpolated if
it has shrunk by less that this factor [0.5]
q(17) If the desired grid limits change by less than this fraction of
the grid length, they will be left unaltered. [0.1]
q(18) Look 2-ahead when choosing probe [1]
q(19) Outlier probability; probe results with lower probability than
this will be ignored [0.01]
q(20) The ranges of SRT and log slope will not be reduced per iteration below this factor [0.5]0001 function [xx,ii,m,v,mr,vr]=v_psycest(iq,x,r,xp,lf) 0002 %V_PSYCEST estimate multiple psychometric functions 0003 % 0004 % Usage: [xx,ii,m,v]=v_psycest(-n,p,q,xp,lf) % initialize n models 0005 % [xx,ii,m,v]=v_psycest(i,x,r) % supply a trial result to v_psycest 0006 % [xx,ii,m,v]=v_psycest(i,x,r,o) % supply a trial result to v_psycest and plot 0007 % v_psycest(i,o) % plot pdf of model i with plot options o 0008 % [xx,ii,m,v,mr,vr]=v_psycest(i,o) % Also determine robust outputs which exclude outliers [takes longer] 0009 % [p,q,msr]=psycest(0) % output model parameters (or print them if no outputs) 0010 % 0011 % Inputs (see usage examples for argument positions): 0012 % -n minus the number of models 0013 % p(:,n) parameters for each model 0014 % 1 thresh [0.75] 0015 % 2 miss [0.04] 0016 % 3 guess [0.1] 0017 % 4 SNR min [-20] 0018 % 5 SNR max [20] 0019 % 6 Slope min (but over-ridden if <q.sl) [0] 0020 % 7 slope max [0.5] 0021 % q(:) parameters common to all models (vector or struct) 0022 % 1 q.nx number of SNR values in pdf [40] 0023 % 2 q.ns number of slope values in pdf [21] 0024 % 3 q.nh number of probe SNR values to evaluate [30] 0025 % 4 q.cs weighting of slope relative to SRT in cost function [1] 0026 % 5 q.dh minimum step size in dB for probe SNRs [0.2] 0027 % 6 q.sl min slope at threshold (must be >0) [0.005] 0028 % 7 q.kp number of std deviations of the pdfs to keep [4] 0029 % 8 q.hg amount to grow expected gains in ni trials [1.3] 0030 % 9 q.cf cost function: 1=variance, 2=v_entropy [2] 0031 % 10 q.pm psychometric model: 1=logistic, 2=cumulative gaussian [1] 0032 % 11 q.lg use log slope in pdf: 0=no, 1=yes [1] 0033 % 12 q.pp Number of prior standard deviations in initial semi-range [1] 0034 % 13 q.pf Probability floor (integrated over entire grid) [0.0001] 0035 % 14 q.ts Number of std devs to explore [2] 0036 % 15 q.dp Maximum probe SNR shift (dB) [10] 0037 % 16 q.it Grid interpolation threshold [0.5] 0038 % 17 q.at Axis change threshold [0.1] 0039 % 18 q.la Look 2-ahead when choosing probe [1] 0040 % 19 q.op Outlier probability [0.01] 0041 % 20 q.rx Minimum range factor per iteration [0.5] 0042 % xp{n}(:) list of available probe SNR values 0043 % lf log file ID 0044 % i model probed 0045 % x probe SNR value used 0046 % r response: 0=wrong, 1=right. 0047 % o plot option string: 0048 % 'p' Plot pdf 0049 % 'h' Plot probe history 0050 % 'x' Plot expected cost function for probe SNRs 0051 % 'c' Plot cost function evolution 0052 % 0053 % Outputs: 0054 % xx recommended probe SNR 0055 % ii recommended model to probe next 0056 % m(2,n,3) estimated srt and slope of all models 0057 % m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates 0058 % v(3,n) estimated covariance matrix entries: 0059 % [var(srt) cov(srt,slope) var(slope)]' 0060 % mr(2,n,3)robust estimated srt and slope of all models (exclusing outliers) 0061 % m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates 0062 % vr(3,n) robust estimated covariance matrix entries (exclusing outliers): 0063 % [var(srt) cov(srt,slope) var(slope)]' 0064 % msr(:,7) List probe snrs and results: [model probe-snr result xe se xv sv] 0065 % 0066 % Algorithm parameters: 0067 % 0068 % The algorithm estimates the psychometric function for a number of models simultaneously. 0069 % A psychometric function specifies the probability of correct recognition as a function of 0070 % SNR and is completely specified by two free parameters: the SNR (in dB) and the slope (in 1/dB) 0071 % at a specified target recognition probability (e.g. 0.5 or 0.75). The p(:,n) parameters specify 0072 % some of the details of the psychometric function and can be different for each model: 0073 % p(1,n) gives the target recognition probability 0074 % p(2,n) gives the probability of incorrect recognition at very good SNRs (the "miss" probability). 0075 % If this value is made too low, a single unwarrented recognition error will have a large 0076 % effect of the estimated parameters and greatly lengthen the time to convergence. The default 0077 % value is 4%. 0078 % p(3,n) gives the probability of correct recognition at very poor SNRs (the "guess" probabiity). 0079 % This should be set to 1/N where N is the number of possible responses to a stimulus. 0080 % p(4:5,n) These give the initial min and max SNR at the target recognition probability. They will 0081 % be adjusted by the algorithm if necessary but wrong values will delay convergence. 0082 % p(6:7,n) These given the initial min and max slope (in probability per dB) at the target 0083 % recognition probability. 0084 % The remaining parameters are shared between all the models and control how the algorithm operates. 0085 % q(1:2) These parameters determine the sampling size of the joint pdf in SNR and Slope respectively. 0086 % q(3) This parameter specifies how many potential values of probe SNR to evaluate in order to 0087 % determine which is likely to give the most improvement in the parameter estimates. 0088 % q(4) This parameter gives the relative weight of SNR and Slope in the cost function. Increasing 0089 % its value will improve the slope estimate (or log(slope) estimate) at the expense of the 0090 % SNR estimate. Actually its value is not all that critical. 0091 % q(5) At each iteration v_psycest evaluates several potential probe SNR levels to see which is 0092 % likely to give the most improvement to the parameter estimates. This parameter fixes the 0093 % minimum spacing between these potential probe values. It only has an effect when the variances 0094 % of the parameter estimates are very small. 0095 % q(6) To prevent the routine getting lost, this parameter specifies the smallest reasonable value 0096 % of the Slope parameter at the target recognition probability. Under normal circumstances, it 0097 % has no effect. 0098 % q(7) For each model, the routine maintains a sampled version of the joint pdf of the SNR and slope. 0099 % The sampling grid is reclculated after each iteration to encompass this number of standard 0100 % deviations in each dimension. 0101 % q(8) At each iteration, v_psycest advises which model to probe next on the basis of which 0102 % gives the greatest expected reduction in cost function. To ensure that all models 0103 % are periodically sampled, this expected reduction of each unprobed model is multiplied 0104 % by this parameter at each iteration. Thus a model will eventually be probed even if 0105 % its expected cost factor improvement is small. 0106 % q(9) This determines whether the recommended probe SNR to use at the next iteration is chosen to 0107 % minimize the expected variance or the v_entropy of the estimated SNR and slope distributions. 0108 % My experience is that v_entropy (the default) gives faster convergence. 0109 % q(10) This selects whether the underlying model is a logistic function (1) or a cumulative 0110 % gaussian (2). The choice makes little difference unless the "miss" or "guess" probabilities 0111 % are very very small. 0112 % q(11) This selects whether the internal pdf samples "slope" (0) or "log(slope)" (1). It is recommended 0113 % that you stick to the default of log(slope) since this results in more symmetrical distributions. [1] 0114 % q(12) This sets the number of std deviations of the prior corresponding to the limits of the 0115 % initial distribution set in p(4:7,n). A small value corresponds to a very weak prior. [1] 0116 % q(13) Probability floor (integrated over entire grid) [0.0001] 0117 % q(14) Number of std devs to explore when evaluating possible probe SNRs [2] 0118 % q(15) Maximum probe SNR shift beyond the range of previous probes (dB) [10] 0119 % q(16) Grid interpolation threshold. The grid is only interpolated if 0120 % it has shrunk by less that this factor [0.5] 0121 % q(17) If the desired grid limits change by less than this fraction of 0122 % the grid length, they will be left unaltered. [0.1] 0123 % q(18) Look 2-ahead when choosing probe [1] 0124 % q(19) Outlier probability; probe results with lower probability than 0125 % this will be ignored [0.01] 0126 % q(20) The ranges of SRT and log slope will not be reduced per iteration below this factor [0.5] 0127 0128 % Copyright (C) Mike Brookes and Clement Doire 2009-2016 0129 % Version: $Id: v_psycest.m 10865 2018-09-21 17:22:45Z dmb $ 0130 % 0131 % VOICEBOX is a MATLAB toolbox for speech processing. 0132 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0133 % 0134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0135 % This program is free software; you can redistribute it and/or modify 0136 % it under the terms of the GNU General Public License as published by 0137 % the Free Software Foundation; either version 2 of the License, or 0138 % (at your option) any later version. 0139 % 0140 % This program is distributed in the hope that it will be useful, 0141 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0142 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0143 % GNU General Public License for more details. 0144 % 0145 % You can obtain a copy of the GNU General Public License from 0146 % http://www.gnu.org/copyleft/gpl.html or by writing to 0147 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0149 0150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0151 % Bugs/Suggestions: 0152 % (2)use a structure for input parameters 0153 % (3) v_entropy seems to change signifiacntly when the grid is rescaled 0154 % (4)add a forgetting factor for old measurements 0155 % (8) use quadratic interpolation to find cost function minimum 0156 % (10) optionally output the whole pdf + its axis values 0157 % (11) optionally output all model probe values 0158 % (13) Should perhaps estimate and return the mean and slope compensated for the guess rate 0159 % i.e. if you want p, you actually test at guess+ p*(1-guess-miss)/(1-miss) 0160 % (17) remember probe snrs, models and results and recalculate the entire 0161 % array when changing precision 0162 % (18) could force probe SNRs to be a multiple of minimum step size 0163 % (19) use a non-uniform prior e.e. 1/(1+x^2) 0164 % (20) possibly use different parametrization (e.g. 1/slope or a second SRT threshold) 0165 % (22) save inputs so that pdfs can be recalculated when rescaling 0166 % (24) Parametrize slope as x=(100s^2-1)/20s to make the distribution more uniform 0167 % inverse is s=(sqrt(x^2+1)-x)/10; alternatively use log(slope) 0168 % (25) optionally have no prior (i.e. maximum likelihood) 0169 % (26) Check why the estimated variances are too small 0170 % (27) Adapt range of potential probes if optimum is at the limit 0171 % (28) Resample the pdf on the basis of the marginal distributions 0172 % (29) Select the probe range on the basis of marginal distributions 0173 % (30) Resample the pdfs by a factor of 2 always 0174 % (31) use uneven samples in pdf concentrated in the middle 0175 % (32) Use a parametric mixture including one-sided constants for the pdf 0176 % (33) Selection of subset of prescribed probe SNRs is not optimum 0177 % (34) use quadratic interpolation to select the probe SNR unless using fixed values 0178 % (35) expand the pdf grid based on the effective number of samples (like particle filters) 0179 0180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0181 % persistent variables 0182 % sizes: nx=SNR values in pdf 0183 % ns=slope values in pdf (1/std dev of cumulative gaussian) 0184 % nh=potential test SNRs 0185 % ni=simultaneous models being estimated 0186 % wq(ns*nx,ni) = log prob of model; each column is a vectorized ns*nx matrix 0187 % nr(1,10) = parameter values: [SNR-pdf slope-pdf SNR-probe ...] 0188 % pr(7,ni) = input model parameters 0189 % qr(4,ni) = derived model parameters 0190 % xq(nr(1),ni) = SNR values in pdf 0191 % sq(nr(2),ni) = slope values in pdf 0192 % mq(2,ni,3) = estimated srt and slope of all models 0193 % vq(3,ni) = estimated covariance matrix entries:[var(srt) cov(srt,slope) var(slope)]' 0194 % xn(1,ni) = next probe value to use 0195 % hn(1,ni) = expected decrease in cost function after next probe 0196 % hfact 0197 % xz 0198 % res(nres,7) = results [model-number, probe-SNR, result] 0199 % nres = total number of probe results 0200 % nresq(1,ni) = number of probe results for each model 0201 % xmm(2,ni) = [min(probe-SNR); max(probe-SNR)] 0202 % mq0(2,ni) = prior means [x-mean; s-mean] 0203 % pq0(2,ni) = scale factors for prior distribution 0204 % wfl = floor relative to peak of log-pdf array 0205 % sqmin = minimum value of slope or log-slope 0206 % LOG = file ID of logging file 0207 % mqr(2,ni,3) = robust mean estimates (excluding outliers) 0208 % vqr(3,ni) = robust variance estimates (excluding outliers) 0209 % nresr = nres value at which robust estimates were last calculated 0210 0211 persistent wq xq sq nr pr qr mq vq xn hn hfact xz res nres nresq xmm mq0 pq0 wfl sqmin LOG mqr vqr nresr xlim 0212 0213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0214 % Initialization (iq<0) % 0215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0216 if iq<0 % initialization 0217 if nargin>4 && ~isempty(lf) 0218 LOG=lf; 0219 fprintf(LOG,'******************************\npsycest Initialization: %s\n******************************\n',datestr(now)); 0220 else 0221 LOG=0; 0222 end 0223 ni=-iq; % number of models 0224 nres=0; % total number of probe-results so far 0225 nresq=zeros(1,ni); % number of probe-results for each model 0226 res=zeros(1,7); % array for saving results [expands as needed] 0227 pr=repmat([0.75 0.04 0.1 -20 20 0 0.5]',1,ni); % default parameters 0228 if nargin>1 0229 if size(x,2)>1 0230 if size(x,2)>ni 0231 error('initialization parameter argument has too many columns'); 0232 end 0233 pr(1:size(x,1),1:size(x,2))=x; 0234 else 0235 pr(1:size(x,1),:)=repmat(x,1,ni); 0236 end 0237 end 0238 nr=[40 21 30 1 0.2 0.02 4 1.3 2 1 1 1 0.0001 2 10 0.5 0.1 1 0.01 0.5]'; % default parameter values 0239 nrf={'nx','ns','nh','cs','dh','sl','kp','hg','cf','pm','lg','pp', 'pf', 'ts', 'dp', 'it','at','la','op','rx'}; % parameter field names 0240 numnr=length(nr); 0241 if nargin>2 && ~isempty(r) % replace defaults with any parameter values specified in the call 0242 if isstruct(r) 0243 fnn=fieldnames(r); 0244 for i=1:length(fnn) 0245 mk=strcmp(fnn{i},nrf); 0246 if any(mk) 0247 nr(mk)=r.(fnn{i}); 0248 end 0249 end 0250 else 0251 nr(1:min(numnr,numel(r)))=r(:); % if parameters are specified as a vector, copy them across 0252 end 0253 nr(1:3)=round(nr(1:3)); % first three parameters must be integers 0254 end 0255 pr(6,:)=max(pr(6,:),nr(6)); % low limit of slope in prob/dB 0256 nxq=nr(1); 0257 nsq=nr(2); 0258 nsxq=nxq*nsq; 0259 xq=(0:nxq-1)'*(pr(5,:)-pr(4,:))/(nxq-1)+repmat(pr(4,:),nxq,1); 0260 if nr(11) % if log slope 0261 sqmin=log(nr(6)); % low limit for axis 0262 sq=(1-nsq:0)'*(log(pr(7,:))-log(max(pr(6,:),nr(6))))/(nsq-1)+repmat(log(pr(7,:)),nsq,1); 0263 else 0264 sqmin=nr(6); % low limit for axis 0265 sq=(1-nsq:0)'*(pr(7,:)-max(pr(6,:),nr(6)))/(nsq-1)+repmat(pr(7,:),nsq,1); 0266 end 0267 mq0=[mean(xq,1); mean(sq,1)]; % find mean of prior 0268 pq0=-2*nr(12)^2*[xq(end,:)-xq(1,:); sq(end,:)-sq(1,:)].^(-2); % scale factor for prior distribution 0269 wq=repmat(pq0(2,1)*(sq(:,1)-mq0(2,1)).^2,1,nxq)+repmat(pq0(1,1)*(xq(:,1)-mq0(1,1))'.^2,nsq,1); % log probs of prior (same for all models) 0270 wfl=log(nr(13)/nsxq); % floor relative to peak of wq array 0271 wq=repmat(wq(:),1,ni); % initialize to +-q.pp std deviations of gaussian prior 0272 qr=zeros(5,ni); 0273 qr(1,:)=1-pr(2,:)-pr(3,:); % prob range covered by cumulative gaussian 0274 if any(qr(1,:)<=0) 0275 [dum,i]=min(qr(1,:)); 0276 error('Model %d: guess prob (%.2f) + miss prob (%.2f) is >=1',i,pr(3,i),pr(2,:)); 0277 end 0278 qr(2,:)=(pr(1,:)-pr(3,:))./qr(1,:); % cumulative gaussian probability at threshold 0279 if any(abs(qr(2,:)-0.5)>=0.5) 0280 [dum,i]=max(qr(2,:)-0.5); 0281 error('Model %d: target SRT threshold (%.2f) must lie in the range guess prob (%.2f) to (1-miss prob) (%.2f)',i,pr(1,i),pr(3,i),1-pr(2,:)); 0282 end 0283 switch nr(10) % switch on model type 0284 case 1 % logistic model 0285 qr(3,:)=log(qr(2,:)./(1-qr(2,:))); % x position of target in std measure 0286 qr(4,:)=qr(1,:).*qr(2,:).*(1-qr(2,:)); % slope*"stddev" at threshold 0287 case 2 % cumulative gaussian model 0288 qr(3,:)=norminv(qr(2,:)); % x position of target in std measure 0289 qr(4,:)=qr(1,:).*normpdf(qr(3,:)); % slope*stddev at threshold 0290 otherwise 0291 error('Unrecognised psychometric model selection'); 0292 end 0293 if nr(9)< 1 || nr(9)>2 0294 error('Unrecognised cost function option'); 0295 end 0296 mq=repmat(mq0,[1,1,3]); % initial means, joint mode and marginal mode all are equal 0297 vq=[var(xq,1,1); zeros(1,ni); var(sq,1,1)]; % initial variances (actually ignores prior probs) 0298 mqr=mq; % robust means and modes 0299 vqr=vq; % robust variances 0300 nresr=0; % robust calculation time 0301 xlim=[-Inf(1,ni); Inf(1,ni)]; % SNR limits for outliers 0302 % hn=[1-nr(4) 0 nr(4)]*vq; % artificially high value of cost function ensures all models are probed early on 0303 hn=Inf(1,ni); % very high initial cost function 0304 hfact=nr(8)^(1/ni); % growth factor to ensure that no model is neglected for too long 0305 xn=mq0(1,:); % start at mean value 0306 xmm=[xn; xn]; % min/max of probe values for each model 0307 if nargin>=4 && ~isempty(xp) 0308 if iscell(xp) 0309 xz=xp; 0310 else 0311 xz=repmat(num2cell(xp(:)',2),ni,1); % else replicate for each model 0312 end 0313 for i=1:ni 0314 [dum,j]=min(abs(xz{i}-mq0(1,i))); % select the closest available probe to the mean 0315 xn(i)=xz{i}(j); 0316 end 0317 else 0318 xz=cell(ni,1); % make empty cells 0319 end 0320 elseif iq>0 0321 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0322 % iq>0 means update or plot model iq % 0323 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0324 % first check the plotting options % 0325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0326 if nargin<3 0327 if nargin==2 0328 po=x; % plot options are 2nd argument 0329 elseif ~nargout 0330 po='p'; 0331 else 0332 po=''; 0333 end 0334 else 0335 if nargin>3 0336 po=xp; % plot options are 4th argument 0337 else 0338 po=''; 0339 end 0340 end 0341 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0342 % now get parameters of current model (iq) % 0343 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0344 nxq=nr(1); % number of SNR values 0345 nsq=nr(2); % number of slope values 0346 nxh=nr(3); % number of probe-SNR values 0347 nsxq=nxq*nsq; % size of log pdf array 0348 thresh=pr(1,iq); % target probability at threshold 0349 guess=pr(3,iq); % guess rate (1/choices) 0350 pscale=qr(1,iq); % prob range left after subtracting miss and guess probs 0351 xtstd=qr(3,iq); % x position of target in std measure 0352 xqi=xq(:,iq); % SRT values of the pdf array 0353 sqi=sq(:,iq); % slope (or log slope) values in PDF 0354 wqi=wq(:,iq); % log pdf array 0355 mqi= mq(:,iq,1); % [xe; se] means 0356 vqi=vq(:,iq); % [xv; sxv; sv] covariance matrix 0357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0358 % rescale the pdfs if necessary % 0359 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0360 ssd=sqrt(vqi(3)); % std deviation of slope 0361 xsd=sqrt(vqi(1)); % std deviation of SRT 0362 % determine range of new grid 0363 xqrange=xqi(end)-xqi(1); % range of old SRT grid 0364 sqrange=sqi(end)-sqi(1); % range of old slope grid 0365 if nresq(iq)<2 % keep old limits if nresq(iq)<2 0366 xq2lim=xqi([1 end]); 0367 sq2lim=sqi([1 end]); 0368 else 0369 xqsemirange=max(nr(7)*xsd,0.5*nr(20)*xqrange); 0370 xq2lim=[mqi(1)-xqsemirange, mqi(1)+xqsemirange]; 0371 sqsemirange=max(nr(7)*ssd,0.5*nr(20)*sqrange); 0372 sq2lim=[max(sqmin,mqi(2)-sqsemirange),mqi(2)+sqsemirange]; 0373 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0374 % To avoid unnecessary changes, keep old limits if the change % 0375 % is less than nr(17) times the previous range [nr(17)=0.1] % 0376 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0377 if abs(xq2lim(1)-xqi(1))<nr(17)*xqrange 0378 xq2lim(1)=xqi(1); 0379 end 0380 if abs(xq2lim(2)-xqi(end))<nr(17)*xqrange 0381 xq2lim(2)=xqi(end); 0382 end 0383 0384 if abs(sq2lim(1)-sqi(1))<nr(17)*sqrange 0385 sq2lim(1)=sqi(1); 0386 end 0387 if abs(sq2lim(2)-sqi(end))<nr(17)*sqrange 0388 sq2lim(2)=sqi(end); 0389 end 0390 end 0391 xq2=linspace(xq2lim(1),xq2lim(2),nxq)'; % new x axis values 0392 sq2=linspace(sq2lim(1),sq2lim(2),nsq)'; % new s axis values 0393 wqup=2; % update flag 0394 if xq2(1)<xqi(1) || xq2(end)>xqi(end) || sq2(1)<sqi(1) || sq2(end)>sqi(end) 0395 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0396 % if extrapolating, recalculate log-pdfs from saved data % 0397 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0398 if LOG, fprintf(LOG,'N=%d:%d, extrapolate: x(%.2g,%.2g)->(%.2g,%.2g) and s(%.2g,%.2g)->(%.2g,%.2g)\n',iq,nresq(iq),xqi([1 end]),xq2([1 end]),sqi([1 end]),sq2([1 end])); end 0399 wq2=repmat(pq0(2,iq)*(sq2-mq0(2,iq)).^2,1,nxq)+repmat(pq0(1,iq)*(xq2-mq0(1,iq))'.^2,nsq,1); % log probs of prior for model iq 0400 ires=find(res(1:nres,1)==iq); % find the results for this model 0401 if nr(11) % if log slope 0402 sqis=exp(sq2)/qr(4,iq); % inverse std dev of gaussian (proportional to slope) 0403 else 0404 sqis=sq2/qr(4,iq); % inverse std dev of gaussian (proportional to slope) 0405 end 0406 switch nr(10) % switch on model type 0407 case 1 0408 for i=1:length(ires) 0409 j=ires(i); 0410 r0=res(j,3)==0; 0411 wq2=wq2+log(r0+(1-2*r0)*(guess+pscale*((1+exp(sqis*xq2'-xtstd-repmat(sqis,1,nxq)*res(j,2))).^(-1)))); % P(l | r,x) 0412 end 0413 case 2 0414 for i=1:length(ires) 0415 j=ires(i); 0416 r0=res(j,3)==0; 0417 wq2=wq2+log(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,1,nxq)*res(j,2)-sqis*xq2'+xtstd))); % P(l | r,x) 0418 end 0419 end 0420 else % possibly do quadratic interpolation 0421 wq2=max(reshape(wqi,nsq,nxq)-max(wqi(:)),wfl); % turn into a normalized, clipped matrix for easy interpolation 0422 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0423 % use quadratic interpolation in SRT axis % 0424 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0425 if (xq2(end)-xq2(1))/(xqrange)>nr(16) % if range has shrunk by < nr(16), leave the SRT axis unchanged 0426 xq2=xqi; % copy the old SRT axis 0427 wqup=1; % update flag 0428 else % else do quadratic interpolation 0429 if LOG, fprintf(LOG,'N=%d:%d, interpolate x: x(%.2g,%.2g)->(%.2g,%.2g)\n',iq,nresq(iq),xqi([1 end]),xq2([1 end])); end 0430 xqf=(xq2-xqi(1))/(xqi(2)-xqi(1)); 0431 xqj=ceil(xqf); 0432 xqf=xqj-xqf; 0433 xqg=1-xqf; 0434 xqh=0.25*xqf.*xqg; % Quadratic coefficient 0435 xqf((xqj<=0) | (xqj>nxq))=0; 0436 xqg((xqj<0) | (xqj>=nxq))=0; 0437 xqh((xqj<1) | (xqj>=nxq-1))=0; 0438 wq2=wq2(:,min(max(xqj,1),nxq)).*repmat(xqf'+xqh',nsq,1)+wq2(:,min(max(xqj+1,1),nxq)).*repmat(xqg'+xqh',nsq,1) - ... 0439 (wq2(:,min(max(xqj-1,1),nxq))+wq2(:,min(max(xqj+2,1),nxq))).*repmat(xqh',nsq,1); 0440 end 0441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0442 % use quadratic interpolation in slope axis % 0443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0444 if (sq2(end)-sq2(1))/(sqrange)>nr(16) % if range has shrunk by < nr(16), leave the slope axis unchanged 0445 sq2=sqi; % copy the old slope axis 0446 wqup=wqup-1; % update flag 0447 else 0448 if LOG, fprintf(LOG,'N=%d:%d, interpolate s: s(%.2g,%.2g)->(%.2g,%.2g)\n',iq,nresq(iq),sqi([1 end]),sq2([1 end])); end 0449 sqf=(sq2-sqi(1))/(sqi(2)-sqi(1)); 0450 sqj=ceil(sqf); 0451 sqf=sqj-sqf; 0452 sqg=1-sqf; 0453 sqh=0.25*sqf.*sqg; % Quadratic coefficient 0454 sqf((sqj<=0) | (sqj>nsq))=0; 0455 sqg((sqj<0) | (sqj>=nsq))=0; 0456 sqh((sqj<1) | (sqj>=nsq-1))=0; 0457 wq2=wq2(min(max(sqj,1),nsq),:).*repmat(sqf+sqh,1,nxq)+wq2(min(max(sqj+1,1),nsq),:).*repmat(sqg+sqh,1,nxq) - ... 0458 (wq2(min(max(sqj-1,1),nsq),:)+wq2(min(max(sqj+2,1),nsq),:)).*repmat(sqh,1,nxq); 0459 end 0460 end 0461 if wqup>0 % check if we have updated either axis 0462 wq2=max(wq2(:)-max(wq2(:)),wfl); % turn back into a normalized, clipped vector 0463 sqi=sq2; % update slope (or log slope) values in PDF 0464 xqi=xq2; % update SRT values of the pdf array 0465 wqi=wq2; % log pdf array 0466 sq(:,iq)=sqi; % save new s axis (log slope or slope) 0467 xq(:,iq)=xqi; % save new x axis (SRT) 0468 wq(:,iq)=wqi; % save log-pdf 0469 end 0470 if nr(11) % if log slope 0471 sqis=exp(sqi)/qr(4,iq); % inverse std dev of gaussian (proportional to slope) 0472 else 0473 sqis=sqi/qr(4,iq); % inverse std dev of gaussian (proportional to slope) 0474 end 0475 if nargin>=3 0476 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0477 % update pdfs with a new probe result % 0478 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0479 if nres>=size(res,1) % ensure there is space for a new result 0480 res=[res; zeros(nres,7)]; % double the size of the results array 0481 end 0482 nres=nres+1; % increment the total number of results 0483 nresq(iq)=nresq(iq)+1; % increment the number of results for this model 0484 res(nres,1:3)=[iq,x,r]; % save the latest result 0485 xmm(:,iq)=[min(x,xmm(1,iq)); max(x,xmm(2,iq))]; % update min/max probe values 0486 hn=hn*hfact; % increase v_entropy gains to ensure all models get a chance 0487 0488 % disp(res(1:nres,:)); % ********** diagnostic 0489 0490 % 0491 % update log probabilities with the previous test result 0492 % If r==1, we multiply (add in log) by a horizontally reflected version of the psychometric function that equals 0493 % the threshold (e.g. 0.5 or 0.75) at the probe SNR, x, for all slopes. 0494 % If r==0, we multiply (add in log) by 1- this reflected version which therefore equals (1-thresh) at x. 0495 % 0496 r0=r==0; 0497 switch nr(10) % switch on model type 0498 case 1 0499 wqi=wqi+log(r0+(1-2*r0)*(guess+pscale*((1+exp(reshape(sqis*xqi'-xtstd,nsxq,1)-repmat(sqis,nxq,1)*x)).^(-1)))); % P(l | r,x) 0500 case 2 0501 wqi=wqi+log(r0+(1-2*r0)*(guess+pscale*normcdf(repmat(sqis,nxq,1)*x-reshape(sqis*xqi'-xtstd,nsxq,1)))); % P(l | r,x) 0502 otherwise 0503 error('Unrecognised psychometric model selection'); 0504 end 0505 wq(:,iq)=wqi; % save updated probabilities 0506 end 0507 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0508 % Calculate mean and covariance and v_entropy % 0509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0510 ewqi=exp(wqi-max(wqi)); % unnormalized probability vector 0511 ewqi=ewqi/sum(ewqi); 0512 wqsx=reshape(ewqi,nsq,nxq); % normalized probabilities 0513 px=sum(wqsx,1); % p(x0) 0514 ps=sum(wqsx,2); % p(s0) 0515 xe=px*xqi; % E(x0) 0516 se=ps'*sqi; % E(s0) 0517 0518 [pxpk,xm]=max(px); % marginal mode in x 0519 if xm>1 && xm<nxq % use quadratic interpolation in log prob if possible 0520 [dum,xm2]=quadpeak(log(px(xm-1:xm+1))'); 0521 xm=xm+xm2-2; 0522 end 0523 xm=(2-xm)*xqi(1)+(xm-1)*xqi(2); % marginal mode(x0) 0524 0525 [pspk,sm]=max(ps); 0526 if sm>1 && sm<nsq % use quadratic interpolation in log prob if possible 0527 [dum,sm2]=v_quadpeak(log(ps(sm-1:sm+1))); 0528 sm=sm+sm2-2; 0529 end 0530 sm=(2-sm)*sqi(1)+(sm-1)*sqi(2); % marginal mode(s0) 0531 0532 [wqpk,j]=max(ewqi); 0533 i=1+floor((j-1)/nsq); 0534 j=j-nsq*(i-1); 0535 if i>1 && i<nxq && j>1 && j<nsq % use quadratic interpolation in log prob if possible 0536 [dum,ji]=quadpeak(wqi(repmat((j-1:j+1)',1,3)+repmat(nsq*(i-2:i),3,1))); 0537 i=i+ji(2)-2; 0538 j=j+ji(1)-2; 0539 end 0540 xj=(2-i)*xqi(1)+(i-1)*xqi(2); % joint mode x 0541 sj=(2-j)*sqi(1)+(j-1)*sqi(2); % joint mode: s 0542 xv=px*(xqi.^2)-xe^2; % Var(x0) 0543 sv=ps'*(sqi.^2)-se^2; % Var(s0) 0544 sxv=ewqi'*(repmat(sqi-se,nxq,1).*reshape(repmat(xqi'-xe,nsq,1),nsxq,1)); % Cov(s0*x0) 0545 mq(:,iq,1)=[xe; se]; % save means 0546 mq(:,iq,2)=[xj; sj]; % save joint mode 0547 mq(:,iq,3)=[xm; sm]; % save marginal modes 0548 vq(:,iq)=[xv; sxv; sv]; % save covariance matrix 0549 xh=(px*log(px)')*(xqi(1)-xqi(2)); % differential v_entropy h(x0) 0550 sh=(ps'*log(ps))*(sqi(1)-sqi(2)); % differential v_entropy h(s0) 0551 if nargin==3 0552 switch nr(9) % cost function: 1=variance, 2=v_entropy 0553 case 1 0554 res(nres,4:7)=[xe se xv sv]; % save info for plotting history 0555 case 2 0556 res(nres,4:7)=[xe se xh sh]; % find the minimum of cost function 0557 end 0558 end 0559 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0560 % now calculate robust mean and covariance for all models % 0561 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0562 if nargout>4 || ~isempty(po) % need to update robust statistics if >4 outputs or if plotting 0563 if nr(19)==0 % if no outliers allowed in any case 0564 mr=mq; 0565 vr=vq; 0566 elseif nresr==nres % if robust statistics are all up-to-date 0567 mr=mqr; 0568 vr=vqr; 0569 else % update robust statistics for all models 0570 for jq=1:length(nresq) % loop for each model 0571 if any(res(nresr+1:nres,1)==jq) % if we have had any recent probe results for this model 0572 guessj=pr(3,jq); % guess rate (1/choices) 0573 pscalej=qr(1,jq); % prob range left after subtracting miss and guess probs 0574 xtstdj=qr(3,jq); % x position of target in std measure 0575 xqj=xq(:,jq); % SRT values of the pdf array 0576 sqj=sq(:,jq); % slope (or log slope) values in PDF 0577 % first determine the SNR range to include 0578 xej=mq(1,jq,1); % mean SRT 0579 sej=mq(2,jq,1); % mean slope or log slope 0580 resj=res(res(:,1)==jq,:); % results for this model 0581 if nr(11) % if log slope 0582 sqisj=exp(sej)/qr(4,jq); % inverse std dev of gaussian (proportional to slope) at mean slope 0583 else 0584 sqisj=sej/qr(4,jq); % inverse std dev of gaussian (proportional to slope) at mean slope 0585 end 0586 % find upper bound for SNR with r=0 0587 switch nr(10) % switch on model type 0588 case 1 0589 xlim(2,jq)=xej-(xtstdj+log(nr(19)/(pscalej-nr(19))))/sqisj; % outlier if r=0 and x > xlim 0590 case 2 0591 xlim(2,jq)=xej+(norminv((pscalej-nr(19))/pscalej)-xtstdj)/sqisj; % outlier if r=1 and x < xlim 0592 end 0593 mou=resj(:,2)>=xlim(2,jq) & resj(:,3)==0; % find high outliers 0594 if nr(19)>guessj % find lower bound for SNR with r=1 0595 switch nr(10) % switch on model type 0596 case 1 0597 xlim(1,jq)=xej-(xtstdj+log(pscalej/(nr(19)-guessj)-1))/sqisj; % outlier if r=1 and x < xlim 0598 case 2 0599 xlim(1,jq)=xej+(norminv((nr(19)-guessj)/pscalej)-xtstdj)/sqisj; % outlier if r=1 and x < xlim 0600 end 0601 mou=mou | resj(:,2)<=xlim(1,jq) & resj(:,3)==1; % add in low outliers 0602 end 0603 if any(mou) % if any results must be excluded 0604 % now recalculate pdf excluding bad points 0605 if nr(11) % if log slope 0606 sqisj=exp(sqj)/qr(4,jq); % inverse std dev of gaussian (proportional to slope) for all grid points 0607 else 0608 sqisj=sqj/qr(4,jq); % inverse std dev of gaussian (proportional to slope) for all grid points 0609 end 0610 wqj=repmat(pq0(2,jq)*(sqj-mq0(2,jq)).^2,1,nxq)+repmat(pq0(1,jq)*(xqj-mq0(1,jq))'.^2,nsq,1); % log probs of prior 0611 switch nr(10) % switch on model type 0612 case 1 % nr(10)=1: logistic model 0613 for j=find(mou==0)' 0614 r0=resj(j,3)==0; 0615 wqj=wqj+log(r0+(1-2*r0)*(guessj+pscalej*((1+exp(sqisj*xqj'-xtstdj-repmat(sqisj,1,nxq)*resj(j,2))).^(-1)))); % P(l | r,x) 0616 end 0617 case 2 % nr(10)=2: cumulative gaussian model 0618 for j=find(mou==0)' 0619 r0=resj(j,3)==0; 0620 wqj=wqj+log(r0+(1-2*r0)*(guessj+pscalej*normcdf(repmat(sqisj,1,nxq)*resj(j,2)-sqisj*xqj'+xtstdj))); % P(l | r,x) 0621 end 0622 end 0623 % now calculate robust means and modes 0624 wqj=wqj(:); % turn 2D array into a column vector 0625 ewqj=exp(wqj-max(wqj)); % unnormalized probability vector 0626 ewqj=ewqj/sum(ewqj); 0627 wqsxr=reshape(ewqj,nsq,nxq); % normalized probabilities 0628 pxr=sum(wqsxr,1); % p(x0) 0629 psr=sum(wqsxr,2); % p(s0) 0630 xer=pxr*xqj; % E(x0) 0631 ser=psr'*sqj; % E(s0) 0632 [pxpk,xmr]=max(pxr); % marginal mode in x 0633 if xmr>1 && xmr<nxq % use quadratic interpolation in log prob if possible 0634 [dum,xm2]=quadpeak(log(pxr(xmr-1:xmr+1))'); 0635 xmr=xmr+xm2-2; 0636 end 0637 xmr=(2-xmr)*xqj(1)+(xmr-1)*xqj(2); % marginal mode(x0) (requires the xqi to be uniformly spaced) 0638 [pspk,smr]=max(psr); 0639 if smr>1 && smr<nsq % use quadratic interpolation in log prob if possible 0640 [dum,sm2]=v_quadpeak(log(psr(smr-1:smr+1))); 0641 smr=smr+sm2-2; 0642 end 0643 smr=(2-smr)*sqj(1)+(smr-1)*sqj(2); % marginal mode(s0) 0644 [wqpk,j]=max(ewqj); % find max of 1-D array 0645 i=1+floor((j-1)/nsq); % convert to 2-D index (i,j) 0646 j=j-nsq*(i-1); 0647 if i>1 && i<nxq && j>1 && j<nsq % use quadratic interpolation in log prob if possible 0648 [dum,ji]=quadpeak(wqj(repmat((j-1:j+1)',1,3)+repmat(nsq*(i-2:i),3,1))); 0649 i=i+ji(2)-2; 0650 j=j+ji(1)-2; 0651 end 0652 xjr=(2-i)*xqj(1)+(i-1)*xqj(2); % joint mode x 0653 sjr=(2-j)*sqj(1)+(j-1)*sqj(2); % joint mode: s 0654 xvr=pxr*(xqj.^2)-xer^2; % Var(x0) 0655 svr=psr'*(sqj.^2)-ser^2; % Var(s0) 0656 sxvr=ewqj'*(repmat(sqj-ser,nxq,1).*reshape(repmat(xqj'-xer,nsq,1),nsxq,1)); % Cov(s0*x0) 0657 mqr(:,jq,:)=reshape([xer xjr xmr; ser sjr smr],[2 1 3]); % save means, joint modes and marginal modes 0658 vqr(:,jq)=[xvr; sxvr; svr]; % save covariance matrix 0659 else % if there are no outliers, just copy non-robust results 0660 mqr(:,jq,:)= mq(:,jq,:); % use non-robust means 0661 vqr(:,jq)= vq(:,jq); % and non-robust variances 0662 end 0663 end 0664 end 0665 nresr=nres; % all robust means/modes are now up-to-date 0666 mr=mqr; % send to output arguments 0667 vr=vqr; 0668 end 0669 if nr(11) % if using log-slope 0670 mr(2,:,:)=exp(mr(2,:,:)); % convert to real slope 0671 vr(2,:)=vr(2,:).*mr(2,:,1); % correct the covariance 0672 vr(3,:)=vr(3,:).*mr(2,:,1).^2; % and the slope variance 0673 end 0674 end 0675 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0676 % now estimate the next probe SNR % 0677 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0678 if ~numel(xz{iq}) % if no list of probe SNRs was specified 0679 ytry=exp(0.25i*pi*(0:7))'; % points around the circle 0680 ytry=[real(ytry) imag(ytry)]; 0681 [vtry,dtry]=eig([xv sxv; sxv sv]); % eigendecomposition of covariance matrix 0682 tryxs=repmat([xe,se],8,1)+nr(14)*ytry*sqrt(dtry)*vtry'; 0683 pmin=0.05; % target probe success probability 0684 if nr(11) 0685 tryxs(:,2)=qr(4,iq)*exp(-tryxs(:,2)); % convert log(slope) to std dev 0686 else 0687 tryxs(:,2)=qr(4,iq)*tryxs(:,2).^(-1); % convert slope to std dev 0688 end 0689 dp=nr(15); % maximum shift of probe value outside previous range 0690 switch nr(10) % switch by model type 0691 case 1 % logistic 0692 qmax=min(xmm(2,iq)+dp,max(tryxs(:,1)+(log((1-pmin)/pmin)-xtstd)*tryxs(:,2))); 0693 qmin=max(xmm(1,iq)-dp,min(tryxs(:,1)+(log(pmin/(1-pmin))-xtstd)*tryxs(:,2))); 0694 case 2 % cumulative gaussian 0695 qmax=min(xmm(2,iq)+dp, max(tryxs(:,1)+(norminv(1-pmin)-xtstd)*tryxs(:,2))); 0696 qmin=max(xmm(1,iq)-dp,min(tryxs(:,1)+(norminv(pmin)-xtstd)*tryxs(:,2))); 0697 end 0698 dxt=max(nr(5),(qmax-qmin)/nxh); % minimum step size of nr(5) [0.2 dB] 0699 xt=(qmin+qmax)/2+((1:nxh)-(1+nxh)/2)*dxt; 0700 else % if a specific list of probe SNRs exists 0701 xzi=xz{iq}; % xzi is the list of available probe SNRs 0702 if numel(xzi)<=nxh % use all available probe SNRs if there are not too many 0703 xt=xzi; 0704 else 0705 [xt,ixt]=min(abs(xzi-xe)); % find the closest one to xe ** not necessarily optimum *** 0706 ixt=max(1,min(1+numel(xzi)-nxh,ixt-floor((1+nxh)/2))); % arrange symmetrically around xt 0707 xt=xzi(ixt:min(ixt+nxh-1,numel(xzi))); 0708 end 0709 end 0710 nxhp=length(xt); % xt are the potential probe SNRs 0711 % Now find the probe value that minimizes the expected value of the cost function 0712 % In the following: l = parameters of psychometric function, x = the probe SNR and r = the probe result 0713 switch nr(10) 0714 case 1 0715 prt=guess+pscale*((1+exp(repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)-repmat(sqis,nxq,1)*xt)).^(-1)); % P(r=1 | l,x) 0716 case 2 0717 prt=guess+pscale*normcdf(repmat(sqis,nxq,1)*xt-repmat(reshape(sqis*xqi'-xtstd,nsxq,1),1,nxhp)); % P(r=1 | l,x) 0718 end 0719 wqt=repmat(ewqi,1,nxhp); 0720 hminr=zeros(2,1); % space for minimum expected cost function for each r0 0721 hminj=zeros(1,nxhp); % space for minimum expected cost function for each x0 0722 if nr(18) % if doing look 2-ahead 0723 pl1=prt.*wqt; % posterior prob given success = p(l | x0,r0=1) unnormalized 0724 pl0=wqt-pl1; % posterior prob given failure = p(l | x0,r0=0) unnormalized 0725 pr0=sum(pl1,1); % p(r0=1 | x0)=Sum{P(r0=1 | l,x0)*P(l)} [note each column of wqt is normalized] (row vector) 0726 plxr=reshape([pl0./repmat(1-pr0,nsxq,1) pl1./repmat(pr0,nsxq,1)],nsxq,nxhp,2); % posterior prob p(l | x0,r0) column-normalized 0727 nx0=nxhp; % outer loop for each x0 0728 nr0=2; % inner loop for each r0 0729 else % if only doing look 1-ahead 0730 nx0=1; % only execute outer loop once 0731 nr0=1; % only execute inner loop once 0732 pr0=0; % dummy value (used at end of outer loop) 0733 end 0734 hx2=zeros(nx0,nxhp,nr0); % space for square array of expected cost functions (for possible plotting only) 0735 for j=1:nx0 % loop for each possible probe SNR, x0, (or only once is look-1-ahead) 0736 for jr=1:nr0 % loop for each possible test result, r0=jr-1 (or only once is look-1-ahead) 0737 if nr(18) % if doing look 2-ahead 0738 wqt=repmat(plxr(:,j,jr),1,nxhp); % posterior prob p(l | x0=xt(j),r0=jr-1) column-normalized 0739 end 0740 pl1=prt.*wqt; % posterior prob given success = p(l | x,r=1) unnormalized 0741 pl0=wqt-pl1; % posterior prob given failure = p(l | x,r=0) unnormalized 0742 prh=sum(pl1,1); % p(r | x)=Sum{P(r | l,x)*P(l)} [note each column of wqtjr is normalized] (row vector) 0743 pl1=pl1./repmat(prh,nsxq,1); % posterior prob given success = p(l | x,r=1) normalized 0744 pl0=pl0./repmat(1-prh,nsxq,1); % posterior prob given failure = p(l | x,r=0) normalized 0745 0746 px1=squeeze(sum(reshape(pl1,nsq,nxq,[]),1)); % p(x0 | x,r=1) 0747 px0=squeeze(sum(reshape(pl0,nsq,nxq,[]),1)); % p(x0 | x,r=0) 0748 ps1=squeeze(sum(reshape(pl1,nsq,nxq,[]),2)); % p(s0 | x,r=1) 0749 ps0=squeeze(sum(reshape(pl0,nsq,nxq,[]),2)); % p(s0 | x,r=0) 0750 xet1=xqi'*px1; % E(x0 | x,r=1) 0751 xvt1=(xqi.^2)'*px1-xet1.^2; % Var(x0 | x,r=1) 0752 xet0=xqi'*px0; % E(x0 | x,r=0) 0753 xvt0=(xqi.^2)'*px0-xet0.^2; % Var(x0 | x,r=0) 0754 xvt=xvt1.*prh+xvt0.*(1-prh); % E(Var(x0 | x )) 0755 set1=sqi'*ps1; % E(s0 | x,r=1) 0756 svt1=(sqi.^2)'*ps1-set1.^2; % Var(s0 | x,r=1) 0757 set0=sqi'*ps0; % E(s0 | x,r=0) 0758 svt0=(sqi.^2)'*ps0-set0.^2; % Var(s0 | x,r=0) 0759 svt=svt1.*prh+svt0.*(1-prh); % E(Var(s0 | x )) 0760 xht1=sum(log(px1).*px1,1); % -H(x0 | x,r=1) 0761 xht0=sum(log(px0).*px0,1); % -H(x0 | x,r=0) 0762 xht=(xht1.*prh+xht0.*(1-prh))*(xqi(1)-xqi(2)); % h(x0 | x) 0763 sht1=sum(log(ps1).*ps1,1); % -H(s0 | x,r=1) 0764 sht0=sum(log(ps0).*ps0,1); % -H(s0 | x,r=0) 0765 sht=(sht1.*prh+sht0.*(1-prh))*(sqi(1)-sqi(2)); % h(s0 | x) 0766 switch nr(9) % cost function: 1=variance, 2=v_entropy 0767 case 1 0768 hx=(xvt + nr(4)*svt)/(1+nr(4)); % expected cost function for each possible test SNR 0769 case 2 0770 hx=(xht + nr(4)*sht)/(1+nr(4)); % expected cost function for each possible test SNR % find the minimum of cost function 0771 end 0772 hx2(j,:,jr)=hx; % save expected cost function for possible plotting 0773 [hminr(jr),ix]=min(hx); % find the minimum of cost function for each value of r (0 or 1) 0774 % fprintf('Probe range: %.3g %.3g; choose %.3g\n',xt(1),xt(end),xt(ix)); 0775 end 0776 hminj(j)=[1-pr0(j) pr0(j)]*hminr; % expected cost function for each x0 assuming x1 is chosen optimally 0777 end 0778 if nr(18) % if doing look 2-ahead 0779 [hminr(1),ix]=min(hminj); % find the minimum of cost function 0780 end 0781 % now ix indexes the optimal probe choice and hminr(1) is the resultant expected cost function 0782 xn(iq)=xt(ix); % save the next probe snr for this model 0783 % calculate the expected decrease in cost function (to decide which model to probe next) 0784 switch nr(9) % cost function: 1=variance, 2=v_entropy 0785 case 1 % minimize variance 0786 hn(iq)=(xv + nr(4)*sv)/(1+nr(4))-hminr(1); % save the expected decrease in cost function for this model 0787 case 2 % minimize v_entropy 0788 hn(iq)=(xh + nr(4)*sh)/(1+nr(4))-hminr(1); % save the expected decrease in cost function for this model 0789 end 0790 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0791 % now do plotting % 0792 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0793 if ~isempty(po) 0794 h=gcf; 0795 if isreal(h) 0796 fig0=round(h); % get figure number 0797 else 0798 fig0=get(h,'number'); % in new versions of matlab use this method 0799 end 0800 axm=[3 -2; 2 -1]; % matrix to extend an axis beackwards by two steps 0801 for i=1:length(po) 0802 poi=po(i); 0803 switch poi 0804 case 'p' % plot pdf 0805 nxq=nr(1); 0806 nsq=nr(2); 0807 xqj=[axm*xqi(1:2);xqi]; 0808 sqj=[axm*sqi(1:2);sqi]; 0809 figure(fig0); 0810 imagesc(xqj,sqj,[zeros(1,2) px/pxpk; zeros(1,nxq+2);ps/pspk zeros(nsq,1) wqsx/wqpk]); 0811 hold on 0812 plot(xe,se,'ko',xm,sm,'k^',xj,sj,'k*'); 0813 plot(xqj(2),se,'wo',xqj(2),sm,'w^',xqj(2),sj,'w*'); 0814 plot(xe,sqj(2),'wo',xm,sqj(2),'w^',xj,sqj(2),'w*'); 0815 if any(mqr(:,iq,:)~=mq(:,iq,:)) 0816 plot(mqr(1,iq,1),mqr(2,iq,1),'go',mqr(1,iq,3),mqr(2,iq,3),'g^',xj,sj,'g*'); 0817 plot(xqj(end),mqr(2,iq,1),'go',xqj(end),mqr(2,iq,3),'g^',xqj(end),mqr(2,iq,2),'g*'); 0818 plot(mqr(1,iq,1),sqj(end),'go',mqr(1,iq,3),sqj(end),'g^',mqr(1,iq,2),sqj(end),'g*'); 0819 title('Joint pdf: o mean, * mode, \Delta marg mode (green=robust)'); 0820 else 0821 title('Joint pdf: o mean, * mode, \Delta marginal mode'); 0822 end 0823 t=linspace(0,2*pi,200); 0824 xcir=cos(t); 0825 scir=sin(t); 0826 vcir=sqrt((sv*xcir.^2+xv*scir.^2-2*sxv*xcir.*scir)/(xv*sv-sxv^2)); 0827 plot(xe+xcir./vcir,se+scir./vcir,'w-'); 0828 hold off 0829 axis 'xy'; 0830 % colorbar; 0831 % v_cblabel('Relative Probability'); 0832 xlabel(sprintf('SNR @ %d%%SRT (dB)',round(pr(1)*100))); 0833 if nr(11) 0834 ylabel('Log psychometric Slope at threshold (prob/dB)'); 0835 else 0836 ylabel('Psychometric Slope at threshold (prob/dB)'); 0837 end 0838 fig0=fig0+1; 0839 case 'c' % plot cost evolution 0840 resi=res(res(:,1)==iq,:); 0841 nresi=size(resi,1); 0842 if nresi>0 0843 figure(fig0); 0844 plot(1:nresi,resi(:,6)/(1+nr(4)),'--r',1:nresi,(resi(:,6) + nr(4)*resi(:,6))/(1+nr(4)),'-b'); 0845 v_axisenlarge([-1 -1.03]); 0846 xlabel('Trial Number'); 0847 switch nr(9) % cost function: 1=variance, 2=v_entropy 0848 case 1 % minimize variance 0849 ylabel('Weighted Variance: SRT and Total'); 0850 case 2 % minimize v_entropy 0851 ylabel('Weighted Entropy: SRT and Total'); 0852 end 0853 title ('Cost function evolution'); 0854 fig0=fig0+1; 0855 end 0856 case 'h' % plot history 0857 resi=res(res(:,1)==iq,:); 0858 nresi=size(resi,1); 0859 if nresi>0 0860 if nr(11) 0861 resi(:,5)=exp(resi(:,5)); % convert to real slope 0862 ylim=xe+([-2 0 thresh 1 3]-thresh)/exp(se); 0863 else 0864 ylim=xe+([-2 0 thresh 1 3]-thresh)/se; 0865 end 0866 figure(fig0); 0867 plot(1:nresi,resi(:,[4 4 4])+resi(:,5).^(-1)*([0 thresh 1]-thresh),'-b'); 0868 hold on; 0869 % could plot outliers in a different colour by checking mou(res(:,1)==iq) 0870 msky=resi(:,3)==1; 0871 plot(find(msky),resi(msky,2),'r+'); 0872 msky=resi(:,3)==0; 0873 plot(find(msky),resi(msky,2),'bo'); 0874 gcaylim=get(gca,'ylim'); 0875 poutb=0; 0876 if nr(19)>0 && xlim(2,iq)<gcaylim(2) % plot upper outlier bound 0877 plot([1;nresi],[1;1]*xlim(2,iq),'-r'); 0878 v_texthvc(nresi,xlim(2,iq),'max o','rtr'); 0879 end 0880 if nr(19)>guess && xlim(1,iq)>gcaylim(1) % plot lower outlier bound 0881 plot([1;nresi],[1;1]*xlim(1,iq),'-r'); 0882 v_texthvc(nresi,xlim(1,iq),'min +','rbr'); 0883 end 0884 hold off; 0885 xlabel('Trial Number'); 0886 ylabel('Probe SNR (dB) [o,x=result]'); 0887 v_axisenlarge(-1.02); 0888 title('SRT mean + slope intercepts'); 0889 fig0=fig0+1; 0890 end 0891 case 'x' % plot probe choice 0892 figure(fig0); 0893 if nr(18) % if doing look 2-ahead 0894 xtaxm=xt(1:2)*axm'; 0895 hx2e=hx2(:,:,1).*repmat(1-pr0',1,nxhp)+hx2(:,:,2).*repmat(pr0',1,nxhp); 0896 [dum,ix2m]=min(hx2,[],2); % ix2 indices to minimize hx2 0897 % ========== note that the following line assumes that the xt() 0898 % are uniformly spaced (not necessarily true) 0899 imagesc([xtaxm xt],xt,[hminj' repmat(max(hx2e(:)),nxhp,1) hx2e]); 0900 axis xy; 0901 colorbar; 0902 hold on 0903 [hnmin,ii]=max(hn); % find the next model to probe 0904 plot(xtaxm(1),xn(iq),'>w'); 0905 if iq==ii 0906 plot(xt(1:2)*axm(1,:)',xn(iq),'*w'); % circle the next model to be probed 0907 end 0908 plot([xtaxm(1) xt(end)],xn(iq)*[1 1],':w'); 0909 plot(xt(ix2m(:,1,1)),xt,'ow'); 0910 plot(xt(ix2m(:,1,2)),xt,'+w'); 0911 hold off 0912 xlabel('SNR probe 2 (dB) [+,o=best for R1=1,0]'); 0913 ylabel(sprintf('SNR probe 1 (dB) [> = %.2f]',xn(iq))); 0914 title('Expected Cost after 2 probes'); 0915 else % if doing look 1-ahead only 0916 plot(xt,hx2(1,:,1),'b'); 0917 v_axisenlarge([-1 -1.03]); 0918 caylim=get(gca,'ylim'); 0919 hold on 0920 plot(xn([iq iq]),caylim,':k'); 0921 xlabel(sprintf('Next probe SNR (dB) [min @ %.2f]',xn(iq))); 0922 ylabel('Cost') 0923 title('Expected Cost after next probe'); 0924 end 0925 fig0=fig0+1; 0926 end 0927 end 0928 end 0929 end 0930 if iq~=0 0931 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0932 % now select the appropriate model to probe next % 0933 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0934 [hnmin,ii]=max(hn); % chose model with the biggest expected decrease 0935 xx=xn(ii); 0936 m=mq; 0937 v=vq; 0938 if nr(11) % if using log-slope 0939 m(2,:,:)=exp(m(2,:,:)); % convert to real slope 0940 v(2,:)=v(2,:).*m(2,:,1); % correct the covariance 0941 v(3,:)=v(3,:).*m(2,:,1).^2; % and the slope variance 0942 end 0943 else 0944 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0945 % iq=0 means output model parameters and probe results % 0946 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0947 xx=pr; 0948 ii=nr; 0949 m=res(1:nres,:); 0950 if ~nargout % print them if no output arguments to call 0951 fid=LOG+(LOG==0); % use standard out if LOG==0 0952 pdesc={'Threshold','Miss prob','Guess Prob','Min SNR','Max SNR','Min Slope','Max Slope'}; 0953 fprintf(LOG,'\n*********\nModel-specific Parameters\n'); 0954 for i=1:7 0955 fprintf(LOG,'%4d) %s: ',i,pdesc{i}); 0956 if size(pr,2)>1 0957 fprintf(LOG,'%.5g, ',pr(i,1:size(pr,2)-1)); 0958 end 0959 fprintf(LOG,'%.5g\n',pr(i,size(pr,2))); 0960 end 0961 qdesc={'nx SNR values in PDF', ... 0962 'ns Slope values in PDF', ... 0963 'nh Probe SNR values to evaluate', ... 0964 'cs Weighting of slope relative to SRT in cost function', ... 0965 'dh Min step size in dB for probe SNRs', ... 0966 'sl Min slope at threshold', ... 0967 'kp Std deviations of the pdfs to keep', ... 0968 'hg Amount to grow expected gains in ni trials', ... 0969 'cf Cost function', ... 0970 'pm Psychometric model', ... 0971 'lg Use log slope as parameter', ... 0972 'pm Prior std devs in semi-width', ... 0973 'pf Integrated probability floor', ... 0974 'ts Number of std devs to explore', ... 0975 'dp Max shift in probe value', ... 0976 'it Grid interpolation threshold', ... 0977 'at Grid boundary tolerance', ... 0978 'la Look 2-ahead when choosing probe', ... 0979 'op Oulier probability threshold', ... 0980 'rx min grid shrink factor'}; 0981 qoptf=[9 10]; % fields with options 0982 qoptval={'variance','entropy'; ... 0983 'logistic','cumulative gaussian'}; 0984 fprintf(LOG,'\nShared Parameters\n'); 0985 for i=1:length(nr) 0986 fprintf(LOG,'%4d) %s: ',i,qdesc{i}); 0987 j=find(qoptf==i,1); 0988 if numel(j) 0989 fprintf(LOG,'%d=%s\n',nr(i),qoptval{j,nr(i)}); 0990 else 0991 fprintf(LOG,'%.5g\n',nr(i)); 0992 end 0993 end 0994 fprintf(LOG,'\n'); 0995 end 0996 end