stationary: Precision of stationary distribution for discrete MCMC variables

Description

Transdimensional MCMC methods include a discrete model-indicator variable \(z\) with a fixed but unknown stationary distribution with probabilities \(\pi\) (i.e., the model posterior probabiltiies). The function stationary draws posterior samples to assess the estimation uncertainty of \(\pi\).

Usage

stationary(
  z,
  N,
  labels,
  sample = 1000,
  epsilon = "1/M",
  cpu = 1,
  method = "arma",
  digits = 6,
  progress = TRUE,
  summary = TRUE
)

Value

default: a summary for the posterior distribution of the model posterior probabilities (i.e., the fixed but unknown stationary distribution of z). If summary=FALSE, posterior samples for \(pi\) are returned.

Arguments

z

MCMC output for the discrete indicator variable with numerical, character, or factor labels (can also be a mcmc.list or a matrix with one MCMC chain per column).

N

the observed transitions matrix (if supplied, z is ignored). A quadratic matrix with sampled transition frequencies (N[i,j] = number of switches from z[t]=i to z[t+1]=j).

labels

optional: vector of labels for complete set of models (e.g., models not sampled in the chain z). If epsilon=0, this does not affect inferences due to the improper Dirichlet(0,..,0) prior.

sample

number of posterior samples to be drawn for the stationary distribution \(\pi\).

epsilon

prior parameter for the rows of the estimated transition matrix \(P\): \(P[i,]\) ~ Dirichlet\((\epsilon, ..., \epsilon)\). The default epsilon="1/M" (with M = number of sampled models) provides estimates close to the i.i.d. estimates and is numerically stable. The alternative epsilon=0 minimizes the impact of the prior and renders non-sampled models irrelevant. If method="iid" (ignores dependencies), a Dirichlet prior is assumed on the stationary distribution \(\pi\) instead of the rows of the transition matrix \(P\).

cpu

number of CPUs used for parallel sampling. Will only speed up computations for large numbers of models (i.e., for large transition matrices).

method

how to compute eigenvectors:

"arma" (default): Uses RcppArmadillo::eig_gen.
"base": Uses base::eigen, which might be more stable, but also much slower than "arma" for small transition matrices.
"eigen": Uses package RcppEigen::EigenSolver
"armas": Uses sparse matrices with RcppArmadillo::eigs_gen, which can be faster for very large number of models if epsilon=0 (might be numerically unstable).
"iid": Assumes i.i.d. sampling of the model indicator variable z. This is only implemented as a benchmark, because results cannot be trusted if the samples z are correlated (which is usually the case for transdimensional MCMC output)

digits

number of digits that are used for checking whether the first eigenvalue is equal to 1 (any difference must be due to low numerical precision)

progress

whether to show a progress bar (not functional for cpu>1)

summary

whether the output should be summarized. If FALSE, posterior samples of the stationary probabilities are returned.

Details

The method draws independent posterior samples of the transition matrix \(P\) for the discrete-valued indicator variable z (usually, a sequence of sampled models). For each row of the transition matrix, a Dirichlet\((\epsilon,...,\epsilon)\) prior is assumed, resulting in a conjugate Dirichlet posterior. For each sample, the eigenvector with eigenvalue 1 is computed and normalized. These (independent) posterior samples can be used to assess the estimation uncertainty in the stationary distribution \(pi\) of interest (e.g., the model posterior probabilities) and to estimate the effective sample size (see summary.stationary).

Examples

Run this code

# data-generating transition matrix
P <- matrix(c(.1,.5,.4,
              0, .5,.5,
              .9,.1,0), ncol = 3, byrow=TRUE)

# input: sequence of sampled models
z <- rmarkov(500, P)
stationary(z)

# input: transition frequencies
N <- transitions(z)
samples <- stationary(N = N, summary = FALSE)

# summaries:
best_models(samples, k = 3)
summary(samples)