%Lecture_052312.m

%model parameters
p.u  = .1;  %drift rate
p.s =  1;   %s.d. of  drift rate
p.a = 2;  %units
p.b = 2; %units

dt = .01;  %step size in time (sampling rate, seconds)
maxT = 1; %seconds

t = (dt:dt:maxT); %time vector
nt = length(t);

%% 1. continuous step sizes

% generate a bunch of random walks

nWalks = 1000;

dy = p.u*dt +   p.s*sqrt(dt)*randn(nt,nWalks);

y = cumsum(dy);

figure(1)
clf
hold on

plot([0,maxT],p.a*[1,1],'k-');

plot([0,maxT],-p.b*[1,1],'k-');

plot([0,maxT],[0,0],'k:');


plot(t,y)

% histogram of last values of y

figure(2)
clf
hist(y(end,:))
title(sprintf('mean: %5.2f std. %5.2f',mean(y(end,:)),std(y(end,:))));

%% fixed y step size

yStep = p.s*sqrt(dt);

prob = .5*(p.u*sqrt(dt)/p.s+1);

dy = yStep*sign( rand(nt,nWalks)+prob-1);



y = cumsum(dy);

figure(3)
clf
hold on

plot([0,maxT],p.a*[1,1],'k-');

plot([0,maxT],-p.b*[1,1],'k-');

plot([0,maxT],[0,0],'k:');


stairs(t,y)

% histogram of last values of y

figure(4)
clf
hist(y(end,:))
title(sprintf('mean: %5.2f std. %5.2f',mean(y(end,:)),std(y(end,:))));


%% implement random walk model with boundaries

alive = true(1,nWalks);
RT = zeros(1,nWalks);
response = zeros(1,nWalks);

y = zeros(1,nWalks);
t =0;
while(sum(alive))
    dy = yStep*sign( rand(1,sum(alive))+prob-1);
    y(alive) = y(alive)+dy;
    t=t+dt;
    correct = y>=p.a;
    if sum(correct)  %deal with correct responses
        RT(correct) = t;
        response(correct) = 1;
        y(correct) = NaN;
        alive(correct) =false;
    end
    
        incorrect = y<=-p.b;
    if sum(incorrect)  %deal with correct responses
        RT(incorrect) = t;
        response(incorrect) = -1;
        y(incorrect) = NaN;
        alive(incorrect) =false;
    end
    
    
end

%%

%percent correct

bins = .1:.5:10

disp(sprintf('Percent correct: %5.2f',sum(response==1)/nWalks));

figure(5)
clf
subplot(1,2,1)
hist(RT(response==1),bins);
title('Correct responses');
set(gca,'XLim',[0,10]);
subplot(1,2,2)
hist(RT(response==-1),bins);
title('Incorrect responses');

set(gca,'XLim',[0,10]);

    
%expected percent correct

 ePC = expectedPC(p);
disp(sprintf('Expected percent correct: %5.2f',ePC));


%mean RTs


[eCorrect,eIncorrect] = expectedMeanRT(p);
fprintf('Correct meanRT: %5.4f, predicted: %5.4f\n',mean(RT(response==1)),eCorrect);
fprintf('Incorrect meanRT: %5.4f, predicted: %5.4f',mean(RT(response==-1)),eIncorrect);

%pdfs of RTs

[pdfCorrect,pdfIncorrect] = expectedRTpdf(p,bins);

pdfCorrect = pdfCorrect*sum(response==1)*(bins(2)-bins(1));
pdfIncorrect = pdfIncorrect*sum(response==-1)*(bins(2)-bins(1));

figure(5)
subplot(1,2,1)
hold on
plot(bins,pdfCorrect,'r-','LineWidth',2);
subplot(1,2,2)
hold on
plot(bins,pdfIncorrect,'r-','LineWidth',2);





