algstat (version 0.0.2)

metropolis: Markov Basis Metropolis-Hastings Algorithm

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)

Value

a list

Details

See Algorithm 1.1.13 in LAS, the reference below.

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)