I’m currently playing around with DifferentialEquations.jl and Catalyst.jl for simulating the paths of Markov Jump Processes, with an end goal of eventually using it to prototype some Bayesian inference algorithms.
Having followed the documentation here, I’m able to simulate paths of my MJP of interest, and plot their solution with Plots.jl, as below.
using Catalyst, DifferentialEquations lotvol_model = @reaction_network begin c1, x1 --> 2x1 c2, x1 + x2 --> 2x2 c3, x2 --> 0 end c1 c2 c3 rates = (1.0, 0.005, 0.6) prob = DiscreteProblem(lotvol_model, [70, 80], (0.0, 50.0), rates) jump_prob = JumpProblem(lotvol_model, prob, Direct()) sol = solve(jump_prob, SSAStepper()) using Plots; plot(sol)
and everything is looking great so far.
What I’d like to do next is take my solution, subsample it at integer times, and be able to plot that in the same format (as this is the sort of data upon which I will eventually be running inference algorithms).
Thus far, I’ve tried both
sol_int_times = [ sol(t) for t in 1:50]
sol_int_times = [(t, sol(t)) for t in 1:50]
With the former attempt, I can plot it, but it forgets the temporal structure of the solution, and so is not really what I’m looking for. The issue with the latter attempt is that it won’t plot at all. Upon inspecting
sol_int_times, the time and position components seem to be nested unevenly, or something like this - I don’t completely understand, but clearly I’m transforming things wrongly somehow.
Anyways, I would like to be able to subsample the solution to the MJP at the specified times, and have it stored in a format which will plot cleanly (and hopefully, be easy to work with downstream). What’s the correct syntax for doing so?