%generate predictions from diffusion model, by simulation and by analytic approximation
%all times are in seconds

%parameters for diffusion model
z = .05;
mu = .15;
a = .1;
sigma = .1;

tau = .001; %timestep for approximating continuous diffusion process
maxRT = 1; %how far to let each trial run
ntraj = 30; %number of sample trajectories
r = zeros(ntraj,1); %set of simulated responses
RT = zeros(ntraj,1); %set of simulated RTs

%plot sample trials
figure(1),clf,hold on
RTmax = 0; %for tracking maximum of simulated RTs
for i=1:ntraj
    E = z + (tau:tau:maxRT)*mu + cumsum(randn(1,maxRT/tau))*sigma*sqrt(tau); %response process for this trial
    T = find(E<=0|E>=a,1); %index of first-passage time
    r(i) = 1 + (E(T)<=0); %response for this trial; 1 for upper threshold, 2 for lower
    RT(i) = T*tau; %RT for this trial
    plot(tau*(1:T),E(1:T),'Color',ones(1,3)*rand(1)) %plot trajectory for this trial in random greyscale
end
plot([0 maxRT],[a a],'k') %upper threshold (response 1)
plot([0 maxRT],[0 0],'k') %lower threshold (response 2)

%get predictions analytically, using the Fourier-type approximation
g = zeros(2,maxRT/tau); %joint density of response and RT
threshold = 1e-40; %for each t, sum over k continues until guaranteed |summand|/k < threshold
for tindex=1:maxRT/tau
    t = tindex*tau;
    A = pi*sigma^2/a^2;
    B = exp((2*mu*[a-z;-z]-t*mu^2)/(2*sigma^2)); %one value for response 1 (a-z), one for response 2 (-z)
    C = -pi^2*sigma^2*t/(2*a^2);
    D = pi*[a-z;z]/a; %one value for response 1 (a-z), one for response 2 (z)
    maxk = sqrt(1/C*log(threshold)); %exp(maxkn^2*C) = threshold
    temp = 0;
    for k=1:maxk
        temp = temp + k*sin(k*D)*exp(k^2*C);
    end
    g(:,tindex) = temp.*(A*B);
end

%plot analytic predictions
densityScale = a/max(max(g)); %for vertical scale of RT densities; greatest peak will have height equal to a
plot(tau:tau:maxRT,a+g(1,:)*densityScale,'k') %joint RT density for upper response
plot(tau:tau:maxRT,-g(2,:)*densityScale,'k') %joint RT density for lower response
box
axis([0 t -a 2*a])
set(gca,'YTick',[])
set(gca,'FontSize',22) %changes tick labels and any future axis labels
xlabel('Reponse Time')
set(gca,'FontSize',16) %changes tick labels but not existing axis labels
