Analysis and Diagnostics for MCMC

Did anyone ever make any progress on this? I’d like to have an estimator that’s consistent with what is outlined in the Flegal & Jones paper.

You might consider using MCMCChains.jl.

It’s a little unclear to me from the documentation which method it uses for estimating the statistic. Also, I don’t want to switch over entirely to MCMCChains - is there a way to call the autocorrelation time function on a given time series (of scalars)?

I understand why you may not want to commit the a Chain object in MCMCChains. About a year ago, the package maintainers split the diagnostics from MCMCChains and put them in standalone package: GitHub - TuringLang/MCMCDiagnosticTools.jl

I think this should have the features you want.

1 Like

Indeed, this might be what I want, but I’m unsure of how to use it without falling back to the chain datastructure. For instance, how would I get MCMCDiagnosticTools to work with the following?

samples = cumsum(randn(10000));
ess_hat(samples)

just generates an error:

ERROR: MethodError: no method matching ess_rhat(::Vector{Float64})
Closest candidates are:
  ess_rhat(::AbstractArray{<:Union{Missing, Real}, 3}; method, maxlag) at ~/.julia/packages/MCMCDiagnosticTools/kKCOf/src/ess.jl:208
Stacktrace:
 [1] top-level scope
   @ REPL[12]:1
1 Like

According to the documentation, the function expects a three dimensional array.

help?> ess_rhat
search:

  ess_rhat(
      samples::AbstractArray{<:Union{Missing,Real},3}; method=ESSMethod(), maxlag=250
  )

  Estimate the effective sample size and the potential scale reduction of the samples of shape (draws, parameters, chains) with the
  method and a maximum lag of maxlag.

  See also: ESSMethod, FFTESSMethod, BDAESSMethod

Example

using MCMCDiagnosticTools: ess_rhat

x = rand(1000, 10, 3)

ess,rhat = ess_rhat(x)

In your example, you do the following:

using MCMCDiagnosticTools: ess_rhat

samples = cumsum(randn(10000, 1, 1), dims=1);

ess,rhat = ess_rhat(samples)

You can also add more dimensions outside of cumsum:

samples = [samples ;;;]
ess,rhat = ess_rhat(samples)

I had the same question as

Although this example was extremely helpful:

I still find this instruction in the documentation confusing,

in that I am not sure what cumsum(randn(10000, 1, 1), dims=1) is supposed to represent given randn(10000, 1, 1) at hand, which to me already looks like a chain of random variables of length 10000…