metropolis: Markov Basis Metropolis-Hastings Algorithm

In algstat: Algebraic statistics in R

Description

Given a starting table (as a vector) and a loglinear model matrix A, compute the Markov basis of A with 4ti2 and then run the Metropolis-Hastings algorithm starting with the starting table.

Usage

metropolis(init, moves, iter = 1000, burn = 1000, thin = 10,
  engine = c("Cpp", "R"))

Arguments

init	the initial step
moves	the markov basis (the negatives will be added). see ?markov
iter	number of chain iterations
burn	burn-in
thin	thinning
engine	C++ or R? (C++ yields roughly a 20-25x speedup)

Details

See Algorithm 1.1.13 in LAS, the reference below.

Value

a list

Author(s)

David Kahle

References

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

Examples

## Not run: 



data(handy)

exp   <- loglin(handy, as.list(1:2), fit = TRUE)$fit
e <- unname(tab2vec(exp))
h <- t(t(unname(tab2vec(handy))))
chisq <- algstat:::computeChisqsCpp(h, e)

out <- hierarchical(~ Gender + Handedness, data = handy)
chisqs <- algstat:::computeChisqsCpp(out$steps, e)

mean(chisqs >= chisq)
fisher.test(handy)$p.value





A <- hmat(c(2,2), as.list(1:2))
moves <- markov(A)
outC <- metropolis(tab2vec(handy), moves, 1e4, engine = "Cpp")
str(outC)
outR <- metropolis(tab2vec(handy), moves, 1e4, engine = "R", thin = 20)
str(outR)

# showSteps(out$steps)


library(microbenchmark)
microbenchmark(
  metropolis(tab2vec(handy), moves, engine = "Cpp"),
  metropolis(tab2vec(handy), moves, engine = "R")
)

# cpp ~ 20-25x faster













showSteps <- function(steps){
  apply(steps, 2, function(x){
    x <- format(x)
    tab <- vec2tab(x, dim(handy))
    message(
      paste(
        apply(tab, 1, paste, collapse = " "),
        collapse = " "
      )
    )
    message("
", appendLF = F)
  })
  invisible()
}
# showSteps(out$steps)









## End(Not run)

algstat documentation built on May 29, 2017, 10:34 p.m.