@@ -202,72 +202,5 @@ psis_n_eff.matrix <- function(w, r_eff = NULL, ...) {
202202# ' @return MCMC effective sample size based on RStan's calculation.
203203# '
204204ess_rfun <- function (sims ) {
205- if (is.vector(sims )) dim(sims ) <- c(length(sims ), 1 )
206- chains <- ncol(sims )
207- n_samples <- nrow(sims )
208- acov <- lapply(1 : chains , FUN = function (i ) posterior :: autocovariance(sims [,i ]))
209- acov <- do.call(cbind , acov )
210- chain_mean <- colMeans(sims )
211- mean_var <- mean(acov [1 ,]) * n_samples / (n_samples - 1 )
212- var_plus <- mean_var * (n_samples - 1 ) / n_samples
213- if (chains > 1 )
214- var_plus <- var_plus + var(chain_mean )
215- # Geyer's initial positive sequence
216- rho_hat_t <- rep.int(0 , n_samples )
217- t <- 0
218- rho_hat_even <- 1
219- rho_hat_t [t + 1 ] <- rho_hat_even
220- rho_hat_odd <- 1 - (mean_var - mean(acov [t + 2 , ])) / var_plus
221- rho_hat_t [t + 2 ] <- rho_hat_odd
222- while (t < nrow(acov ) - 5 && ! is.nan(rho_hat_even + rho_hat_odd ) &&
223- (rho_hat_even + rho_hat_odd > 0 )) {
224- t <- t + 2
225- rho_hat_even = 1 - (mean_var - mean(acov [t + 1 , ])) / var_plus
226- rho_hat_odd = 1 - (mean_var - mean(acov [t + 2 , ])) / var_plus
227- if ((rho_hat_even + rho_hat_odd ) > = 0 ) {
228- rho_hat_t [t + 1 ] <- rho_hat_even
229- rho_hat_t [t + 2 ] <- rho_hat_odd
230- }
231- }
232- max_t <- t
233- # this is used in the improved estimate
234- if (rho_hat_even > 0 )
235- rho_hat_t [max_t + 1 ] <- rho_hat_even
236-
237- # Geyer's initial monotone sequence
238- t <- 0
239- while (t < = max_t - 4 ) {
240- t <- t + 2
241- if (rho_hat_t [t + 1 ] + rho_hat_t [t + 2 ] >
242- rho_hat_t [t - 1 ] + rho_hat_t [t ]) {
243- rho_hat_t [t + 1 ] = (rho_hat_t [t - 1 ] + rho_hat_t [t ]) / 2 ;
244- rho_hat_t [t + 2 ] = rho_hat_t [t + 1 ];
245- }
246- }
247- ess <- chains * n_samples
248- # Geyer's truncated estimate
249- # tau_hat <- -1 + 2 * sum(rho_hat_t[1:max_t])
250- # Improved estimate reduces variance in antithetic case
251- tau_hat <- - 1 + 2 * sum(rho_hat_t [1 : max_t ]) + rho_hat_t [max_t + 1 ]
252- # Safety check for negative values and with max ess equal to ess*log10(ess)
253- tau_hat <- max(tau_hat , 1 / log10(ess ))
254- ess <- ess / tau_hat
255- ess
256- }
257-
258-
259- fft_next_good_size <- function (N ) {
260- # Find the optimal next size for the FFT so that
261- # a minimum number of zeros are padded.
262- if (N < = 2 )
263- return (2 )
264- while (TRUE ) {
265- m = N
266- while ((m %% 2 ) == 0 ) m = m / 2
267- while ((m %% 3 ) == 0 ) m = m / 3
268- while ((m %% 5 ) == 0 ) m = m / 5
269- if (m < = 1 )
270- return (N )
271- N = N + 1
272- }
205+ posterior :: ess_basic(sims )
273206}
0 commit comments