Home > voicebox > psychofunc.m

psychofunc

PURPOSE ^

Calculate psychometric functions: trial success probability versus SNR

SYNOPSIS ^

function p=psychofunc(m,q,x,r)

DESCRIPTION ^

 Calculate psychometric functions: trial success probability versus SNR

 Usage: p=psychofunc('',q,x)       % calculate probabilities
        b=psychofunc('r',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
                      ['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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function p=psychofunc(m,q,x,r)
0002 % Calculate psychometric functions: trial success probability versus SNR
0003 %
0004 % Usage: p=psychofunc('',q,x)       % calculate probabilities
0005 %        b=psychofunc('r',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 %                      ['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: psychofunc.m 9302 2017-01-18 16:19:20Z 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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003