mcse.initseq {mcmcse}

Multivariate Monte Carlo standard errors for expectations with the initial sequence method of Dai and Jones (2017)

Description

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.

Usage

mcse.initseq(x, g = NULL, adjust = FALSE, blather = FALSE)

Arguments

x

A matrix or data frame of Markov chain output. Number of rows is the Monte Carlo sample size.

g

A function that represents features of interest. g is applied to each row of x and thus g should take a vector input only. If g is NULL, g is set to be identity, which is estimation of the mean of the target density.

adjust

Logical; if TRUE, an adjustment is made to increase slightly the eigenvalues of the initial sequence estimator. The default is FALSE.

blather

if TRUE, outputs under the hood information about the function.

Value

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 adjust=TRUE.

eigen_values

eigen values of the estimate cov.

method

method used to calculate matrix cov.

est

estimate of g(x).

nsim

number of rows of the input x. Only if blather = TRUE.

Adjustment_Used

logical of whether an adjustment was made to the initial sequence estimator. Only if blather = TRUE.

References

Dai, N and Jones, G.L. (2017) Multivariate initial sequence estimators in Markov chain Monte Carlo, Journal of Multivariate Analysis, 159, 184-199.

See Also

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.

Examples

## Bivariate Normal with mean (mu1, mu2) and covariance sigma
n <- 1e3
mu <- c(2, 50)
sigma <- matrix(c(1, 0.5, 0.5, 1), nrow = 2)
out <- BVN_Gibbs(n, mu, sigma)

out.mcse <- mcse.initseq(x = out)
out.mcse.adj <- mcse.initseq(x = out, 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 = out, g = g)