added two scripts to paper directory for simulated vs. CN and MC and added some text

Author: briandepasquale

paper/paper.md

@@ -44,12 +44,20 @@ bibliography: paper.bib
 
 # Summary
 
-Drift diffusion models (DDMs) are a popular model class for modeling a unobserved process that determines an subject's choice during a decision-making task [@Bogacz2006]. 
+Drift diffusion models (DDMs) are a popular model class for modeling a unobserved process that determines an subject's choice during a decision-making task [@Bogacz2006]. Mathematically, in their simplest form, they are equivalent of an [Ornstein-Uhlenbeck](https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process) process, a type of mean-reverting stochastic process similar to Brownian motion, described by the following stochastic differential equation
 
 ```math
- dz = \lambda zdt + u(t)dt + \sigma dW
+ dz = \lambda zdt + u(t)dt + \sigma dW \tag{1}
 ```
 
+These dynamics can be equivalently expressed as a partial differential equation, which describes the motion of the probability distribution of $z$
+
+```math
+\frac{\partial P(z(t))}{\partial t} = \frac{\sigma}{2}\frac{\partial^2 P}{\partial z^2} - \frac{\partial(\lambda zP)}{\partial z} - \frac{\partial(u(t)P)}{\partial z}. \tag{2}
+```
+
+In neuroscience, as mass of this distribution moves, this can be considered to be a simple model of the internal process by which evidence is accumlated and weighed between options. In cases where evidence is received continuously in time, an external input $u(t)$ is included. This external input forces our PDE to be solved numerically. 
+
 [@Brunton2013]. 
 
 # Statement of need

paper/test.jl

@@ -0,0 +1,32 @@
+using LinearAlgebra, PyPlot
+
+# Parameters
+dt = 1e-2
+B = 10.
+λ = -4.
+σ2 = 10.
+n = 53
+
+xc, dx = PulseInputDDM.bins(B, n);
+P0 = PulseInputDDM.P0(dx^2, n, dx, xc, dt) 
+M = PulseInputDDM.transition_M(σ2*dt, λ, 0., dx, xc, n, dt);
+
+# Calculate coefficients
+cDiff = σ2 * dt / (2*dx^2)
+cDrft = λ * dt / (2 * dx)
+Ddff = diagm(0 => ones(n) * -2cDiff, 1 => ones(n-1) * cDiff, -1 => ones(n-1) * cDiff)
+Ddff[:, [1, n]] .= 0
+Dder = diagm(0 => zeros(n), 1 => ones(n-1) * cDrft, -1 => ones(n-1) * -cDrft)
+Dder[:, [1, n]] .= 0
+C = exp(-Dder * diagm(0 => xc) + Ddff)
+
+#plot(xc, P0, label="P0")
+#plot(xc, M^(1/dt) * P0, label="M")
+#plot(xc, C^(1/dt) * P0, label="C")
+#legend()
+
+plot(sort(real(eigvals(M))), label="realM")
+plot(sort(abs.(imag(eigvals(M)))), label="imagM")
+plot(sort(real(eigvals(C))),  label="realC")
+plot(sort(abs.(imag(eigvals(C)))), label="imagC")
+legend()

paper/test2.jl

