Chapter 5: dark adaptation

This problem would be essentially impossible using the ‘brute-force’ method of computing probability distributions we have so far been using, if it were not for the massive computational simplification afforded us by our prior information (as encoded by the indicator function, .

The first step is to load the dataset. There are two datasets. In the first, experimental light intensities used to measure threshold at each time-duration are chosen between the minimum and maximum thresholds as defined by our background information. In the second dataset, light intensities are chosen adaptively (so that if a correct response is given for one intensity, the next trial will try a lower intensity, and vice versa). You should try both.

load darkdatU %uniform light intensities

load darkdatA %adaptively chosen light intensities

The next step is to define the relevant functions. We will need a function to describe the exponentially decaying thresholds (f), a function to compute the correct-responding rates for a given underlying threshold (R), and a function to compute the log-posterior (L) associated with each possible tau value.

f=@(taunow) 5*exp(-Dt/taunow)+1;

R=@(Ivec,Ihat,Sin) .5./(1+(Ivec./Ihat).^-Sin)+.5;

L=@(Rin,Din) sum(Din.*log(Rin)+(1-Din).*log(1-Rin));

The next step is to define the list of tau and shape parameters to examine, and compute the log-probability of the result.

SP=linspace(0,10,102); SP=SP(2:end);

tau=linspace(1,20,303);

LT=nan(length(tau),length(SP),6);

for nT=1:length(tau), Iguess=f(tau(nT));

for n=1:6, for nSP=1:length(SP), SPnow=SP(nSP);

Rnow=R(Idat(:,n),Iguess(n),SPnow);

LT(nT,nSP,n)=L(Rnow,dat(:,n)); end, end, end

Notice that the loop above defines tau and the 6-element Iguess vector before cycling through possible shape parameters for each time-duration. This is part of the simplification that the indicator function affords us. If we had not used this simplification, we would have had to examine a range of possible thresholds for each time-duration independently, and there would have been six for-loops, one for each of the possible thresholds at each time-duration. Finally, the combination of these would be compared to the tau values predicted by each possible combination of thresholds.

To plot the marginalized log-posterior, we first take the product of probabilities across the 6 time-durations (sum of the log of the last dimension of the LT matrix), and then marginalize over the unknown shape parameter values (logsum of the second dimension).

figure; LS=logsum(sum(LT,3),2,[diff(SP(1:2))]); plot(tau,LS-max(LS),'k.'); box off

figure; hold on; box off

tauhat=logpeakval([LS(:) tau(:)]);

plot([0 30],6*[1 1],'k--','LineWidth',1.8)

plot([0 30],1*[1 1],'k--','LineWidth',1.8)

plot(dt,5*exp(-dt/tauhat)+1,'k.','LineWidth',2);

plot(Dt,Iprime,'ko','MarkerFaceColor',.6*[1 1 1],'MarkerSize',15,'LineWidth',2);

axis([-.05 30.05 -.05 6.05])