# Stochastic Differential Equations -- basic question

**URL:** <https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200>\
**Category:** General Usage\
**Created:** [August 3, 2020, 3:31pm UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200 "2020-08-03T15:31:07Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 3, 2020, 3:31pm UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/1 "2020-08-03T15:31:07Z")

</div>

I’m making my first attempt of _Stochastic Differential Equations_ using a SIR model, and follow the stochastic formulation described in Keeling and Rohani (2008): Modeling Infectious Diseases in Humans and Animals, Princeton University Press.

My questions (at the **bottom**!) relates to the description of the noise + the choice of solver. First, the deterministic model (extensive variables, with immigration and emigration – sorry for my partially “private” notation):  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/b/eba97bb97d246e3eb341a35c40a61c76fade23f7.png)  
When neglecting immigration and emigration and using data from Martcheva’s 2015 Springer book (An Introduction to Mathematical Epidemiology), the deterministic model gives:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/f/1f5fab25ef19ac2a03ed6dc18c8c438bf16985e9.png)

Assuming Poisson distribution for the events leading to infection and recovery, Keeling and Rohani’s description gives the Stochastic Differential Equation:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/7/974ee6278451fb9e23a631f8ce18262a76c57924.png)

Here, f() is the drift/“deterministic” term while G() is the diffusion/“stochastic” term, while dw holds two independent Wiener processes:

```julia
# Drift function for SIR model
function sir_drift!(f,x,p,t)
    S,I,R = x
    tau_i,tau_r,N = p
    #
    f[1] = -I*S/tau_i/N
    f[2] = I*S/tau_i/N - I/tau_r
    f[3] = I/tau_r
end
# Diffusion model for SIR model
function sir_diffuse!(G,x,p,t)
    S,I,R = x
    tau_i,tau_r,N = p
    #
    G[1,1] = -sqrt(I*S/tau_i/N)
    G[1,2] = 0.
    G[2,1] = sqrt(I*S/tau_i/N)
    G[2,2] = - sqrt(I/tau_r)
    G[3,1] = 0.
    G[3,2] = sqrt(I/tau_r)
end

```

Next, I define the problem; elements of the description:

```julia
#
prob = SDEProblem(sir_drift!,sir_diffuse!,x3,tspan,p,noise_rate_prototype=zeros(3,2))
ensembleprob = EnsembleProblem(prob)
sol = solve(ensembleprob,EnsembleThreads(),trajectories=100)

```

I can compute some statistics, and find:

![image](https://global.discourse-cdn.com/julialang/original/3X/d/d/dd134336c2302f5c73fcf4f87cd5faadb7852b45.png)

**Question 1** : I have specified the noise as `noise_rate_prototype=zeros(3,2)`. Is this correct for saying that there are two independent Wiener processes for the 3 states? Does it matter whether I use `zeros(3,2)` or `ones(3,2)` or something else, or is the important thing that the matrix has correct dimension? [And why not just say that I have a _vector_ of Wiener processes with 2 elements??]

**Question 2** : I have not specified the solver. I tried to specify solver `SRWI1`, but that gave an error message.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 3, 2020, 11:29pm UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/2 "2020-08-03T23:29:54Z")

</div>

> [@BLI](#):
>
> **Question 1** : I have specified the noise as `noise_rate_prototype=zeros(3,2)` . Is this correct for saying that there are two independent Wiener processes for the 3 states? Does it matter whether I use `zeros(3,2)` or `ones(3,2)` or something else, or is the important thing that the matrix has correct dimension? [And why not just say that I have a _vector_ of Wiener processes with 2 elements??]

It just needs to have the correct type. The reason why we want the matrix, the prototype, is that this allows structures like sparse, Tridiagonal, etc. for specializing the computation.

> [@BLI](#):
>
> **Question 2** : I have not specified the solver. I tried to specify solver `SRWI1` , but that gave an error message.

SRIW1 cannot handle this one. LambaEM is probably the right call here, and probably what it uses.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 4, 2020, 9:26am UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/3 "2020-08-04T09:26:26Z")

</div>

Thanks! If I understand you correctly, if I set `noise_rate_prototype = A` where `A = sparse([1 0;1 1;0 1])` (in my case), the solver can take advantage of the structure of G(). [In my simple case, there is probably little to win by that.]

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 4, 2020, 9:45am UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/5 "2020-08-04T09:45:19Z")