@@ -0,0 +1,114 @@
+using LinearAlgebra, SparseArrays, Random, Distributions, Plots, StatsBase
+
+# Clear all variables
+T = 1.0          # time of simulation
+dt = 2e-2        # Fokker-Planck approximations timestep
+nMC = 1e4        # number of Monte Carlo particles
+dtMC = 1e-4      # Monte Carlo timestep
+
+B = 1.0          # bound height
+vara = 1e-2      # diffusion variance
+lambda = -1e-3   # drift
+mu0 = 0.0 * B    # initial mean of distribution
+n = 53           # number of bins
+
+dx = 2 * B / (n - 2)   # bin width
+dx2 = dx^2             # dx^2
+vari = dx2             # make initial distribution as large as dx^2
+
+# Finite difference
+
+xe = vcat(-(B+dx), range(-B, -dx/2, ceil(Int, (n-3)/2+1))..., 
+          range(dx/2, B, floor(Int, (n-3)/2+1))..., B+dx) # edges of the bins
+xc = (xe[1:n] .+ dx/2)'  # bin centers
+
+cDiff = vara / (dx^2 * 2)  # scale factor for diffusion
+cDrft = lambda / (dx * 2)  # scale factor for drift
+
+# Diffusion matrix
+Ddff = diagm(0 => -2 * cDiff * [0; ones(n-2); 0]) +
+       diagm(-1 => cDiff * [0; ones(n-2)]) +
+       diagm(1 => cDiff * [ones(n-2); 0])
+
+# Drift matrix
+Dder = diagm(-1 => -cDrft * [0; ones(n-2)]) +
+       diagm(1 => cDrft * [ones(n-2); 0])
+
+# Multiply by a values and add
+DD = -Dder * diagm(0 => xc') + Ddff
+
+# Matrix exponential
+M = exp(DD * dt)
+
+# Brunton method
+
+# You would replace this with your implementation of `make_F` in Julia
+M_B = PulseInputDDM.transition_M(vara*dt, lambda, 0., dx, xc', n, dt);
+
+# Initialize
+
+Pa = pdf.(Normal(mu0, sqrt(vari)), xc) .* dx
+Pa /= sum(Pa)  # Fin. Diff.
+Pa = collect(Pa')
+Pa_B = copy(Pa)  # Brunton
+
+# For saving
+PA = zeros(n, round(Int, T/dt))
+PA_B = similar(PA)
+
+# Propagate
+
+for t in 1:round(Int, T/dt)
+    global Pa, Pa_B  # Ensure these refer to the global variables
+    PA[:, t] = Pa
+    PA_B[:, t] = Pa_B
+    
+    Pa = M * Pa  # Fin. Diff.
+    Pa_B = M_B * Pa_B  # Brunton (Uncomment this once you define `M_B`)
+    
+end
+
+# Monte Carlo
+
+# For saving
+PMC = fill(NaN, size(PA))
+
+a = mu0 .+ sqrt(vari) .* randn(Int(nMC))  # Initialize Gaussian
+
+for t in 1:round(Int, T/dtMC)
+
+    global a
+
+    a[a .< -(B+dx)] .= -(B+dx/2)
+    a[a .> B+dx] .= B+dx/2
+    
+    # Bin the individual particles
+    if mod(t-1, round(Int, dt/dtMC)) + 1 == 1
+        # Bin the particles and normalize by the number of Monte Carlo particles
+        counts = StatsBase.fit(Histogram, a, xe).weights
+        PMC[:, ceil(Int, dtMC*(t)/dt)] = (1/nMC) * counts
+        #PMC[:, ceil(Int, dtMC*(t)/dt)] = (1/nMC) * histcounts(a, xe)
+    end
+    
+    go = (a .< B) .& (a .> -B)  # Only integrate those that haven't crossed the boundary
+    
+    # OU process: 2 terms, 1-drift and 2-diffusion
+    a[go] .= a[go] .+ (lambda * dtMC) .* a[go] .+ sqrt(vara * dtMC) .* randn(sum(go))
+end
+
+# Plot over time
+
+p = plot()
+plot!(p, xc', PA[:, 1], label="Fin.Diff.", color="red", lw=2, marker=:o)
+plot!(p, xc', PA_B[:, 1], label="Brunton", color="green", lw=2, marker=:star)
+plot!(p, xc', PMC[:, 1], label="Monte Carlo", color="blue", lw=2, marker=:star)
+display(p)
+
+for i in 1:size(PA, 2)
+    p = plot()
+    plot!(p, xc', PA[:, i], label="Fin.Diff.")
+    plot!(p, xc', PA_B[:, i], label="Brunton")
+    plot!(p, xc', PMC[:, i], label="Monte Carlo")
+    sleep(0.1)
+    display(p)
+end