Chapter 5: measuring exponentially decaying sensory-motor error

%%% Sensory-Motor %%%

%%%%%%%%%%%%%%%%

%%% Example 1 %%

%%% Create DATA

T=26; nTau=751; nSIG=201; t1=[-T/2:-1]'; t2=[0:T]';

A=100; tau=[4]; SD=7;

D1=SD*randn([T/2 1]);

D2pre=SD*randn([T+1 1]); D2pre=D2pre-mean(D2pre);

dprime=A*exp(-t2./tau);

D2=dprime+D2pre;

figure(1); clf;

subplot(2,1,1); hold on;

plot([0 T],A*[1 1],'k--','LineWidth',2)

plot([0 T],[0 0],'k--','LineWidth',2);

plot([t1; t2],[D1; D2],'ko','MarkerFaceColor',.4*[1 1 1],'MarkerSize',10,'LineWidth',2.25);

plot(t1,zeros(1,T/2),'--','Color',.7*[1 1 1],'LineWidth',1.2);

plot(t2,dprime,'--','Color',.7*[1 1 1],'LineWidth',1.2);

axis([-15 30 -10 110]); xlabel('reach number (pert relative)'); ylabel('reach error (mm)')

Tlist=linspace(1,10,nTau+1); Tlist=Tlist(2:end); L=zeros(nTau,nSIG);

Slist=linspace(0,25,nSIG+1); Slist=Slist(2:end);

fD2=A*[exp(-[t2*ones(1,nTau)]./[ones(T+1,1)*Tlist])];

D1now=D1; D2now=D2*ones(1,nTau);

for nSig=1:nSIG, Snow=Slist(nSig);

L(:,nSig)=-(T+1)*log(Snow)-.5*(sum((D1now).^2)./Snow^2)-.5*(sum((fD2-D2now).^2)./Snow^2); end

LT=logsum(L-max(max(L)),2,diff(Slist(1:2)));

%tau estimate

subplot(2,1,2); hold on

plot(Tlist,exp(LT-max(LT)),'--','Color',[.2 .4 .7],'LineWidth',2); axis([1 10 0 1.01])

xlabel('time constant (tau)'); ylabel('p')

LPVt=logpeakval([LT(:)-max(LT) Tlist(:)]);

%Convert Tau to Gamma

Ntrials=20;

gamma=linspace(.1,.6,501); tau=nan(size(gamma));

elist=1./linspace(1,10,1501); elist=elist(end:-1:1);

es=exp(-[0:Ntrials]'*elist);

for n=1:length(gamma),

err=ones(Ntrials+1,1); plot(0,err,'k.')

for nn=2:(Ntrials+1),

err(nn)=err(nn-1)-err(nn-1)*gamma(n); end

%plot([0:Ntrials]',err,'ko');

%plot(0:Ntrials,es(:,i),'--')

L=sum((err(:)*ones(1,length(elist))-es).^2,1);

i=find(isnear(L,min(L)));

tau(n)=1/mean(elist(i)); end

figure(2); clf; subplot(2,1,2); hold on

plot(tau,gamma,'k.'); xlabel('time constant (tau)'); ylabel('update proportion (gamma)')

%%%%%%%%%%%%%%%%

%%% Example 2 %%

%%% Create DATA

T=26; nTau=91; nAlpha=351; nBeta=41; nSIG=101; t1=[-T/2:-1]'; t2=[0:T]';

A=100; alpha=[.75]; beta=[-.05]*A; tau=[1.75]; SD=7;

D1=beta+(SD*randn([T/2 1]));

D2pre=SD*randn([T+1 1]); D2pre=D2pre-mean(D2pre);

dprime=A*alpha.*exp(-t2./tau)+[A.*(1-alpha)]+beta;

D2=dprime+D2pre;

figure(3); clf;

subplot(2,1,1); hold on;

plot([0 T],A*[1 1],'k--','LineWidth',2)

plot([0 T],[0 0],'k--','LineWidth',2);

plot([t1; t2],[D1; D2],'ko','MarkerFaceColor',.4*[1 1 1],'MarkerSize',10,'LineWidth',2.25);

plot(t1,beta*ones(1,T/2),'--','Color',.7*[1 1 1],'LineWidth',1.2);

plot(t2,dprime,'--','Color',.7*[1 1 1],'LineWidth',1.2);

Tlist=linspace(.3,3,nTau+1); Tlist=Tlist(2:end); L=zeros(nAlpha,nBeta,nTau,nSIG);

Alist=linspace(.25,1.25,nAlpha)'; Slist=linspace(0,25,nSIG+1); Slist=Slist(2:end);

Blist=linspace(-.25*A,.25*A,nBeta)';

fD1=beta;

fD2=@(Anow,Bnow) Bnow+A*[1-Anow*(1-exp(-[t2*ones(1,nTau)]./[ones(T+1,1)*Tlist]))];

D1now=D1; D2now=D2*ones(1,nTau);

for nA=1:nAlpha, Anow=Alist(nA);

for nB=1:nBeta, Bnow=Blist(nB);

for nSig=1:nSIG, Snow=Slist(nSig);

L(nA,nB,:,nSig)=-(T+1)*log(Snow)-.5*(sum((fD1-D1now).^2)./Snow^2)-.5*(sum((fD2(Anow,Bnow)-D2now).^2)./Snow^2);

end, end, end

L2=logsum(L-max(max(max(max(L)))),[2 4],[diff(Blist(1:2)) diff(Slist(1:2))]);

%amplitude estimate

subplot(2,1,2); hold on

LA=logsum(L2-max(max(L2)),[2],[diff(Tlist(1:2))]);

plot(Alist,exp(LA-max(LA)),'--','Color',[.2 .4 .7],'LineWidth',2); %axis([LPV+[-.1 .1] 0 1.01])

LPVa=logpeakval([LA-max(LA) Alist]);

axis([.1 1 0 1.01])

%tau estimate

figure(4); clf

subplot(2,1,2); hold on

LT=logsum(L2-max(max(L2)),[1],[diff(Alist(1:2))]);

plot(Tlist,exp(LT-max(LT)),'--','Color',[.2 .4 .7],'LineWidth',2); axis([.3 3 0 1.01])

LPVt=logpeakval([LT(:)-max(LT) Tlist(:)]);