mschauer / Bridge.jl

A statistical toolbox for diffusion processes and stochastic differential equations. Named after the Brownian Bridge.
Other
111 stars 19 forks source link

Feedback and Contributing #12

Open mschauer opened 7 years ago

mschauer commented 7 years ago

This issue is held open for general feedback, feature requests and to coordinate contributions to the package; or just to say "hello".

getzdan commented 7 years ago

"hello" (sorry couldn't help it)

sdwfrost commented 3 years ago

Hi @mschauer

I'm trying to port this example to Bridge.jl. Can I check that the noise term is correct (ie matching my other example)?

function Bridge.σ(t, u, P::SIR)
    (S,I,R) = u
    N = S + I + R
    ifrac = P.β*P.c*I/N*S
    rfrac = P.γ*I
    return @SMatrix Float64[
     sqrt(ifrac)      0.0
    -sqrt(ifrac)  -sqrt(rfrac)
     0.0   sqrt(rfrac)
    ]
end

or do I need to remove the square roots?

In my StochasticDiffEq example, sometimes the states go below zero, so I have a callback to stop the solver then. Is there something comparable in Bridge.jl? I can run a single run, but the benchmark current fails because of this.

Any tips on parameter inference using Bridge.jl when data is the cumulative number of cases (see here) would be great, and would allow me to add another example.

Full gist so far here

mschauer commented 3 years ago

Bridge doesn't support callbacks. I would just simulate and reject simulated trajectories which become negative, for that to work, you can use sqrt(abs(ifrac)) to avoid triggering an error message. But can we reparametrize the SDE by using Ito's formula on log(X(t))? That would be convenient from a statistical point of view swell.

I can give you an example for inference from the cumulative number of cases, perhaps using https://github.com/mschauer/MitosisStochasticDiffEq.jl with @frankschae because that is what we are currently working on/working with.

sdwfrost commented 3 years ago

Thanks @mschauer. I have a reparameterization from Fintzi et al. that transforms to consider log(X(t)+1), but your workaround (which is obvious now that you mention it) is great for now.

https://github.com/epirecipes/sir-julia/blob/master/markdown/mbp/mbp.md

Thanks for the pointer to Mitosis.jl!