Copy information from Mamba.jl

LICENSE

@@ -1,6 +1,10 @@
 MIT License
 
-Copyright (c) 2021 David Widmann
+Copyright (c) 2021 Brian J Smith, David Widmann and contributors:
+
+https://github.com/brian-j-smith/Mamba.jl/contributors
+
+https://github.com/devmotion/MCMCDiagnosticTools.jl/contributors
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

src/gelmandiag.jl

@@ -1,5 +1,3 @@
-#################### Gelman, Rubin, and Brooks Diagnostics ####################
-
 function _gelmandiag(psi::AbstractArray{<:Real,3}; alpha::Real=0.05)
     niters, nparams, nchains = size(psi)
     nchains > 1 || error("Gelman diagnostic requires at least 2 chains")
@@ -56,7 +54,13 @@ end
 """
     gelmandiag(chains::AbstractArray{<:Real,3}; alpha::Real=0.95)
 
-Compute the Gelman, Rubin and Brooks diagnostics.
+Compute the Gelman, Rubin and Brooks diagnostics [^Gelman1992] [^Brooks1998].
+Values of the diagnostic’s potential scale reduction factor (PSRF) that are close to one suggest convergence.
+As a rule-of-thumb, convergence is rejected if the 97.5 percentile of a PSRF is greater than 1.2.
+
+[^Gelman1992]: Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472.
+
+[^Brooks1998]: Brooks, S. P., & Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of computational and graphical statistics, 7(4), 434-455.
 """
 function gelmandiag(chains::AbstractArray{<:Real,3}; kwargs...)
     estimates, upperlimits = _gelmandiag(chains; kwargs...)

src/gewekediag.jl

@@ -1,7 +1,15 @@
 """
     gewekediag(x::AbstractVector{<:Real}; first::Real=0.1, last::Real=0.5, kwargs...)
 
-Compute the Geweke diagnostic from the `first` and `last` proportion of samples `x`.
+Compute the Geweke diagnostic [^Geweke1991] from the `first` and `last` proportion of samples `x`.
+
+The diagnostic is designed to asses convergence of posterior means estimated with autocorrelated samples.
+It computes a normal-based test statistic comparing the sample means in two windows containing proportions of the first and last iterations.
+Users should ensure that there is sufficient separation between the two windows to assume that their samples are independent.
+A non-significant test p-value indicates convergence.
+Significant p-values indicate non-convergence and the possible need to discard initial samples as a burn-in sequence or to simulate additional samples.
+
+[^Geweke1991]: Geweke, J. F. (1991). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments (No. 148). Federal Reserve Bank of Minneapolis.
 """
 function gewekediag(x::AbstractVector{<:Real}; first::Real=0.1, last::Real=0.5, kwargs...)
     0 < first < 1 || throw(ArgumentError("`first` is not in (0, 1)"))
@@ -11,8 +19,7 @@ function gewekediag(x::AbstractVector{<:Real}; first::Real=0.1, last::Real=0.5,
     n = length(x)
     x1 = x[1:round(Int, first * n)]
     x2 = x[round(Int, n - last * n + 1):n]
-    z =
-        (Statistics.mean(x1) - Statistics.mean(x2)) /
+    z = (Statistics.mean(x1) - Statistics.mean(x2)) /
         hypot(mcse(x1; kwargs...), mcse(x2; kwargs...))
     p = SpecialFunctions.erfc(abs(z) / sqrt(2))
 

src/heideldiag.jl

@@ -3,7 +3,12 @@
         x::AbstractVector{<:Real}; alpha::Real=0.05, eps::Real=0.1, start::Int=1, kwargs...
     )
 
-Compute the Heidelberger and Welch diagnostic.
+Compute the Heidelberger and Welch diagnostic [^Heidelberger1983].
+This diagnostic tests for non-convergence (non-stationarity) and whether ratios of estimation interval halfwidths to means are within a target ratio.
+Stationarity is rejected (0) for significant test p-values.
+Halfwidth tests are rejected (0) if observed ratios are greater than the target, as is the case for `s2` and `beta[1]`.
+
+[^Heidelberger1983]: Heidelberger, P., & Welch, P. D. (1983). Simulation run length control in the presence of an initial transient. Operations Research, 31(6), 1109-1144.
 """
 function heideldiag(
     x::AbstractVector{<:Real}; alpha::Real=0.05, eps::Real=0.1, start::Int=1, kwargs...

src/rafterydiag.jl

@@ -1,9 +1,19 @@
-#################### Raftery and Lewis Diagnostic ####################
-
 """
-    rafterydiag(x::AbstractVector{<:Real}; q, r, s, eps, range)
+    rafterydiag(x::AbstractVector{<:Real}; q=0.025, r=0.005, s=0.95, eps=0.001, range=1:length(x))
+
+Compute the Raftery and Lewis diagnostic [^Raftery1992].
+This diagnostic is used to determine the number of iterations required to estimate a specified quantile `q` within a desired degree of accuracy.
+The diagnostic is designed to determine the number of autocorrelated samples required to estimate a specified quantile $\theta_q$, such that $\Pr(\theta \le \theta_q) = q$, within a desired degree of accuracy.
+In particular, if $\hat{\theta}_q$ is the estimand and $\Pr(\theta \le \hat{\theta}_q) = \hat{P}_q$ the estimated cumulative probability, then accuracy is specified in terms of $r$ and $s$, where $\Pr(q - r < \hat{P}_q < q + r) = s$.
+Thinning may be employed in the calculation of the diagnostic to satisfy its underlying assumptions.
+However, users may not want to apply the same (or any) thinning when estimating posterior summary statistics because doing so results in a loss of information.
+Accordingly, sample sizes estimated by the diagnostic tend to be conservative (too large).
+
+Furthermore, the argument `r` specifies the margin of error for estimated cumulative probabilities and `s` the probability for the margin of error.
+`eps` specifies the  tolerance within which the probabilities of transitioning from initial to retained iterations are within the equilibrium probabilities for the chain.
+This argument determines the number of samples to discard as a burn-in sequence and is typically left at its default value.
 
-Compute the Raftery and Lewis diagnostic.
+[^Raftery1992]: A L Raftery and S Lewis. Bayesian Statistics, chapter How Many Iterations in the Gibbs Sampler? Volume 4. Oxford University Press, New York, 1992.
 """
 function rafterydiag(
     x::AbstractVector{<:Real}; q=0.025, r=0.005, s=0.95, eps=0.001, range=1:length(x)