# MCMC landscape
**URL:** https://discourse.julialang.org/t/mcmc-landscape/25654
**Category:** Statistics
**Tags:** question
**Created:** [June 25, 2019, 3:23pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654 "2019-06-25T15:23:28Z")
**Posts on this page:** 20
**Page:** 5
### Author:  [@aplavin](https://discourse.julialang.org/u/aplavin)
#### Post date: [August 11, 2019, 5:47pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/81 "2019-08-11T17:47:44Z")
Thanks for suggesting to help! I’m definitely not an expert in inference methods myself, so not going to implement them in Julia. What I’m possibly able to do is integrating already existing libraries in other languages (e.g. PolyChord I mentioned above) to take advantage of the Turing infrastructure of defining models.
Looking at the importance sampling implementation as the simplest reference I still don’t completely grasp how samplers should work in Turing. For example, the polychord library takes two separate functions as the model specification, and I get confused how to get the first one already. The `prior` function should take N values - the number of scalar parameters in the model - and return the same number of values converted from the unit cube to the model prior distribution. As I understand Turing interface, I would need to run the `Model` object once and set up `assume` method to gather prior distributions for all variables in one array. Then the `prior` function I give to the library just calls `invcdf` for all distributions in that array. Does that sound at least generally the right way to go?
---
### Author:  [@trappmartin](https://discourse.julialang.org/u/trappmartin)
#### Post date: [August 12, 2019, 2:51pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/82 "2019-08-12T14:51:13Z")
So the idea behind the assume and the observe statements is as follows.
Assume is called for every parameter of the model and allows you to specify what happens for model parameters. In the case of importance sampling, we sample the parameter value from the prior distribution.
See:
[https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/is.jl#L60-L64](https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/is.jl#L60-L64)
In the case of HMC, we retrieve the new parameter value from the internal data structure and evaluate the log pdf of the prior given in the transformed space.
See:
[https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/hmc.jl#L606-L612](https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/hmc.jl#L606-L612)
The observe function is used for observations and usually returns the log pdf or log pmf of an observation for the observation model. See:
[https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/is.jl#L66-L69](https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/is.jl#L66-L69)
As you see in the importance sampling example, this is implemented using multi-dispatch. So to implement an interface to a custom sampler you need to define a new type which sub-types `InferenceAlgorithm` and a `Sampler` function which returns a new sampler, see: [https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/is.jl#L38](https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/inference/is.jl#L38)
Then you can implement the `sample` function which is called by the user. We will soon have a better API for this.
I’m not 100% sure I understand how to use the PolyChord library (how do you tell the library which prior distributions to use?), but here is an attempt on how you can interface the PolyChord library in your `sample` function. (Note that we don’t have an API to get the number of parameters as of yet)
```nohighlight
vi = VarInfo(model)
# get the parameter values
initial_params = mapreduce(m -> m.vals, vcat, vi.metadata)
# Have a look here: https://github.com/TuringLang/Turing.jl/blob/abbc8c011f7a51192378d8681208083b649dc027/src/core/RandomVariables.jl#L118
# call PolyChord
params = PolyChord.prior(rand(length(initial_params)))
# set the new parameter values
c = 1
for m in vi.metadata
vn = m.vns # This will likely break cause issues if the user uses the same symbol in different contexts of the same model
L = length(m.vals)
vi[vn] = params[c:(c+L-1)]
c += L
end
```
Does that make sense to you?
Apologies if I misunderstood you, I should really read the paper first. 😅
---
### Author:  [@aplavin](https://discourse.julialang.org/u/aplavin)
#### Post date: [August 12, 2019, 10:11pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/83 "2019-08-12T22:11:20Z")
Well, I guess I just described it in a little bit confusing way. The workflow for polychord library is basically like this:
```julia
# `prior` converts from uniform cube to model parameters
# e.g. if we have two parameters, one with Normal(0, 1) prior and the other with Normal(5, 2):
prior(cube::Vector{Float})::Vector{Float} = invcdf.([Normal(0, 1), Normal(5, 2)], cube)
# computes likelihood for given model parameters
loglike(params::Vector{Float})::Float = ...
run_polychord(prior, loglike, options)
```
It looks like this approach and the one in Turing have “inverted” directions of control flow, in a sense.
---
### Author:  [@trappmartin](https://discourse.julialang.org/u/trappmartin)
#### Post date: [August 13, 2019, 9:20am UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/84 "2019-08-13T09:20:06Z")
I see I think this could be done much easier than in my example. I also had a look at the nested sampling paper and I agree it is a very nice should add this to Turing. Could you open an issue on Turing so that we can discuss the details there?
You can also have a look at [https://github.com/TuringLang/Turing.jl/issues/886](https://github.com/TuringLang/Turing.jl/issues/886) and contribute some ideas on how to improve the API from a user perspective.
---
### Author:  [@mlanghinrichs](https://discourse.julialang.org/u/mlanghinrichs)
#### Post date: [February 27, 2020, 1:04pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/85 "2020-02-27T13:04:24Z")
I would like to revive this chat a bit, if possible. I’m coming from Python’s MCMC environment and want to try out Julia for some new projects. I have used the great `emcee` and `dynesty` packages for sampling a self-defined loglikelihood function and nested sampling, respectively, and I’m looking for similar options in Julia. So it kind of fits @davidanthoff initial questions and also the later discussion on nested sampling. Are there any news/developements in the meantime?
More specifically, I was wondering how I could use Turing.jl to sample from a user-defined log likelihood function. If there were recent developments in this direction and, in case, if there is a minimal example somewhere?
Any help is really appreciated!
---
### Author:  [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)
#### Post date: [February 27, 2020, 1:31pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/86 "2020-02-27T13:31:50Z")
> [@mlanghinrichs](#):
>
> I would like to revive this chat a bit, if possible. I’m coming from Python’s MCMC environment and want to try out Julia for some new projects. I have used the great `emcee` and `dynesty` packages for sampling a self-defined loglikelihood function and nested sampling, respectively, and I’m looking for similar options in Julia. So it kind of fits @davidanthoff initial questions and also the later discussion on nested sampling. Are there any news/developements in the meantime?
[NestedSamplers.jl](https://github.com/mileslucas/NestedSamplers.jl) exists. It’s not yet integrated into Turing, but its relatively well designed in that it could be built into Turing with some effort. I’m also looking into building a Julia/emcee bridge that would allow for the use of `emcee` in Julia.
> [@mlanghinrichs](#):
>
> More specifically, I was wondering how I could use Turing.jl to sample from a user-defined log likelihood function. If there were recent developments in this direction and, in case, if there is a minimal example somewhere?
I mean, technically all Turing models are “user-defined”. More seriously, if you are asking about how you get a posterior if you have some likelihood function `f`, you don’t need Turing. You can go upstream to the inference packages behind Turing, like [AdvancedHMC](https://github.com/TuringLang/AdvancedHMC.jl) or [AdvancedMH](https://github.com/TuringLang/AdvancedMH.jl). You can also use [DynamicHMC](https://github.com/tpapp/DynamicHMC.jl).
Turing’s main strength is defining the model. The breadth of inference methods and modularity I think are great, but they ultimately pale in comparison to the ability to write super simple models that just work. If you already have a likelihood function and don’t need to write the model specification, you probably don’t need Turing.
---
### Author:  [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [February 27, 2020, 1:53pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/87 "2020-02-27T13:53:19Z")
To build upon cpfiffer’s response, you can create a custom log likelihood function by extending logpdf and creating a new distribution subtype from Distributions.jl. Here is an example from the [Turing documentation](https://turing.ml/dev/docs/using-turing/advanced).
---
### Author:  [@mlanghinrichs](https://discourse.julialang.org/u/mlanghinrichs)
#### Post date: [February 27, 2020, 3:43pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/88 "2020-02-27T15:43:03Z")
Thanks to you both and great work! I will take a look at the different options, AdvancedMH might be what I need since my problem is not AD-able (I guess). But I will also try to do it within Turing and get a bit used to this kind of probabilistic programming.
And nice to see that there is also some active development on nested sampling. I saw NestedSampling.jl before, but missed out NestedSamplers.jl so far. Thanks.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [February 27, 2020, 4:16pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/89 "2020-02-27T16:16:13Z")
> [@mlanghinrichs](#):
>
> my problem is not AD-able (I guess)
Obtaining gradients will have the highest payoff for convergence in a nontrivial application since it allows you to use HMC variant methods, so you may want to explore this first.
---
### Author:  [@mlanghinrichs](https://discourse.julialang.org/u/mlanghinrichs)
#### Post date: [February 27, 2020, 5:24pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/90 "2020-02-27T17:24:57Z")
> [@Tamas\_Papp](#):
>
> Obtaining gradients will have the highest payoff for convergence in a nontrivial application since it allows you to use HMC variant methods, so you may want to explore this first.
yes. Maybe you guys can help me. I’m trying to fit mathematical models to data (D), typically stochastic Markov models with discrete states and continuous time (“Gillespie” algorithm). So given a model M and parameters θ, I can simulate correct stochastic realisations from the model/stochastic process (you can think of finance stock market or so). Then I compute a (log)likelihood L = p(D | θ, M) based on the collection of many model realisations and some statistics that apply both to data and model. Then MCMC to get posterior p(θ | D, M). Might automatic differentiation be possible here?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [February 27, 2020, 5:36pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/91 "2020-02-27T17:36:49Z")
Sounds like you are using indirect inference?
If the simulation results are differentiable in the parameters, there should be no reason you cannot make it work with AD. Just use [constant random numbers](https://en.wikipedia.org/wiki/Variance_reduction#Common_Random_Numbers_(CRN)).
If this is not the case, then frankly, I would transform the model until it is 😉
---
### Author:  [@mlanghinrichs](https://discourse.julialang.org/u/mlanghinrichs)
#### Post date: [February 27, 2020, 5:43pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/92 "2020-02-27T17:43:09Z")
> [@Tamas\_Papp](#):
>
> Sounds like you are using indirect inference?
Might be, right away I cannot say for sure 😃 My field is quantitative biology, not economics, but concepts might be similar.
> [@Tamas\_Papp](#):
>
> If the simulation results are differentiable in the parameters, there should be no reason you cannot make it work with AD. Just use [constant random numbers](https://en.wikipedia.org/wiki/Variance_reduction#Common_Random_Numbers_(CRN)).
>
> If this is not the case, then frankly, I would transform the model until it is 😉
Thanks a lot! I will look into this.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [February 27, 2020, 5:52pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/93 "2020-02-27T17:52:15Z")
> [@mlanghinrichs](#):
>
> My field is quantitative biology, not economics, but concepts might be similar.
Yeah, I recall the many quantitative biology papers I read (as an economist) when learning about this methodology 😉
You may be interested in
> **[ABC and Indirect Inference](https://arxiv.org/abs/1803.01999)**
>
> This chapter will appear in the forthcoming Handbook of Approximate Bayesian Computation (2018).
> Indirect inference (II) is a classical likelihood-free approach that pre-dates the main developments of ABC and relies on simulation from a parametric...
---
### Author:  [@mlanghinrichs](https://discourse.julialang.org/u/mlanghinrichs)
#### Post date: [February 27, 2020, 6:01pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/94 "2020-02-27T18:01:56Z")
> [@Tamas\_Papp](#):
>
> Yeah, I recall the many quantitative biology papers I read (as an economist) when learning about this methodology 😉
All the same it seems, just different data 😃
Thanks for the link, I’ll try to process the stuff and see how far I come!
EDIT: There is actually a topic related to my problem ([Differentiating through a Jump Problem](https://discourse.julialang.org/t/differentiating-through-a-jump-problem/31496)).
EDIT2: Just for completeness, my current conclusions:
1. Following the linked topic above, Gillespie simulations with discrete states and discrete random selection of reactions do _not_ seem to be (automatic) differentiable. You could _make_ them continuous as a SDE for AD to work (see linked topic for more info), but that’s not the path I want to take. So I’m using samplers now without the need of a gradient; i.e. [AdvancedMH](https://github.com/TuringLang/AdvancedMH.jl) as mentioned by @cpfiffer.
If there are more non-gradient samplers in Julia, I would be glad to know! On a first glance, the currently implemented MH sampler generally works, but can have low number of effective samples.
2. Approximate Bayesian (ABC) methods should be applicable for these problems; they do not need the evaluation of a potentially costly likelihood, but are – as the name suggests – only an approximation of the posterior (how good, has to be checked case-specific). I will try them out once I run into runtime limitations. Potential Julia implementations: [ApproxBayes](https://github.com/marcjwilliams1/ApproxBayes.jl), [GpABC](https://github.com/tanhevg/GpABC.jl), [ApproximateBayesianComputing](https://github.com/eford/ApproximateBayesianComputing.jl), [ABC](https://github.com/eford/ABC.jl/blob/master/docs/src/index.md) (and maybe more).
---
### Author:  [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [February 28, 2020, 6:22pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/95 "2020-02-28T18:22:24Z")
> [@mlanghinrichs](#):
>
> If there are more non-gradient samplers in Julia, I would be glad to know! On a first glance, the currently implemented MH sampler generally works, but can have low number of effective samples.
I have found that Differential Evolution MCMC works well for non-gradient based models. I implemented a [version](https://github.com/itsdfish/AdvancedPS.jl) of DE-MCMC that might go into one of Turing’s library for particle based samplers. There is a folder with examples. It’s not as efficient as NUTS but it’s much better than AdvancedMH (Edit: MH methods in general, not AdvancedMH in particular).
---
### Author:  [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)
#### Post date: [February 29, 2020, 12:39am UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/96 "2020-02-29T00:39:25Z")
> [@Christopher\_Fisher](#):
>
> It’s not as efficient as NUTS but it’s much better than AdvancedMH.
Do you mean MH methods specifically or the implementation of it? I’m always looking to make AdvancedMH better.
---
### Author:  [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)
#### Post date: [February 29, 2020, 10:06am UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/97 "2020-02-29T10:06:28Z")
Sorry for the confusion. I was referring to MH as a method in general. Your implementation looks good to me. However, if there is something inefficient in your implementation, I probably wouldn’t be the one to identify it 🙂
---
### Author:  [@cpfiffer](https://discourse.julialang.org/u/cpfiffer)
#### Post date: [February 29, 2020, 4:30pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/98 "2020-02-29T16:30:19Z")
> [@Christopher\_Fisher](#):
>
> However, if there is something inefficient in your implementation, I probably wouldn’t be the one to identify it
Crap, I really needed someone to tell me what’s wrong 😉.
---
### Author:  [@monty](https://discourse.julialang.org/u/monty)
#### Post date: [July 2, 2021, 3:05pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/99 "2021-07-02T15:05:21Z")
I just want to add our:
> **[GitHub - madsjulia/AffineInvariantMCMC.jl: Affine Invariant Markov Chain...](https://github.com/madsjulia/AffineInvariantMCMC.jl)**
>
> Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble sampler - GitHub - madsjulia/AffineInvariantMCMC.jl: Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble sampler
---
### Author:  [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [July 2, 2021, 4:04pm UTC](https://discourse.julialang.org/t/mcmc-landscape/25654/100 "2021-07-02T16:04:16Z")
(edit: sorry, this is in reponse to @mlanghinrichs, up above) Some work I have done is a related to this, I think. I estimated a continuous time diffusion model by doing MCMC using the asymptotic Gaussian likelihood of a vector of statistics as the “quasi” likelihood. The mean and the covariance matrix of the Gaussian likelihood at given parameters are estimated by simulations from the model. The dimension of the statistics is reduced using a previously trained neural net. The method uses plain MH MCMC. This works well, contrary to what might be expected, because the theory behind the method leads to a very effective proposal density. At any rate, it is an example of MCMC for a continuous time model, based on statistics rather than the full sample. The example is at [https://github.com/mcreel/SNM/tree/master/examples/JD](https://github.com/mcreel/SNM/tree/master/examples/JD). A working paper referenced at the top directory of the repo gives explanations.
There are also simpler models in the examples. I have found this method to work quite well for all of the examples I have explored. The limitation is the time needed to sample from the model. Because a neural net must be trained, if it’s expensive to sample from the model, then it can be a little time consuming to estimate a model. The jump diffusion model took 2 or 3 days, maybe 4, on an old 32 core machine
[Previous page](https://discourse.julialang.org/t/mcmc-landscape/25654.md?page=4)