# Which MCMC library for continuous time stochastic process inference?

**URL:** <https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422>\
**Category:** Statistics\
**Created:** [January 11, 2017, 5:57pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422 "2017-01-11T17:57:43Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![nbren12](https://avatars.discourse-cdn.com/v4/letter/n/ac91a4/32.png) [@nbren12](https://discourse.julialang.org/u/nbren12)\
**Post date:** [January 11, 2017, 5:57pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/1 "2017-01-11T17:57:43Z")

</div>

I am just getting started in julia, and there seem to be many options for doing bayesian inference. I have been looking at Mamba and Klara, but I am not sure which is better for what I want to do: inference using continuous time Markov processes observed at discrete intervals.

In this scenario, writing the likelihood function requires a matrix exponential function. HMC would then require computing the derivatives of the matrix exponential in a numerically stable fashion. This is apparently not a trivial thing to do (see this STAN issue: github ([https://github.com/stan-dev/math/issues/32](https://github.com/stan-dev/math/issues/32)). Has anybody here used julia for my application?

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [January 11, 2017, 6:42pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/3 "2017-01-11T18:42:14Z")

</div>

Can we include reasons why posts are withdrawn?

---

<div class="post-metadata">

**Author:** ![Tem\_Pl](https://avatars.discourse-cdn.com/v4/letter/t/8dc957/32.png) [@Tem\_Pl](https://discourse.julialang.org/u/Tem_Pl)\
**Post date:** [January 11, 2017, 6:44pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/4 "2017-01-11T18:44:05Z")

</div>

Can sequential monte carlo work for this? Turing.jl has a SMC sampler

There is also this : [https://github.com/JuliaPOMDP](https://github.com/JuliaPOMDP)

---

<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:** [January 11, 2017, 6:50pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/5 "2017-01-11T18:50:07Z")

</div>

> [@johnmyleswhite](#):
>
> Can we include reasons why posts are withdrawn?

I was short-sighted and pointed him to a continuous-time stochastic process library

> **[GitHub - sdwfrost/Gillespie.jl: Stochastic Gillespie-type simulations using...](https://github.com/sdwfrost/Gillespie.jl)**
>
> Stochastic Gillespie-type simulations using Julia. Contribute to sdwfrost/Gillespie.jl development by creating an account on GitHub.

But that’s not what he needed. Not sure if this response is cluttering the discussion now, you can remove it if you like.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 11, 2017, 6:57pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/6 "2017-01-11T18:57:22Z")

</div>

Hi, I am currently working on such a model. There are various issues and workarounds:

1. For the matrix exponential, you can remove the balancing code from the Julia implementation, and do the rest in native Julia. This should be OK if your model transition matrix is not too unbalanced. This works with `ForwardDiff.jl`, and presumably reverse diff, but I could not get that working [because of an issue](https://github.com/JuliaStats/Klara.jl/issues/152).
2. I made a library for HMC likelihood calculations, available [here](https://github.com/tpapp/HiddenMarkovChains.jl). A test fails, I will plan to fix it next week when I am less busy. Otherwise it is fine, uses log probabilities for numerical stability. Works with AD.
3. Mamba is more mature, Klara is more hands-on. Mamba has no AD, Klara does, but since I cannot get reverse diff working at the moment, I am using Mamba, since forwarddiff AD is actually not much of an improvement.

I am happy to answer questions about details.

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [January 11, 2017, 7:20pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/7 "2017-01-11T19:20:17Z")

</div>

> [@ChrisRackauckas](#):
>
> I was short-sighted and pointed him to a continuous-time stochastic process library

I wouldn’t worry about it. If anything, my preference is that nothing is ever deleted under any circumstances.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 11, 2017, 7:32pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/8 "2017-01-11T19:32:11Z")

</div>

Forgot to mention two important differences between Klara and Mamba:

1. in Klara, you can work directly with the likelihood function, while for Mamba, you [need to use a trick](https://github.com/brian-j-smith/Mamba.jl/issues/110#issuecomment-263771452), since it prefers DAGs for model construction. But he trick is simple.
2. Mamba has much better facilities for convergence analysis.
3. Neither have tools for posterior predictive checks, so I am working on a library that does it (you can find it among my repos, but it is heavily WIP).

BTW, I (think I) fixed the heisenbug in the tests, the tests are still running:  
[https://github.com/tpapp/HiddenMarkovChains.jl](https://github.com/tpapp/HiddenMarkovChains.jl)

---

<div class="post-metadata">

**Author:** ![nbren12](https://avatars.discourse-cdn.com/v4/letter/n/ac91a4/32.png) [@nbren12](https://discourse.julialang.org/u/nbren12)\
**Post date:** [January 11, 2017, 8:01pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/9 "2017-01-11T20:01:16Z")

</div>

Thanks for your tips. Right now I am just trying to get a simple two-state discrete time Markov chain working in mamba. But I am having some difficulty figuring out how to use it.

Do you know if it possible to implement `gradexpm` by hand following the strategy used by theano, and then use automatic differentiation for everything _except_ `exmp`?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 11, 2017, 8:24pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/10 "2017-01-11T20:24:30Z")

</div>

Sorry, I don’t know the approach used by Theano and I have no time at the moment to read up on it. You can always use AD for steps implemented in pure Julia.

---

<div class="post-metadata">

**Author:** ![nbren12](https://avatars.discourse-cdn.com/v4/letter/n/ac91a4/32.png) [@nbren12](https://discourse.julialang.org/u/nbren12)\
**Post date:** [January 11, 2017, 8:58pm UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/11 "2017-01-11T20:58:36Z")

</div>

I was just wondering if I can use AD and plugin a hand-implemented gradient for a certain function.

---

<div class="post-metadata">

**Author:** ![nbren12](https://avatars.discourse-cdn.com/v4/letter/n/ac91a4/32.png) [@nbren12](https://discourse.julialang.org/u/nbren12)\
**Post date:** [January 12, 2017, 4:24am UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/12 "2017-01-12T04:24:39Z")

</div>

just to answer my own question, it looks like a custom derivative function in ForwardDiff.jl can be implemented by doing something similar to the code at the bottom of [dual.jl](https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/src/dual.jl).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 12, 2017, 5:56am UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/13 "2017-01-12T05:56:43Z")

</div>

In principle yes, and it is simple for ℝ → ℝ functions, but last time I checked there was an issue with matrix to matrix functions, in particular, it was my impression that you need to wait for [this](https://github.com/JuliaDiff/ForwardDiff.jl/pull/165). But ask the `ForwardDiff.jl` crowd, they are really friendly and helpful.

That said, forward AD is not what you want for likelihoods/posteriors, which are ℝⁿ → ℝ. You want reverse mode AD.

And finally, if you are working with a 2x2 system, you can calculate `expm` in closed form. This may be a good way to get started and explore the tools available at the moment.

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [January 12, 2017, 8:14am UTC](https://discourse.julialang.org/t/which-mcmc-library-for-continuous-time-stochastic-process-inference/1422/14 "2017-01-12T08:14:42Z")

</div>

This maybe of relevance:

[https://github.com/JuliaDiff/ForwardDiff.jl/pull/165](https://github.com/JuliaDiff/ForwardDiff.jl/pull/165)