Multivariate Monte Carlo standard errors for expectations with the initial sequence method of Dai and Jones (2017).
Function returns the estimate of the covariance matrix in the Markov Chain central limit theorem using initial sequence method. This method is designed to give an asymptotically conservative estimate of the Monte Carlo standard error.
mcse.initseq(x, g = NULL, adjust = FALSE, blather = FALSE)
x |
a matrix or data frame of Markov chain output. Number of rows is the Monte Carlo sample size. |
adjust |
logical; if TRUE, an adjustment is made to increase slightly the eigenvalues of the initial sequence estimator. The default is FALSE. |
g |
a function that represents features of interest. g is applied to each row of |
blather |
if TRUE, outputs under the hood information about the function. |
A list is returned with the following components,
cov |
a covariance matrix estimate using intial sequence method. |
cov.adj |
a covariance matrix estimate using adjusted initial sequence method if the input |
est |
estimate of g(x). |
nsim |
number of rows of the input |
adjust |
logical of whether an adjustment was made to the initial sequence estimator. Only if blather = TRUE |
Dai, N and Jones, G.L. (2017+) Multivariate initial sequence estimators in Markov chain Monte Carlo, Journal of Multivariate Analysis.
initseq{mcmc}, which is a different univariate initial sequence estimator.
mcse, which acts on a vector.
mcse.mat, which applies mcse to each column of a matrix or data frame.
mcse.q and mcse.q.mat, which compute standard errors for quantiles.
mcse.multi, which estimates the covariance matrix in the Markov Chain CLT using batch means or spectral variance methods.
library(mAr) p <- 3 n <- 1000 omega <- 5*diag(1,p) ## Making correlation matrix var(1) model set.seed(100) foo <- matrix(rnorm(p^2), nrow = p) foo <- foo %*% t(foo) phi <- foo / (max(eigen(foo)$values) + 1) dat <- as.matrix(mAr.sim(rep(0,p), phi, omega, N = n)) out.mcse <- mcse.initseq(x = dat) out.mcse.adj <- mcse.initseq(x = dat,adjust = TRUE) # If we are only estimating the mean of the first component, # and the second moment of the second component g <- function(x) return(c(x[1], x[2]^2)) out.g.mcse <- mcse.initseq(x = dat, g = g)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.