V_PSYCHOFUNC Calculate psychometric functions: trial success probability versus SNR Usage: p=v_psychofunc('',q,x) % calculate probabilities b=v_psychofunc('r',q,x) % generate boolean variables with success prob p p=v_psychofunc(m,q,x,r) % Calculate likelihoods for observations r x=v_psychofunc([m 'i'],q,p) % Calculate inverse Inputs: m mode string [may be omitted if not required] 'n' do not normalize likelihoods 'f' do not squeeze output arrays to remove singleton dimensions 'i' calculate inverse function 'r' calculate binary random variables with probability p ['s' calculate sweet points for threshold and slope] ['d' calculate partial derivatives with respect to q(1:5)] 'g' plot graph 'G' plot image 'c' include colourbar q model parameters. Either a column vector with a single model, a matrix with one model per column or a cell array with multiple values for some or all of the parameters 1 probability at threshold [0.5] 2 threshhold [0 dB] 3 slope at threshold [0.1 prob/dB ] 4 miss or lapse probability [0] 5 guess probability [0] 6 psychometric function type [1] 1 = logistic 2 = cumulative Gaussian 3 = Weibull [4 = reversed Weibull] [5 = Gumbell] [6 = reversed Gumbell] x vector of SNR values r test results (0 or 1) corresponding to x p vector of probabilities Outputs: p array of probabilities or random variates ('r' option). p is a squeezed 7-dimensional array whose dimensions correspond to x followed by the 6 model parameter entries. if q is a cell array, singleton dimensions are removed unless the 'f' option is given. x Inverse function gives SNR, x, as a function of p b array of boolean variables
0001 function p=v_psychofunc(m,q,x,r) 0002 %V_PSYCHOFUNC Calculate psychometric functions: trial success probability versus SNR 0003 % 0004 % Usage: p=v_psychofunc('',q,x) % calculate probabilities 0005 % b=v_psychofunc('r',q,x) % generate boolean variables with success prob p 0006 % p=v_psychofunc(m,q,x,r) % Calculate likelihoods for observations r 0007 % x=v_psychofunc([m 'i'],q,p) % Calculate inverse 0008 % 0009 % Inputs: 0010 % m mode string [may be omitted if not required] 0011 % 'n' do not normalize likelihoods 0012 % 'f' do not squeeze output arrays to remove singleton dimensions 0013 % 'i' calculate inverse function 0014 % 'r' calculate binary random variables with probability p 0015 % ['s' calculate sweet points for threshold and slope] 0016 % ['d' calculate partial derivatives with respect to q(1:5)] 0017 % 'g' plot graph 0018 % 'G' plot image 0019 % 'c' include colourbar 0020 % q model parameters. Either a column vector with a single model, 0021 % a matrix with one model per column or a cell array with multiple values for 0022 % some or all of the parameters 0023 % 1 probability at threshold [0.5] 0024 % 2 threshhold [0 dB] 0025 % 3 slope at threshold [0.1 prob/dB ] 0026 % 4 miss or lapse probability [0] 0027 % 5 guess probability [0] 0028 % 6 psychometric function type [1] 0029 % 1 = logistic 0030 % 2 = cumulative Gaussian 0031 % 3 = Weibull 0032 % [4 = reversed Weibull] 0033 % [5 = Gumbell] 0034 % [6 = reversed Gumbell] 0035 % x vector of SNR values 0036 % r test results (0 or 1) corresponding to x 0037 % p vector of probabilities 0038 % 0039 % Outputs: 0040 % p array of probabilities or random variates ('r' option). 0041 % p is a squeezed 7-dimensional array 0042 % whose dimensions correspond to x followed by the 6 model parameter entries. 0043 % if q is a cell array, singleton dimensions are removed unless the 'f' option is given. 0044 % x Inverse function gives SNR, x, as a function of p 0045 % b array of boolean variables 0046 0047 % Copyright (C) Mike Brookes 2009-2010 0048 % Version: $Id: v_psychofunc.m 10865 2018-09-21 17:22:45Z dmb $ 0049 % 0050 % VOICEBOX is a MATLAB toolbox for speech processing. 0051 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0052 % 0053 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0054 % This program is free software; you can redistribute it and/or modify 0055 % it under the terms of the GNU General Public License as published by 0056 % the Free Software Foundation; either version 2 of the License, or 0057 % (at your option) any later version. 0058 % 0059 % This program is distributed in the hope that it will be useful, 0060 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0061 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0062 % GNU General Public License for more details. 0063 % 0064 % You can obtain a copy of the GNU General Public License from 0065 % http://www.gnu.org/copyleft/gpl.html or by writing to 0066 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0068 0069 % first sort out input arguments 0070 minp=0.01; % minimum probability to use for inverse function by default 0071 qq=[0.5 0 0.1 0 0 1]'; % default values for q 0072 if nargin<4 0073 r=[]; 0074 if nargin<3 0075 x=[]; 0076 if nargin<2 0077 q=[]; 0078 if ~nargin 0079 m=''; 0080 end 0081 end 0082 end 0083 end 0084 if ~ischar(m); % mode argument is optional 0085 r=x; 0086 x=q; 0087 q=m; 0088 m=''; 0089 end 0090 sq=size(q); 0091 ckmod=0; 0092 if iscell(q) 0093 nq=ones(1,6); 0094 qax=num2cell([0; qq]); % used for plotting 0095 for i=1:min(numel(q),6) 0096 nq(i)=numel(q{i}); 0097 if nq(i)>=1 0098 nr=size(qq,2); 0099 qax{i+1}=q{i}; 0100 if i<=5 % do not replicate for multiple models 0101 qq=repmat(qq,1,nq(i)); 0102 qq(i,:)=reshape(repmat(q{i}(:)',nr,1),1,nr*nq(i)); 0103 else 0104 qq(i,:)=repmat(q{i}(1),1,nr); 0105 end 0106 end 0107 end 0108 nq=max(nq,1); 0109 nmod=nq(6); 0110 if nmod>1 % list of models to use 0111 modlist=q{6}; 0112 else 0113 modlist=qq(6,1); % default model 0114 end 0115 else 0116 nq=sq(2); 0117 if nq 0118 ql=repmat(qq,1,nq); 0119 ql(1:sq(1),:)=q; 0120 else 0121 ql=qq; 0122 nq=1; 0123 end 0124 modlist=unique(ql(6,:)); 0125 nmod=length(modlist); 0126 ckmod=nmod>1; % need to check model list 0127 qq=ql; 0128 end 0129 % now perform the calculation 0130 nx=numel(x); 0131 npt=50; % number of points 0132 if any(m=='i') % doing inverse 0133 if ~nx 0134 nx=npt; 0135 xlim=[max(qq(5,:)),1-max(qq(4,:))]*[1-minp minp; minp 1-minp]; 0136 x=linspace(xlim(1),xlim(2),nx)'; 0137 end 0138 p=zeros([nx nq]); % space for SNRs 0139 ia=0; 0140 for i=1:nmod % loop for each model type 0141 mod=modlist(i); 0142 if ckmod 0143 qq=ql(:,ql(6,:)==mod); 0144 end 0145 pscale=1-qq(4,:)-qq(5,:); 0146 pstd=(qq(1,:)-qq(5,:))./pscale; % prob target compensating for miss and lapse probs 0147 sstd=qq(3,:)./pscale; % slope compensating for miss and lapse probs 0148 px=x(:)*pscale.^(-1)-repmat(qq(5,:)./pscale,nx,1); % adjust for miss and lapse probs 0149 switch mod 0150 case 1 0151 beta=sstd./(pstd.*(1-pstd)); 0152 % alpha=qq(2,:)+log((1-pstd)./pstd)./beta; 0153 px=repmat(qq(2,:)+log((1-pstd)./pstd)./beta,nx,1)-log(px.^(-1)-1).*repmat(beta.^(-1),nx,1); 0154 case 2 % cumulative Gaussian function 0155 xtstd=norminv(pstd); % x position of target in std measure 0156 sig=normpdf(xtstd)./sstd; 0157 px= repmat(qq(2,:)-sig.*xtstd,nx,1) + repmat(sig,nx,1).*norminv(px); 0158 case 3 0159 wlog=log(1-pstd); 0160 kbeta=sstd./((pstd-1).*wlog); 0161 alpha=qq(2,:)-log(-wlog)./kbeta; 0162 px=repmat(alpha,nx,1)+log(-log(1-px)).*repmat(kbeta.^(-1),nx,1); 0163 otherwise 0164 error('Invalid psychometric model index'); 0165 end 0166 if ckmod 0167 p(:,ql(6,:)==i)=px; 0168 else 0169 ib=ia+numel(p)/nmod; 0170 p(ia+1:ib)=px(:); 0171 ia=ib; 0172 end 0173 end 0174 else % doing forward mapping 0175 if ~nx 0176 ef=2; % expansion factor 0177 nx=npt; 0178 x=linspace(min(qq(2,:)-ef*(qq(1,:)-qq(5,:))./qq(3,:)), ... 0179 max(qq(2,:)+ef*(1-qq(1,:)-qq(4,:))./qq(3,:)),nx)'; 0180 end 0181 p=zeros([nx nq]); % space for probabilities 0182 ia=0; 0183 for i=1:nmod % loop for each model type 0184 mod=modlist(i); 0185 if ckmod 0186 qq=ql(:,ql(6,:)==mod); 0187 end 0188 pscale=1-qq(4,:)-qq(5,:); % prob range excluding miss and lapse probs 0189 pstd=(qq(1,:)-qq(5,:))./pscale; % prob target compensating for miss and lapse probs 0190 sstd=qq(3,:)./pscale; % slope compensating for miss and lapse probs 0191 switch mod 0192 case 1 % logistic function 0193 beta=sstd./(pstd.*(1-pstd)); 0194 % alpha=qq(2,:)+log((1-pstd)./pstd)./beta; 0195 px=(1+exp(repmat(beta.*qq(2,:)+log((1-pstd)./pstd),nx,1)-x(:)*beta)).^(-1); 0196 case 2 % cumulative Gaussian function 0197 xtstd=norminv(pstd); % x position of target in std measure 0198 sigi=sstd./normpdf(xtstd); 0199 px=normcdf(x(:)*sigi-repmat(qq(2,:).*sigi-xtstd,nx,1)); 0200 case 3 0201 wlog=log(1-pstd); 0202 kbeta=sstd./((pstd-1).*wlog); 0203 alpha=qq(2,:)-log(-wlog)./kbeta; 0204 px=1-exp(-exp(x(:)*kbeta-repmat(alpha.*kbeta,nx,1))); 0205 otherwise 0206 error('Invalid psychometric model index'); 0207 end 0208 px=repmat(qq(5,:),nx,1)+repmat(pscale,nx,1).*px; % adjust for miss and lapse probs 0209 if ckmod 0210 p(:,ql(6,:)==i)=px; 0211 else 0212 ib=ia+numel(p)/nmod; 0213 p(ia+1:ib)=px(:); 0214 ia=ib; 0215 end 0216 end 0217 if numel(r) % we are calculating likelihoods 0218 mk=r(:)==0; 0219 p(mk,:)=1-p(mk,:); % invert probability for results that are zero 0220 if nx>1 0221 if any(m=='n') 0222 p=prod(p,1); 0223 else 0224 p=sum(log(p),1); 0225 p=exp(p-max(p(:))); 0226 p=p/sum(p(:)); % normalize to equal 1 0227 end 0228 nx=1; 0229 end 0230 end 0231 0232 end 0233 pg=p; % save unsqueezed p for plotting 0234 if ~any(m=='f') && iscell(q) % remove all singleton dimensions 0235 szp=size(p); 0236 szq=szp(szp>1); 0237 szq=[szq ones(1,max(0,2-numel(szq)))]; 0238 p=reshape(p,szq); 0239 end 0240 if any(m=='r') && ~any(m=='i'); 0241 p=rand(size(p))<p; 0242 end 0243 0244 if ~nargout || any(lower(m)=='g') 0245 clf; 0246 szp=[nx nq]; 0247 czp=sum(szp>1); 0248 if czp>0 % check if there is anything to plot 0249 if iscell(q) 0250 axlab={'Input SNR','Threshold prob','Threshold SNR','Threshold slope','Lapse prob','Guess prob','Sigmoid type'}; 0251 [szs,izs]=sort(szp,'descend'); 0252 pg=permute(pg,izs); 0253 qax{1}=x; 0254 if any(m=='G') || czp>2 % image 0255 ngr=prod(szs(3:end)); 0256 ncol=ceil(sqrt(ngr)); 0257 nrow=ceil(ngr/ncol); 0258 npix=szs(1)*szs(2); 0259 ia=0; 0260 for i=1:ngr 0261 subplot(nrow,ncol,i); 0262 ib=ia+npix; 0263 imagesc(qax{izs(1)},qax{izs(2)},reshape(pg(ia+1:ib),szs(1:2))'); 0264 axis 'xy' 0265 if any(m=='c') 0266 colorbar; 0267 end 0268 if nrow*ncol-i<ncol 0269 xlabel(axlab(izs(1))); 0270 end 0271 if rem(i-1,ncol)==0 0272 ylabel(axlab(izs(2))); 0273 end 0274 ia=ib; 0275 end 0276 else % graph 0277 plot(qax{izs(1)},reshape(permute(pg,izs),szs(1:2)),'-'); 0278 xlabel(axlab{izs(1)}); 0279 end 0280 else 0281 if any(m=='G') % image 0282 imagesc(pg'); 0283 axis 'xy' 0284 if any(m=='c') 0285 colorbar; 0286 end 0287 xlabel('Input SNR (dB)'); 0288 ylabel('Model Index'); 0289 else % graph 0290 if nx>=nq 0291 plot(x,pg,'-'); 0292 xlabel('Input SNR (dB)'); 0293 else 0294 plot(1:nq,pg','-'); 0295 xlabel('Model Index'); 0296 end 0297 end 0298 end 0299 end 0300 end 0301 0302