# v_psychofunc

## PURPOSE

V_PSYCHOFUNC Calculate psychometric functions: trial success probability versus SNR

## SYNOPSIS

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

## DESCRIPTION

    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

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_psycdigit V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli

## SOURCE CODE

    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 Lesser General Public License as published by
    0056 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0063 %
    0064 %   You can obtain a copy of the GNU Lesser General Public License from
    0065 %   https://www.gnu.org/licenses/ .
    0066 %
    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