# Analysis and Diagnostics for MCMC
**URL:** https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104
**Category:** Statistics
**Created:** [April 7, 2017, 2:07pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104 "2017-04-07T14:07:37Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@jacobcvt12](https://discourse.julialang.org/u/jacobcvt12)
#### Post date: [April 7, 2017, 2:07pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/1 "2017-04-07T14:07:37Z")
Is there a package in Julia similar to the `coda` package in R? Most of the functions in `coda` are pretty straightforward, but it’s nice to have a bunch of MCMC diagnostics/analysis functions all grouped together. I see that the documentation for `Klara` has an [MCMC Stats section](http://klarajl.readthedocs.io/en/latest/11.html), but no information.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [April 7, 2017, 3:21pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/2 "2017-04-07T15:21:55Z")
Mamba.jl has some diagnostic functions.
It would be great if in the long run one would factor out
1. a common MCMC format,
2. diagnostic for that format
into (two) separate packages. If you are interested, this would be a great project.
---
### Author:  [@jacobcvt12](https://discourse.julialang.org/u/jacobcvt12)
#### Post date: [May 3, 2017, 1:49pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/3 "2017-05-03T13:49:38Z")
Hi @Tamas_Papp thanks for the reply. Agreed - it would be great for Julia to have a standalone package for an MCMC object and diagnostics and visualizations methods for this type. Would it be acceptable to begin by pulling some of the existing code from Mamba.jl? Or should this be written from scratch?
---
### Author:  [@Nosferican](https://discourse.julialang.org/u/Nosferican)
#### Post date: [March 25, 2018, 8:05pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/4 "2018-03-25T20:05:19Z")
Would be nice to have the MCMC methods play with `StatsBase.StatisticalModel`. I imagine a good integration with solvers (for Maximum a Posterior), and `Distributions.jl`.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 26, 2018, 6:26am UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/5 "2018-03-26T06:26:37Z")
I am not sure how that would work. Typically, one uses Bayesian methods for problems which are more complex than those one would handle in a GLM framework. In my experience.
Instead of integrating things under a single umbrella framework, I prefer a _modular_ approach:
1. just coding the likelihood/posterior, ideally in an AD-compatible manner (potentially using Distributions.jl, of course),
2. optimizing that with [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) or [BlackBoxOptim.jl](https://github.com/robertfeldt/BlackBoxOptim.jl), or whatever is the most convenient,
3. when doing MCMC, just getting vectors from the posterior, for which I can calculate \hat{R} or ESS.
I think this way because when I run into problems with MCMC, it is much easier to debug small independent parts, than the insides of some framework.
BTW, on the original question: since then I have written and registered
[https://github.com/tpapp/MCMCDiagnostics.jl](https://github.com/tpapp/MCMCDiagnostics.jl)
---
### Author:  [@Nosferican](https://discourse.julialang.org/u/Nosferican)
#### Post date: [March 26, 2018, 6:57am UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/6 "2018-03-26T06:57:35Z")
Cool suite of packages! I was actually thinking mostly of Bayesian / GLM problems (e.g., Stata’s `bayes:`). Those could potentially be a low hanging fruit with tools for more complex scenarios.
---
### Author:  [@alee](https://discourse.julialang.org/u/alee)
#### Post date: [June 1, 2018, 7:08pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/7 "2018-06-01T19:08:23Z")
It would be very nice to have a small number of methods for estimating asymptotic variance / ESS associated with particular functions [since these are one-to-one]. Particularly batch means, overlapping batch means, and a spectral variance approach mentioned in
Flegal, J. M., & Jones, G. L. (2010). [Batch means and spectral variance estimators in Markov chain Monte Carlo](https://projecteuclid.org/euclid.aos/1266586622). The Annals of Statistics.
These are relatively simple, and then I would certainly use the MCMCDiagnostics package for some of my work. I would be happy to try and contribute if you are interested, although I don’t know when I will have time.
---
### Author:  [@jacobcvt12](https://discourse.julialang.org/u/jacobcvt12)
#### Post date: [June 1, 2018, 8:32pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/8 "2018-06-01T20:32:01Z")
Picking back up on this, it would be really nice to get this going. I think what @Tamas_Papp has provided in his package is a really nice start. I also really like his idea of having two packages - one for a common MCMC format, and one for diagnostic operating on that format.
I think a first step should be to get a package for storing chains. I also feel like these packages would probably get better traction if they were hosted on JuliaStats in github. Any thoughts there?
I would be happy to start the work - any volunteers to help with it? I’ll start a package and share the link this weekend.
---
### Author:  [@mschauer](https://discourse.julialang.org/u/mschauer)
#### Post date: [June 1, 2018, 9:04pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/9 "2018-06-01T21:04:56Z")
I am interested. Things I use or have use for are: Online statistics, ability to deal with vectors of static arrays, compute credibility bands and uncertainty measures and visualise them.
For me the interface I would need would be to basically create an mcmc chain
```julia
samples = MCSamples(x0, options)
...
while cond
...
push!(samples, vector, accepted)
end
```
Useful options for the MCSampler object would be to to keep only subsamples, to keep only accepted states and their count, or certain coordinates, or to do only online statistics (running mean and variance). Here MCSamples could be an abstract type with different subtypes depending on the options.
I have some preliminary things for mcmc online statistics in
[https://github.com/mschauer/Bridge.jl/blob/master/src/mclog.jl](https://github.com/mschauer/Bridge.jl/blob/master/src/mclog.jl)
---
### Author:  [@alee](https://discourse.julialang.org/u/alee)
#### Post date: [June 1, 2018, 10:17pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/10 "2018-06-01T22:17:47Z")
I think that there will be potentially a very large number of packages implementing different MCMC algorithms. Some of them may be so specialized that they can’t use a simple Diagnostics package. However, for many cases, one ultimately wants asymptotic variances or effective sample sizes associated with one or more univariate functions, e.g. in order to provide asymptotically exact confidence intervals.
So I agree with Tamas’ approach also for the reason that one can keep the scope focused. One could have separate methods for doing online estimation, as these would necessarily have to be somewhat different.
Here is a simple batch means algorithm for the asymptotic variance, with a test.
[https://gist.github.com/awllee/fe8411179d280da65a4dfc5d6af0b120](https://gist.github.com/awllee/fe8411179d280da65a4dfc5d6af0b120)
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 2, 2018, 12:36pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/11 "2018-06-02T12:36:51Z")
> [@alee](#):
>
> Particularly batch means, overlapping batch means, and a spectral variance approach mentioned in
I use
\tau = \frac{\text{effective sample size}}{\text{sample size}}
(from `MCMCDiagnostics.ess_factor_estimate`) for diagnostics, not for quantification of the variance. This is because, asymptotic guarantees aside, it can have a relative large sample variance for chains of a few thousand draws, even with good algorithms. I did some experiments [here](https://tpapp.github.io/post/ess-sampling/).
Because of this, I am distrustful of algorithms that terminate conditional on this statistic, no matter how it is obtained. MCMCDiagnostics.jl uses the relatively simple autocorrelation cutoff method, which has proven to be very robust (references are in the docstrings). When 0.2 \< \tau in difficult problems, I am reasonably happy about the outcome, otherwise I get suspicious and experiment with variable transformations.
> [@alee](#):
>
> I would be happy to try and contribute if you are interested
Contributious are always welcome, but I am not sure this is the right package for online methods. Maybe [OnlineStats.jl](https://github.com/joshday/OnlineStats.jl)?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 2, 2018, 12:44pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/12 "2018-06-02T12:44:52Z")
> [@mschauer](#):
>
> For me the interface I would need would be to basically create an mcmc chain
In practice, I find that I like to work with the following two representations:
1. A vector of vectors, in \mathbb{R}^n. These are the untransformed, “raw” coordinates in the parameter space, before the application of various transformations. I use these for convergence diagnostics, eg \hat{R}.
2. The above transformed into variables. Eg for a variable on \> 0, one would use `exp`, [ContinuousTransformations.jl](https://github.com/tpapp/ContinuousTransformations.jl) implements quite a few of these, also for vectors of such variables. I like to use these for summarizing inference (quantiles, HPD regions), and posterior predictive checks.
A package about these is in the works, but I am waiting for `v0.7` because named tuples provide such a nice interface for (2) above. I will announce it here soon.
Also, I have
> **[GitHub - tpapp/HighestDensityRegions.jl: Julia library for calculating...](https://github.com/tpapp/HighestDensityRegions.jl)**
>
> Julia library for calculating Highest Density Regions - GitHub - tpapp/HighestDensityRegions.jl: Julia library for calculating Highest Density Regions
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 2, 2018, 12:46pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/13 "2018-06-02T12:46:44Z")
> [@mschauer](#):
>
> ability to deal with vectors of static arrays
On this: any reasonable implementation should deal with all `<: AbstractVector`, but I don’t see the advantage of static arrays for MCMC analysis. All the operations are super-cheap and I don’t think it is worth the specialization overhead.
---
### Author:  [@alee](https://discourse.julialang.org/u/alee)
#### Post date: [June 2, 2018, 12:48pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/14 "2018-06-02T12:48:32Z")
To be honest, I am personally not as interested in diagnostics for conditional termination.
But I am interested in having widely used and consistent methods available for computing the integrated autocorrelation time / asymptotic variance. This is anyway what you are using to compute your diagnostic, and is a good measure to provide alongside MCMC estimates and to assess the performance of different Monte Carlo Markov chains.
I am also not that interested in online methods, at least as a first aim. The main thing I would be interested in is the availability of the different basic methods for estimating the asymptotic variance as a post-processing step, which requires very little complicated design or specialization.
---
### Author:  [@mschauer](https://discourse.julialang.org/u/mschauer)
#### Post date: [June 2, 2018, 12:53pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/15 "2018-06-02T12:53:22Z")
I don’t get you. I argue against overspecialisation. MCMC state logging methods should be at least as general as to be able to accept state vectors whose elements are geometric points. If it is restricted to `Vector{Float64}`.
---
### Author:  [@alee](https://discourse.julialang.org/u/alee)
#### Post date: [June 2, 2018, 1:01pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/16 "2018-06-02T13:01:43Z")
In many cases, at least for estimating the asymptotic variance which is the main algorithmic part, specialization should be free. In the gist I posted, the function is
```julia
function estimateAvar(xs::Vector{T}, f::F = x -> x) where F<:Function where T
```
which I have used before with T being an SVector or MVector, and really it could be anything.
Of course, one might want to further specialize beyond Vector{T}, and I assume that is feasible.
Anyway, my main question is whether different asymptotic variance estimators should be part of MCMCDiagnostics.jl or if they belong in another package.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 2, 2018, 1:23pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/17 "2018-06-02T13:23:33Z")
> [@alee](#):
>
> Anyway, my main question is whether different asymptotic variance estimators should be part of MCMCDiagnostics.jl or if they belong in another package.
As I said, I am happy to accept PRs, except for online algorithms (in the style of OnlineStats.jl), which should go there instead IMO.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 2, 2018, 1:27pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/18 "2018-06-02T13:27:56Z")
> [@mschauer](#):
>
> I argue against overspecialisation. MCMC state logging methods should be at least as general as to be able to accept state vectors whose elements are geometric points.
I am sorry I don’t fully understand. If you accept all kinds of inputs, overspecialization (of the generated code) is precisely what you get. Eg I work with `StaticVector` types in `DynamicHMC`, when I implement algorithms that benefit from it. Currently these generate specialized methods for diagnostics.
But an MCMC diagnostics package may just want to (quietly, automatically) convert everything to `Vector{Float64}`, otherwise specialized code is generated for everything, which may not be worth it. In any case, this is not a major cost ATM, especially on `v0.7`, so I may not change it unless someone opens an issue.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 2, 2018, 1:55pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/19 "2018-06-02T13:55:46Z")
Also, a nice (if a bit NUTS-specific) thread on diagnostics recently on the Stan forum:
> **[Best way to identify pathologies from divergent transitions (general step by...](https://discourse.mc-stan.org/t/best-way-to-identify-pathologies-from-divergent-transitions-general-step-by-step-workflow/4358)**
>
> Hello, following the post https://groups.google.com/forum/#!topic/stan-users/Du2vwAI5UD0 I am interested in learning how to do what @betanalpha suggested At that most general level diverges occur when the posterior on the unconstrained scaled...
I plan to implement these in Julia soon.
---
### Author:  [@jacobcvt12](https://discourse.julialang.org/u/jacobcvt12)
#### Post date: [June 2, 2018, 3:18pm UTC](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104/20 "2018-06-02T15:18:51Z")
Just to get the process started, I created package [here](https://github.com/jacobcvt12/Chains.jl) for storing chains. I think the first step is setting up a struct that essentially holds a vector of vectors, as stated. After that, we can start to expand to have things like `push!`. I’d also like to figure out a way to store chains that have birth/death, but I’ve never come across a good method for that. Any ideas?
[Next page](https://discourse.julialang.org/t/analysis-and-diagnostics-for-mcmc/3104.md?page=2)