Skip to content

Commit 58fe5be

Browse files
committed
Merge branch 'master' of https://github.com/stan-dev/loo
2 parents fe83b61 + a2f26de commit 58fe5be

2 files changed

Lines changed: 23 additions & 5 deletions

File tree

DESCRIPTION

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,6 @@ Suggests:
4343
ggplot2,
4444
graphics,
4545
knitr,
46-
posterior,
4746
rmarkdown,
4847
rstan,
4948
rstanarm (>= 2.19.0),

R/effective_sample_sizes.R

Lines changed: 23 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -199,13 +199,17 @@ psis_n_eff.matrix <- function(w, r_eff = NULL, ...) {
199199
#' @return MCMC effective sample size based on RStan's calculation.
200200
#'
201201
ess_rfun <- function(sims) {
202-
if (!requireNamespace("posterior", quietly = TRUE)) {
203-
stop("Please install the 'posterior' package", call. = FALSE)
204-
}
202+
# Compute the effective sample size for samples of several chains
203+
# for one parameter; see the C++ code of function
204+
# effective_sample_size in chains.cpp
205+
#
206+
# Args:
207+
# sims: a 2-d array _without_ warmup samples (# iter * # chains)
208+
#
205209
if (is.vector(sims)) dim(sims) <- c(length(sims), 1)
206210
chains <- ncol(sims)
207211
n_samples <- nrow(sims)
208-
acov <- lapply(1:chains, FUN = function(i) posterior::autocovariance(sims[,i]))
212+
acov <- lapply(1:chains, FUN = function(i) autocovariance(sims[,i]))
209213
acov <- do.call(cbind, acov)
210214
chain_mean <- colMeans(sims)
211215
mean_var <- mean(acov[1,]) * n_samples / (n_samples - 1)
@@ -271,3 +275,18 @@ fft_next_good_size <- function(N) {
271275
N = N + 1
272276
}
273277
}
278+
279+
autocovariance <- function(y) {
280+
# Compute autocovariance estimates for every lag for the specified
281+
# input sequence using a fast Fourier transform approach.
282+
N <- length(y)
283+
M <- fft_next_good_size(N)
284+
Mt2 <- 2 * M
285+
yc <- y - mean(y)
286+
yc <- c(yc, rep.int(0, Mt2-N))
287+
transform <- stats::fft(yc)
288+
ac <- stats::fft(Conj(transform) * transform, inverse = TRUE)
289+
# use "biased" estimate as recommended by Geyer (1992)
290+
ac <- Re(ac)[1:N] / (N^2 * 2)
291+
ac
292+
}

0 commit comments

Comments
 (0)