Markov Basis Metropolis-Hastings Algorithm

Share:

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.

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

 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)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.