# v_psycestu

## PURPOSE

V_PSYCESTU estimate unimodal psychometric function

## SYNOPSIS

function [xx,ii,m,v]=v_psycestu(iq,x,r,xp)

## DESCRIPTION

    V_PSYCESTU estimate unimodal psychometric function

     Usage: [xx,ii,m,v]=v_psycestu(-n,p,q,xp) % initialize n models
            [xx,ii,m,v]=v_psycestu(i,x,r)     % supply a trial result to v_psycest
                        v_psycestu(i)         % plot pdf of model i
                  [p,q]=psychestu(0)        % output model parameters (or print them if no outputs)

     Inputs:
             -n        minus the number of models
              p(:,n)   parameters for each model
                          1  miss [0.04]
                          2  guess [0.1]
                          3,4  SNR at peak: [min max] [-20 20]
                          5,6  normalized semi-width: [min max] [0 20]
                          7,8  low side slope: [min max] [0.03 0.3]
                          9,10  high side slope: [min max] [0.03 0.3]
              q(:)     parameters common to all models (vector or struct)
                          1  q.nxslh  size of pdf array: [nx ns nl nh] [20 21 19 18]
                          2  q.nh  number of probe SNR values to evaluate [30]
                          3  q.cs  weighting of pdf factors [1 1 1 1]
                          5  q.dh  minimum step size in dB for probe SNRs [0.2]
                          6  q.sl  min slopes threshold [0.02]
                          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]
              xp{n}(:) list of available probe SNR values
              i        model probed
              x        probe SNR value used
              r        response: 0=wrong, 1=right.

     Outputs:
              xx       recommended probe SNR
              ii       recommended model to probe next
              m(4,n,3) estimated srt and slope of all models
                       m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
              v(4,n)   estimated diagonal covariance matrix entries:

## CROSS-REFERENCE INFORMATION

This function calls:

- v_quadpeak V_PEAK2DQUAD find quadratically-interpolated peak in a N-D array

This function is called by:

- v_psycdigit V_PSYCDIGIT measures psychometric function using TIDIGITS stimuli

## SUBFUNCTIONS

- function y=linres(x,d,a,b)

