Analysis and Diagnostics for MCMC

Specific Domains Statistics
21

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.

22

You might consider using MCMCChains.jl.

23

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)?

24

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
25

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
26

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)
27

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)
28

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…