Skip to content

Commit

Permalink
Merge pull request #215 from stan-dev/replace-autocovariance
Browse files Browse the repository at this point in the history
autocovariance -> posterior::autocovariance
  • Loading branch information
jgabry authored Mar 25, 2023
2 parents d954501 + df72980 commit 1832085
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 23 deletions.
1 change: 1 addition & 0 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ Suggests:
ggplot2,
graphics,
knitr,
posterior,
rmarkdown,
rstan,
rstanarm (>= 2.19.0),
Expand Down
27 changes: 4 additions & 23 deletions R/effective_sample_sizes.R
Original file line number Diff line number Diff line change
Expand Up @@ -199,17 +199,13 @@ psis_n_eff.matrix <- function(w, r_eff = NULL, ...) {
#' @return MCMC effective sample size based on RStan's calculation.
#'
ess_rfun <- function(sims) {
# Compute the effective sample size for samples of several chains
# for one parameter; see the C++ code of function
# effective_sample_size in chains.cpp
#
# Args:
# sims: a 2-d array _without_ warmup samples (# iter * # chains)
#
if (!requireNamespace("posterior", quietly = TRUE)) {
stop("Please install the 'posterior' package", call. = FALSE)
}
if (is.vector(sims)) dim(sims) <- c(length(sims), 1)
chains <- ncol(sims)
n_samples <- nrow(sims)
acov <- lapply(1:chains, FUN = function(i) autocovariance(sims[,i]))
acov <- lapply(1:chains, FUN = function(i) posterior::autocovariance(sims[,i]))
acov <- do.call(cbind, acov)
chain_mean <- colMeans(sims)
mean_var <- mean(acov[1,]) * n_samples / (n_samples - 1)
Expand Down Expand Up @@ -275,18 +271,3 @@ fft_next_good_size <- function(N) {
N = N + 1
}
}

autocovariance <- function(y) {
# Compute autocovariance estimates for every lag for the specified
# input sequence using a fast Fourier transform approach.
N <- length(y)
M <- fft_next_good_size(N)
Mt2 <- 2 * M
yc <- y - mean(y)
yc <- c(yc, rep.int(0, Mt2-N))
transform <- stats::fft(yc)
ac <- stats::fft(Conj(transform) * transform, inverse = TRUE)
# use "biased" estimate as recommended by Geyer (1992)
ac <- Re(ac)[1:N] / (N^2 * 2)
ac
}

0 comments on commit 1832085

Please sign in to comment.