stationary: Precision of stationary distribution for discrete MCMC...

In danheck/MCMCprecision: Precision of Discrete Parameters in Transdimensional MCMC

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
)

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).

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.

See Also

best_models, summary.stationary

Examples

# 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)

danheck/MCMCprecision documentation built on July 3, 2024, 8:51 p.m.