Copy information from Mamba.jl (#19)

* Copy information from Mamba.jl

* Fix latex and update license

* Put weird formatting back

* Fix rafterydiag

* Fix line length

* Small changes

* Remove whitespace

LICENSE

@@ -1,6 +1,11 @@
 MIT License
 
-Copyright (c) 2021 David Widmann
+Copyright (c) 2021 Brian J Smith, the Turing team and contributors:
+
+https://github.com/brian-j-smith/Mamba.jl/contributors
+https://github.com/TuringLang/MCMCChains.jl/contributors
+https://github.com/devmotion/MCMCDiagnosticTools.jl/contributors
+https://turing.ml/dev/team/
 
 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,14 @@ 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,18 @@
 """
     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)"))

src/heideldiag.jl

@@ -3,7 +3,13 @@
         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,28 @@
-#################### Raftery and Lewis Diagnostic ####################
+@doc raw"""
+    rafterydiag(
+        x::AbstractVector{<:Real}; q=0.025, r=0.005, s=0.95, eps=0.001, range=1:length(x)
+    )
 
-"""
-    rafterydiag(x::AbstractVector{<:Real}; q, r, s, eps, range)
+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)