Chapter 5: complex observation functions in threshold detection

In this problem, more so than in many other examples so far, there are multiple parts to the solution to the measurement problem.

First, let us assume that we know the simulated lapse rate and shape parameter (lambda=0.04 and sigma=4). If these values are known, then the solution for the first dataset is, at base, just computing measurements for binomial rate constants; but whereas simply making 11 such measurements would not provide any conceptual difficulty, here the 11 rates are related by the underlying value of the sigmoid. If we were to draw a sigmoid with a median of 1 atop the circular data, the sigmoid would have an even sharper initial upturn and the rate constants predicted by this pattern would not match the 11 observed rates as well as if the median were closer to 2.5.

There are several steps to creating the plots in Fig. 5.11b. First, we define a few constants:

clear; load Sdat t1 d1

tlist=reshape(t1,101,11); tlist=tlist(1,:); Nlist=sum(reshape(d1,101,11));

MT=linspace(0,10,302)'; MT=MT(2:end); Omt=ones(size(MT)); Ot=ones(size(Nlist)); ONL=Omt*Nlist;

We start the solution with the binomial likelihood, comparing data to a candidate rate constant:

L=@(Rin) ONL.*log(Rin)+(101-ONL).*log(1-Rin);

This equation requires the lapse-adjusted prediction of correct responding, s:

R=@(Din,Lin) (1-Lin)*(.5*Din+.5)+Lin*.5;

and this equation in turn requires the definition of the selected ogive :

D=@(Sin) 1./(1+((Omt*tlist)./(MT*Ot)).^-Sin);

We combine these in a loop that tests possible lapse rates (lambda), shape parameters (SP), and medians (MT):

SP=linspace(0,10,102); SP=SP(2:end); lambda=linspace(0,.1,101); Lr=nan(301,101,21); %SP,MT,lambda

for nlambda=1:21,

for nSP=1:101, Rnow=R(D(SP(nSP)),lambda(nlambda));

Lr(:,nSP,nlambda)=sum(L(Rnow),2)-SP(nSP); end, end

LSr=logsum(Lr,3,diff(lambda(1:2)));

figure(7); mesh(SP,MT,exp(LSr-max(max(LSr)))); colormap bone

LSr=logsum(Lr,2:3,[diff(SP(1:2)) diff(lambda(1:2))]); LSr=LSr-max(LSr);

figure; plot(MT,LSr,'k.'); axis([min(MT) max(MT) min(LSr) 0])

MThat=logpeakval([LSr MT]); ind=findnearestN(MT,2.5,1);

hold on; plot(2.5*[1 1],[min(LSr) LSr(ind)],'k--','LineWidth',2)

ind=findnearestN(MT,MThat,1); plot(MThat,min(LSr),'ko','MarkerFaceColor','k','MarkerSize',8)

The plotting at the end reproduces the measurement for the first dataset. The measurement for the second dataset is left as an exercise. Now use a variant of the D function above to plot the sigmoid corresponding to the change in detection as stimulus strength increases:

t=linspace(0,10,501); t=t(2:end); lambda=0.04; sigma=4;

F=1./(1+(t./MThat).^-sigma);

s=(1-lambda)*(.5*F+.5)+lambda*.5;

figure; plot(t,s,'k.'); hold on;

plot(tlist,Nlist/100,'ko','MarkerFaceColor',.6*[1 1 1],'MarkerSize',13)

which requires that you choose a shape parameter, sigma. Test several and see how the result changes.