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