## SOURCE CODE

    0001 function [xx,ii,m,v]=v_psycestu(iq,x,r,xp)
    0002 %V_PSYCESTU estimate unimodal psychometric function
    0003 %
    0004 % Usage: [xx,ii,m,v]=v_psycestu(-n,p,q,xp) % initialize n models
    0005 %        [xx,ii,m,v]=v_psycestu(i,x,r)     % supply a trial result to v_psycest
    0006 %                    v_psycestu(i)         % plot pdf of model i
    0007 %              [p,q]=psychestu(0)        % output model parameters (or print them if no outputs)
    0008 %
    0009 % Inputs:
    0010 %         -n        minus the number of models
    0011 %          p(:,n)   parameters for each model
    0012 %                      1  miss [0.04]
    0013 %                      2  guess [0.1]
    0014 %                      3,4  SNR at peak: [min max] [-20 20]
    0015 %                      5,6  normalized semi-width: [min max] [0 20]
    0016 %                      7,8  low side slope: [min max] [0.03 0.3]
    0017 %                      9,10  high side slope: [min max] [0.03 0.3]
    0018 %          q(:)     parameters common to all models (vector or struct)
    0019 %                      1  q.nxslh  size of pdf array: [nx ns nl nh] [20 21 19 18]
    0020 %                      2  q.nh  number of probe SNR values to evaluate [30]
    0021 %                      3  q.cs  weighting of pdf factors [1 1 1 1]
    0022 %                      5  q.dh  minimum step size in dB for probe SNRs [0.2]
    0023 %                      6  q.sl  min slopes threshold [0.02]
    0024 %                      7  q.kp  number of std deviations of the pdfs to keep [4]
    0025 %                      8  q.hg  amount to grow expected gains in ni trials [1.3]
    0026 %          xp{n}(:) list of available probe SNR values
    0027 %          i        model probed
    0028 %          x        probe SNR value used
    0029 %          r        response: 0=wrong, 1=right.
    0030 %
    0031 % Outputs:
    0032 %          xx       recommended probe SNR
    0033 %          ii       recommended model to probe next
    0034 %          m(4,n,3) estimated srt and slope of all models
    0035 %                   m(:,:,1:3) are respectively the mean, mode (MAP) and marginal mode estimates
    0036 %          v(4,n)   estimated diagonal covariance matrix entries:
    0037 persistent pp qq fl fh fxs tx ts tl th
    0038 if iq<0
    0039     if iq~=-1
    0040         error('Cannot yet have multiple models');
    0041     end
    0042     pp=x;
    0043     qq=r;
    0044     pp(7:2:9)=max(pp(7:2:9),qq.sl);  % enforce the minimum slope
    0045 end
    0046 nxslh=qq.nxslh;  % make local copies of some useful values
    0047 nx=nxslh(1);
    0048 ns=nxslh(2);
    0049 nl=nxslh(3);
    0050 nh=nxslh(4);
    0051 la=pp(1);
    0052 gu=pp(2);
    0053 if iq<0
    0054     % reserve space for pdfs
    0055     fl=ones(nx,nl,ns);
    0056     fh=ones(nx,nh,ns);
    0057     fxs=ones(nx,1,ns); % marginal x-s distribution
    0058     % define initial ranges
    0059     tx=linspace(pp(3),pp(4),nx)';
    0060     ts=linspace(pp(5),pp(6),ns)';
    0061     tl=linspace(log(pp(7)),log(pp(8)),nl)';
    0062     th=linspace(log(pp(9)),log(pp(10)),nh)';
    0063 elseif iq>0 && nargin==3
    0064     if iq~=1
    0065         error('Cannot yet have multiple models');
    0066     end
    0067     r0=r==0;
    0068     % Update the pdf arrays
    0069     mskx=tx>x; % find which tx values are >x
    0070     nxm=sum(mskx);
    0071     ets=reshape(exp(-ts),[1 1 ns]);
    0072     sd0=1+ets.^2;
    0073     sm0=2*ets;
    0074     if any(mskx)
    0075         fl(mskx,:,:)=fl(mskx,:,:).*(r0+(1-2*r0)*(gu+(1-gu-la)*(repmat(sd0,[nxm,nl,1])+repmat(sm0,[nxm,nl,1]).*repmat(cosh((tx(mskx)-x)*exp(tl')),[1 1 ns])).^(-1)));
    0076     end
    0077     if ~all(mskx)
    0078         fh(~mskx,:,:)=fh(~mskx,:,:).*(r0+(1-2*r0)*(gu+(1-gu-la)*(repmat(sd0,[nx-nxm,nh,1])+repmat(sm0,[nx-nxm,nh,1]).*repmat(cosh((x-tx(~mskx))*exp(th')),[1 1 ns])).^(-1)));
    0079     end
    0080 else % iq=0
    0081     xx=pp;
    0082     ii=qq;
    0083 end
    0084
    0085 % Normalize the pdf arrays
    0086
    0087 sfl=sum(fl,2);  % sum over slope variable
    0088 fl=fl./sfl(:,ones(nl,1),:); % normalize each x-s plane
    0089 sfh=sum(fh,2);
    0090 fh=fh./sfh(:,ones(nh,1),:); % normalize each x-s plane
    0091 fxs=fxs.*sfl.*sfh;          % Marginal x-s distribution
    0092 fxs=fxs/sum(fxs(:));        % Normalize to 1
    0093
    0094 % calculate marginal distributions
    0095
    0096 px=squeeze(sum(fxs,3));
    0097 ps=squeeze(sum(fxs,1));
    0098 pl=reshape(permute(fl,[2 1 3]),nl,nx*ns)*fxs(:);
    0099 ph=reshape(permute(fh,[2 1 3]),nh,nx*ns)*fxs(:);
    0100
    0101 % calculate means, modes and entropies
    0102
    0103 m=zeros(4,1,3);
    0104 m(:,1,1)=[tx'*px ts'*ps tl'*pl th'*ph]';
    0105 vraw=[tx'.^2*px ts'.^2*ps tl'.^2*pl th'.^2*ph]'-m(:,1,1).^2;
    0106 % h=[(tx(1)-tx(2))*px'*log(px) (ts(1)-ts(2))*ps'*log(ps) (tl(1)-tl(2))*pl'*log(pl) (th(1)-th(2))*ph'*log(ph)]';
    0107 % calculate joint mode
    0108 [flm,iflm]=max(fl,[],2);
    0109 [fhm,ifhm]=max(fh,[],2);
    0110 fxsm=fxs.*flm.*fhm;
    0111 [pxpk,i]=max(fxsm(:)); % find global mode
    0112 j=1+floor((i-1)/nx);
    0113 i=i-nx*(j-1);
    0114 jl=iflm(i,1,j);
    0115 jh=ifhm(i,1,j);
    0116 % could use quadratic interpolation but it seems a bit flaky for 4-D
    0117 % if all([i j jl jh]>1) && all([i j jl jh]<[nx ns nl nh])
    0118 %     [dum,ij]=v_quadpeak(repmat(fxs(i-1:i+1,1,j-1:j+1),[1 3 1 3]).*repmat(fl(i-1:i+1,jl-1:jl+1,j-1:j+1),[1 1 1 3]).*permute(repmat(fh(i-1:i+1,jh-1:jh+1,j-1:j+1),[1 1 1 3]),[1 4 3 2]));
    0119 % end
    0120 m(:,1,2)=[tx(i) ts(j) tl(jl) th(jh)]'; % replicate mean for now
    0121 % calculate marginal modes
    0122 [pxpk,i]=max(px);                          % marginal mode in x
    0123 if i>1 && i<nx                           % use quadratic interpolation if possible
    0124     [dum,j]=v_quadpeak(px(i-1:i+1));
    0125     i=i+j-2;
    0126 end
    0127 xmm=(2-i)*tx(1)+(i-1)*tx(2);             % marginal mode in x
    0128
    0129 [pxpk,i]=max(ps);                          % marginal mode in sw
    0130 if i>1 && i<ns                           % use quadratic interpolation if possible
    0131     [dum,j]=v_quadpeak(ps(i-1:i+1));
    0132     i=i+j-2;
    0133 end
    0134 smm=(2-i)*ts(1)+(i-1)*ts(2);             % marginal mode in sw
    0135
    0136 [pxpk,i]=max(pl);                          % marginal mode in l
    0137 if i>1 && i<nl                           % use quadratic interpolation if possible
    0138     [dum,j]=v_quadpeak(pl(i-1:i+1));
    0139     i=i+j-2;
    0140 end
    0141 lmm=(2-i)*tl(1)+(i-1)*tl(2);             % marginal mode in l
    0142
    0143 [pxpk,i]=max(ph);                          % marginal mode in h
    0144 if i>1 && i<nh                           % use quadratic interpolation if possible
    0145     [dum,j]=v_quadpeak(ph(i-1:i+1));
    0146     i=i+j-2;
    0147 end
    0148 hmm=(2-i)*th(1)+(i-1)*th(2);             % marginal mode in h
    0149 m(:,1,3)=[xmm smm lmm hmm]';
    0150 m(3:4,:,:)=exp(m(3:4,:,:));              % convert log slopes to slopes
    0151 mfact=(m(3,1,:)+m(4,1,:))./(m(3,1,:).*m(4,1,:));
    0152 m(2,1,:)=m(2,1,:).*mfact;  % convert normalized semi-width to width
    0153 v=[vraw(1); vraw(2)*mfact(1).^2; vraw(3:4).*m(3:4,1,1).^2]; % calculate variances
    0154
    0155 % figure(21); plot(tx,px,m([1 1]),[0 1.05*max(px)]); title('SNR at peak');
    0156 % figure(22); plot(ts,ps,m([2 2]),[0 1.05*max(ps)]); title('semi-width');
    0157 % figure(23); plot(tl,pl,log(m([3 3])),[0 1.05*max(pl)]); title('lower log slope');
    0158 % figure(24); plot(th,ph,log(m([4 4])),[0 1.05*max(ph)]); title('upper log slope');
    0159 if ~nargout && iq==1  % plot pdf
    0160     subplot(121)
    0161     imagesc(tx,ts,squeeze(fxs)');
    0162     axis 'xy'
    0163     xlabel('Peak position (dB)');
    0164     ylabel('Normalized semi-width (dB)');
    0165     subplot(122)
    0166     imagesc(th,tl,reshape(permute(fl,[2 1 3]),nl,nx*ns)*reshape(permute(fh.*fxs(:,ones(nh,1),:),[1 3 2]),nx*ns,nh));
    0167         axis 'xy'
    0168         xlabel('Ln down slope (ln prob/dB)');
    0169          ylabel('Ln up slope (ln prob/dB)');
    0170 end
    0171 % now determine the next recommended probe SNR
    0172
    0173 if iq~=0
    0174 xt=tx;  % for now we just try all the tx values
    0175 nxt=numel(xt);
    0176 ets=reshape(exp(-ts),[1 1 ns]);
    0177 sd0=1+ets.^2;
    0178 sm0=2*ets;
    0179 hh=zeros(nxt,1);  % store for entropies
    0180 for i=1:nxt
    0181     y=xt(i);        % y = probe value
    0182     mskx=tx>y; % find which tx values are >y
    0183     nxm=sum(mskx);
    0184     flm=gu+(1-gu-la)*(repmat(sd0,[nxm,nl,1])+repmat(sm0,[nxm,nl,1]).*repmat(cosh((tx(mskx)-y)*exp(tl')),[1 1 ns])).^(-1);
    0185     fhm=gu+(1-gu-la)*(repmat(sd0,[nx-nxm,nh,1])+repmat(sm0,[nx-nxm,nh,1]).*repmat(cosh((y-tx(~mskx))*exp(th')),[1 1 ns])).^(-1);
    0186     % calculate probs conditional on r=1
    0187     fl1=fl;
    0188     fh1=fh;
    0189     fl1(mskx,:,:)=fl1(mskx,:,:).*flm;
    0190     fh1(~mskx,:,:)=fh1(~mskx,:,:).*fhm;
    0191     sfl1=sum(fl1,2);  % sum over slope variable
    0192     fl1=fl1./sfl1(:,ones(nl,1),:); % normalize each x-s plane
    0193     sfh1=sum(fh1,2);
    0194     fh1=fh1./sfh1(:,ones(nh,1),:); % normalize each x-s plane
    0195     fxs1=fxs.*sfl1.*sfh1;          % Marginal x-s distribution
    0196     pr=sum(fxs1(:)); % P(r=1 | y)
    0197     fxs1=fxs1/pr;        % Normalize to 1
    0198     px1=squeeze(sum(fxs1,3));
    0199     ps1=squeeze(sum(fxs1,1));
    0200     pl1=reshape(permute(fl1,[2 1 3]),nl,nx*ns)*fxs1(:);
    0201     ph1=reshape(permute(fh1,[2 1 3]),nh,nx*ns)*fxs1(:);
    0202     % calculate probs conditional on r=0
    0203     fl0=fl;
    0204     fh0=fh;
    0205     fl0(mskx,:,:)=fl0(mskx,:,:).*(1-flm);
    0206     fh0(~mskx,:,:)=fh0(~mskx,:,:).*(1-fhm);
    0207     sfl0=sum(fl0,2);  % sum over slope variable
    0208     fl0=fl0./sfl0(:,ones(nl,1),:); % normalize each x-s plane
    0209     sfh0=sum(fh0,2);
    0210     fh0=fh0./sfh0(:,ones(nh,1),:); % normalize each x-s plane
    0211     fxs0=fxs.*sfl0.*sfh0;          % Marginal x-s distribution
    0212     fxs0=fxs0/(1-pr);        % Normalize to 1
    0213     px0=squeeze(sum(fxs0,3));
    0214     ps0=squeeze(sum(fxs0,1));
    0215     pl0=reshape(permute(fl0,[2 1 3]),nl,nx*ns)*fxs0(:);
    0216     ph0=reshape(permute(fh0,[2 1 3]),nh,nx*ns)*fxs0(:);
    0217     % now calculate the entropies
    0218     h1=[(tx(1)-tx(2))*px1'*log(px1) (ts(1)-ts(2))*ps1'*log(ps1) (tl(1)-tl(2))*pl1'*log(pl1) (th(1)-th(2))*ph1'*log(ph1)]';
    0219     h0=[(tx(1)-tx(2))*px0'*log(px0) (ts(1)-ts(2))*ps0'*log(ps0) (tl(1)-tl(2))*pl0'*log(pl0) (th(1)-th(2))*ph0'*log(ph0)]';
    0220     hh(i)=qq.cs*(pr*h1+(1-pr)*h0);
    0221 end
    0222 [hmin,ih]=min(hh);
    0223 xx=xt(ih);
    0224 ii=1;
    0225 end
    0226
    0227 % now rescale the pdf arrays
    0228
    0229 % sraw=sqrt(vraw);  % std deviations
    0230 % clim=[tx(1) tx(end); ts(1) ts(end); tl(1) tl(end); th(1) th(end)]  % current axis limits
    0231 % minlim=clim*[3 1; 0 2]/3; % always keep at least the middle third of the existing array
    0232 % maxlim=clim*[2 0; 1 3]/3;
    0233 % pra=min(max(repmat(m(:,1,1),1,2)+qq.kp*sraw*[-1 1],minlim),maxlim);
    0234
    0235 plow=max(min(0.1*gu,0.001),1e-6);
    0236
    0237 pcx=cumsum(px);
    0238 tmp=linspace(tx(max(find(pcx>plow,1)-1,1)),tx(find(pcx>(1-plow),1)),nx)';
    0239 fxs=linres(fxs,1,tx,tmp);
    0240 fl=linres(fl,1,tx,tmp);
    0241 fh=linres(fh,1,tx,tmp);
    0242 tx=tmp;
    0243 %
    0244 pcx=cumsum(ps);
    0245 tmp=linspace(ts(max(find(pcx>plow,1)-1,1)),ts(find(pcx>(1-plow),1)),ns)';
    0246 fxs=linres(fxs,3,ts,tmp);
    0247 fl=linres(fl,3,ts,tmp);
    0248 fh=linres(fh,3,ts,tmp);
    0249 ts=tmp;
    0250 %
    0251 % For now, we do not update the slope distributions
    0252 %
    0253 % tmp=linspace(pra(3,1),pra(3,2),nl)';
    0254 % fl=linres(fl,2,tl,tmp);
    0255 % tl=tmp;
    0256 %
    0257 % tmp=linspace(pra(4,1),pra(4,2),nh)';
    0258 % fh=linres(fh,2,th,tmp);
    0259 % th=tmp;
    0260
    0261 function y=linres(x,d,a,b)
    0262 % linear resample x along dimension d from grid a to b
    0263 sz=size(x);
    0264 n=sz(d);
    0265 p=[d:numel(sz) 1:d-1];
    0266 sz2=prod(sz)/n;
    0267 z=reshape(permute(x,p),n,sz2);  % put the target dimension in column 1
    0268 na=numel(a);
    0269 nb=numel(b);
    0270 [xx,ix]=sort([a(:); b(:)]);
    0271 jx=zeros(1,na+nb);
    0272 jx(ix)=(1:na+nb);
    0273 % ja=jx(1:na)-(1:na); % last new point < each old point
    0274 jb=jx(na+1:end)-(1:nb); % last old point < each new point
    0275 y=zeros(nb,size(z,2));
    0276 y(jb==0,:)=z(ones(sum(jb==0),1),:);  % replicated entries
    0277 y(jb==na,:)=z(na*ones(sum(jb==na),1),:);
    0278 mk=(jb>0) & (jb<na);  % interpolation mask
    0279 jmk=jb(mk);
    0280 f=((b(mk)-a(jmk))./(a(1+jmk)-a(jmk)));
    0281 fc=1-f;
    0282 y(mk,:)=z(jmk,:).*fc(:,ones(1,sz2))+z(1+jmk,:).*f(:,ones(1,sz2));
    0283 sz(d)=nb;
    0284 y=ipermute(reshape(y,sz(p)),p);
    0285
    0286
    0287
    0288
    0289
    0290