</div>

One more question on manual computing of statistics… it is “trivial” to compute the mean, the median, and the std, e.g.:

```julia
sol_a = Array(sol)
sol_mean = mean(sol_a; dims=3)[:,:,1] |> permutedims
sol_median = median(sol_a; dims=3)[:,:,1] |> permutedims
sol_std = std(sol_a; dims=3)[:,:,1] |> permutedims
plot(sol_mean)
plot!(sol_median)
plot!(sol_mean,ribbon=(sol_std,sol_std))

```

**However** – for some reason, function `quantile` does not seem to allow for the `dims` keyword.

_Question_: Does this mean that `quantile` only accepts _vectors_ as first argument, and that I need to loop through both the number of states (first dimension of `sol`, = 3 in my case) and the number of time points (second dimension of `sol`, = 111 in my case)? In other words, something like the following:

```julia
[quantile(sol_a[i,j,:],0.5) for j in 1:111, i in 1:3]

```

in my case of 3 states, 111 time points, and the median (0.5 quantile)? Or is there a built-in method that handles `dims`?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 4, 2020, 10:46am UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/6 "2020-08-04T10:46:17Z")

</div>

Ensembles and reliability in statistics… I seem to recall that Bradley Efron recommends the following in his book on Bootstrapping [or maybe I read it somewhere else…]:

- For reliable computation of mean and variance, 50-100 trajectories are needed
- For reliable computation of quantiles, etc., 1000+ trajectories are needed

Are these good rules of thumb in general for statistics based on an ensemble of realizations?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 4, 2020, 11:21am UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/7 "2020-08-04T11:21:38Z")

</div>

> [@BLI](#):
>
> _Question_ : Does this mean that `quantile` only accepts _vectors_ as first argument, and that I need to loop through both the number of states (first dimension of `sol` , = 3 in my case) and the number of time points (second dimension of `sol` , = 111 in my case)? In other words, something like the following:

Try using `Array(sol_a)`?

> [@BLI](#):
>
> Are these good rules of thumb in general for statistics based on an ensemble of realizations?

Depends on the amount of nonlinearity and the `dt`. But as a rule of thumb, 10,000 trajectories is generally enough to saturate convergence for mean effects, and you need a lot higher for variance and quantiles to be very accurate from what I’ve seen.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [August 4, 2020, 1:21pm UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/8 "2020-08-04T13:21:04Z")

</div>

With `sol` being the solution of the Ensemble problem:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/c/cc58350180b750e6fd0d046454380d08275cd544.png)

It works exactly the same way if I replace `mean` with `median`, `std`, etc.

But for `quantile` (where I need an extra parameter), the following does _not_ work (with, e.g., the `0.3` quantile):

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/9/39b26a99069da7c03192cf7b74098bc39233b9e2.png)

I wrote my own function for doing what I want, but it is not perfect – it only works for 2 and 3 dim matrices (I think), and there is no proper testing of arguments.

```julia
function quantile_my(A,q=0.5;dims=1)
    sA = size(A)
    nA = length(sA)
    iA = setdiff(1:nA,dims)
    if nA == 2
        return permutedims([quantile(A[i,:],q) for i in 1:sA[iA[1]]])
    elseif nA == 3
        return permutedims([quantile(A[i,j,:],q) for j in 1:sA[iA[2]], i in 1:sA[iA[1]]])
    else
        println("Doesn't work for your array dimension")
    end
end

```

It would be nice if `quantile` were “symmetric” in the sense of offering the same input array as the other statistical functions. Maybe I’m just overlooking the documentation.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 4, 2020, 1:26pm UTC](https://discourse.julialang.org/t/stochastic-differential-equations-basic-question/44200/9 "2020-08-04T13:26:20Z")

</div>

That seems like an issue in Julia’s Base then.