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

 1 2 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.

a list

David Kahle

References

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

Examples

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 ## 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.