Continous Markov Processes

New to Julia
1

Hi,

I just came across the “ContinuousTimeMarkov” written by Tom Pappas. Unfortunately, there are no examples in the repository and no documentation that I could find.

Could somebody, possibly Tom, provide me with one or more examples of how to use the package (without all the @test macros, which confuse me.

Thank you,

Gordon

2

Do you mean the package GitHub - tpapp/ContinuousTimeMarkov.jl: Julia package for Markov chains in continuous time. by Tamas Papp?
Please provide links when referring to packages.

When a package has no examples or docs, it usually means it’s designed to be work in progress that the author is not ready to release yet.

But in fact often the tests are the best way to see how to use the package.
The macro @test <statement> just tests if <statement> is true or not.

3 
julia> using Test

julia> @test 1 == 1
Test Passed

julia> @test 1 == 2
Test Failed at REPL[3]:1
  Expression: 1 == 2
   Evaluated: 1 == 2
ERROR: There was an error during testing
4

Here’s a tutorial on continuous time Markov processes:

https://diffeq.sciml.ai/latest/tutorials/discrete_stochastic_example/

In the context of epidemiology, there’s a repo with a lot of nice examples (GitHub - epirecipes/sir-julia: Various implementations of the classical SIR model in Julia) and https://github.com/epirecipes/sir-julia/blob/master/markdown/rn_mtk/rn_mtk.md is an example of running through the range of models from ODEs, SDEs, and SSAs. The DiffEq benchmarks have quite a few examples as well of large-scale SSAs like:

https://benchmarks.sciml.ai/html/Jumps/Diffusion_CTRW.html

We don’t have any examples on the newest method since it was just released, but we now have post-leap adaptive tau-leaping and regular-jump Euler-Maruyama for jump diffusions, so if that’s what you’re looking for that could be helpful too.

Another nice related package on this topic is Gillespie.jl:

We probably need to increase the discoverability of these tools.

1 Like
5

Thanks, Chris. I have discovered all of these in the past few days and am in the process of getting some of the Julia code in epirecipe to run. They were written for Julia 0.6 I believe. Biosimulation.jl is also a potentially useful package.

By the way, are jump processes the same thing as Continuous Time Markov Chains?

6

Yes, jump processes are continuous-time Markov chains, and Gillespie’s SSA algorithm is one classic method for simulating them.

Which website? The refreshed one should be good since it was just made. The older epirecipes were made for v0.6 but should still work (if you add using ParameterizedFunctions)

7

Which website?

http://epirecip.es/epicookbook/

Thanks for your links,

8

Simon’s newer version is here: https://github.com/epirecipes/sir-julia

9

Thanks. On the original site, there are a few network programs in R. You would not happen to know of any Julia implementations?

Have you ever worked with the Household-Network-workplace model? That is what I am looking at to work with along with a generalization of SEIR to perhaps account for assumptions cases. I would like to use NeuralODE to help fit the data. Seems like a good problem,

10

I read the tutorials. In the context of SIR, I see where the transition rates fill in, but I do not see where the statistics Al properties of the time jumps are set. If the intervals between time events vary according to an exponential distribution with mean \lambda, where is that specified when using the various stochastic components of DifferentialEquations.jl?

Thanks,

11

I do not know that model.

I plan to show in good detail how to do this soon, but it’s not entirely straightforward.

12

I look forward to that day, Chris. But if the user is not setting the value, you must be doing this in the code. In a SIR model doesn’t this parameter of the exponential distribution relate to the average infection length? So you explicitly allow setting the model contact rate, but not the infection rate, which could also,depend on time. Why was this decision made? What source code could I examine to develop some intuition? There is a lot I still do not understand.

This is a whole new world for me.

I found the following callback.jl function in DiffEqJump.jl .

@inline function time_to_next_jump(u,p,t,rate)
randexp()/rate(u,p,t)
#rand(Exponential(1/rate(u,p,t)))
end

ConstantRateJumpCallback(next_jump,c::ConstantRateJump,end_time) = ConstantRateJumpCallback(next_jump,end_time,c.rate,c.affect!,c.save_positions)

DiscreteCallback(c::ConstantRateJumpCallback) = DiscreteCallback(c,c,c.save_positions)

function DiscreteCallback(u,p,t,c::ConstantRateJump,end_time)
next_jump = time_to_next_jump(u,p,t,c.rate)
DiscreteCallback(ConstantRateJumpCallback(next_jump,c,end_time))
end

Looking more carefully at DiffEqJump.jl, I noticed

@inline function (p::ConstantRateJumpCallback)(integrator) # affect!
p.affect!(integrator)
p.next_jump = integrator.t + time_to_next_jump(integrator.t,integrator.u,p.rate)
if p.next_jump < p.end_time
add_tstop!(integrator,p.next_jump)
end
end

I do not understand how the call to time_to_next_jump() can only have three arguments. It definition has four with no default arguments.

13

The parameters are all of the rates for each compound Poisson process. These can depend on time, but then you can no longer use Gillespie’s algorithm so you’d have to use a VaraibleRateJump or tau-leaping.

The call seems to have 4 arguments. What version of the code are you looking at?https://github.com/SciML/DiffEqJump.jl/blob/master/src/aggregators/direct.jl#L89

14

I do not know what verion I am looking up, but using the link you provide, I doubled back to callbacks.jl:
https://github.com/SciML/DiffEqJump.jl/blob/master/src/callbacks.jl

Lines 15: 3 arguments.

p.next_jump = integrator.t + time_to_next_jump(integrator.t,integrator.u,p.rate)

Line 31: 3 arguments

next_jump = time_to_next_jump(u,p,t,c.rate)

Lines 21: 4 arguments

@inline function time_to_next_jump(u,p,t,rate)

Regarding VariableRateJump: The documentation implied that this is what you use when the rate depends on the solution variables, implying that the rate would vary between two Poisson time points since the solutions change between two Poisson time points. Is this true? Please excuse the improper use of words since I am only learning this.
:

15

That part isn’t being used because jumps have to be aggregated. @isaacsas I think we can just delete that portion?

Yes. The other jumps make the assumption that the rates are constant between events, and so you can use an exponential to get the next event time. VariableRateJumps drop that assumption, at a cost.

16

Thanks. It is making more sense now. So the SDE library is a discrete process with fixed time steps?

Where can I read about the aggregation process?

I will stop for now :). I appreciate your patience.

Gordon

17

It’s continuous and discrete (and mixed) with adaptive and fixed time steps.

The basic one is the Gillespie direct method. We really need to do a refresher on our tutorials here, because https://docs.sciml.ai/latest/tutorials/discrete_stochastic_example/ only scratches the surface, and all of https://docs.sciml.ai/latest/types/jump_types/#Constant-Rate-Jump-Aggregators-1 are different SSAs.

This is a good overview of Direct: http://users.abo.fi/ipetre/advcompmod/Lecture_6.pdf

18

I would like to mention PiecewiseDeterministicMarkovProcesses.jl as well

2 Likes
19

Indeed, that’s a good alternative to VariableRateJumps

20

If that is really me (that’s the most creative way I have seen my name mangled yet): the package was used for some exploratory work for research, and is currently dormant. It wasn’t a full suite for CTM processes anyway, just a subset of the relevant functionality.

3 Likes