From 90f502dec5f116f7c12a42df1aa285137dffe6af Mon Sep 17 00:00:00 2001 From: Rik Huijzer Date: Wed, 1 Sep 2021 18:28:57 +0200 Subject: [PATCH] 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 | 7 ++++++- src/gelmandiag.jl | 11 ++++++++--- src/gewekediag.jl | 13 ++++++++++++- src/heideldiag.jl | 8 +++++++- src/rafterydiag.jl | 27 +++++++++++++++++++++++---- 5 files changed, 56 insertions(+), 10 deletions(-) diff --git a/LICENSE b/LICENSE index 3653064c..90124513 100644 --- a/LICENSE +++ b/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 diff --git a/src/gelmandiag.jl b/src/gelmandiag.jl index e080ed0e..87626672 100644 --- a/src/gelmandiag.jl +++ b/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...) diff --git a/src/gewekediag.jl b/src/gewekediag.jl index 80fcc36d..f804b62c 100644 --- a/src/gewekediag.jl +++ b/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)")) diff --git a/src/heideldiag.jl b/src/heideldiag.jl index 936a9a7f..e93291d8 100644 --- a/src/heideldiag.jl +++ b/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... diff --git a/src/rafterydiag.jl b/src/rafterydiag.jl index f580fd88..8f34f0b8 100644 --- a/src/rafterydiag.jl +++ b/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)