metropolis: The Metropolis Algorithm

In dkahle/algstat: Algebraic Statistics

The Metropolis Algorithm

Description

Given a starting table (as a vector) and a collection of moves, run the Metropolis-Hastings algorithm starting with the starting table.

Usage

metropolis(
  init,
  moves,
  iter = 1000L,
  burn = 0L,
  thin = 1L,
  dist = c("hypergeometric", "uniform"),
  engine = c("C++", "R")
)

rawMetropolis(init, moves, iter = 1000, dist = "hypergeometric")

Arguments

init	the initial step
moves	the moves to be used (the negatives will be added); they are arranged as the columns of a matrix.
iter	number of chain iterations
burn	burn-in
thin	thinning
dist	steady-state distribution; "hypergeometric" (default) or "uniform"
engine	"C++" or "R"? (C++ is significantly faster)

Details

See Algorithm 1.1.13 in LAS, the reference below.

Value

a list containing named elements

• steps: an integer matrix whose columns represent individual samples from the mcmc.

• moves: the moves supplied.

• accept_prob: the empirical transition probability of the moves, including the thinned moves.

• accept_count: the numerator of accept_prob.

• total_count: the numerator of total_prob.

Author(s)

David Kahle

References

Drton, M., B. Sturmfels, and S. Sullivant (2009). Lectures on Algebraic Statistics, Basel: Birkhauser Verlag AG.

Examples

## basic use
############################################################

# move up and down integer points on the line y = 100 - x
# sampling from the hypergeometric distribution
# note: negative moves are added internally
init <- c(10L, 90L)
moves <- matrix(c(1,-1), ncol = 1)

# it helps running each of these lines several times to get a feel for things
metropolis(init, moves, iter = 10, burn = 0, thin = 1)
metropolis(init, moves, iter = 10, burn = 0, thin = 1, dist = "uniform")
metropolis(init, moves, iter = 10, burn = 0, thin = 1, engine = "R")
metropolis(init, moves, iter = 10, burn = 0, thin = 1, engine = "R", dist = "uniform")

# a bigger simulation
iter <- 1e4
out <- metropolis(init, moves, iter = iter, burn = 0)
hist(out$steps[1,], breaks = 100, col = "gray20")

# view convergence through trace plot
plot(1:iter, out$steps[1,1:iter])

# look at autocorrelation
acf(out$steps[1,])


## thinning to reduce autocorrelation with the thin argument
############################################################

out <- metropolis(init, moves, iter = iter, thin = 200)
acf(out$steps[1,])
hist(out$steps[1,], breaks = 100, col = "gray20")



## burn in with the burn argument
############################################################

set.seed(1L)
metropolis(init, moves, iter = 10, burn = 0, thin = 1, engine = "R")$steps

set.seed(1L)
metropolis(init, moves, iter = 10, burn = 0, thin = 1, engine = "R")$steps

set.seed(1L)
metropolis(init, moves, iter =  5, burn = 5, thin = 1, engine = "R")$steps


# to do, align these:
set.seed(1L)
metropolis(init, moves, iter = 10, burn = 0, thin = 1)$steps

set.seed(1L)
metropolis(init, moves, iter = 10, burn = 0, thin = 1)$steps

set.seed(1L)
metropolis(init, moves, iter =  5, burn = 5, thin = 1)$steps

dkahle/algstat documentation built on May 23, 2023, 12:29 a.m.