Home > voicebox > psycestu.m

psycestu

PURPOSE ^

psycestu estimate unimodal psychometric function

SYNOPSIS ^

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

DESCRIPTION ^

 psycestu estimate unimodal psychometric function

 Usage: [xx,ii,m,v]=psycestu(-n,p,q,xp) % initialize n models
        [xx,ii,m,v]=psycestu(i,x,r)     % supply a trial result to psycest
                    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: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [xx,ii,m,v]=psycestu(iq,x,r,xp)
0002 % psycestu estimate unimodal psychometric function
0003 %
0004 % Usage: [xx,ii,m,v]=psycestu(-n,p,q,xp) % initialize n models
0005 %        [xx,ii,m,v]=psycestu(i,x,r)     % supply a trial result to psycest
0006 %                    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]=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]=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]=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]=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]=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

Generated on Fri 22-Sep-2017 19:37:38 by m2html © 2003