rfiber: Sample from the fiber of a contingency table

In dkahle/algstat: Algebraic Statistics

Sample from the fiber of a contingency table

Description

Sample from the fiber of a contingency table

Usage

rfiber(
  n,
  ...,
  A,
  b,
  x,
  method = c("sis", "walk", "hybrid"),
  dist = c("uniform", "hypergeometric"),
  parallel = FALSE,
  includeRejections = FALSE,
  thin = 1,
  format = c("mat", "vec", "tab"),
  dim
)

Arguments

n	the number of observations
...	...
A	the configuration matrix of the model defining the fiber
b	integer vector; the vector of sufficient statistics
x	integer vector; a vector in the fiber to be used as an initialization point for the "walk" method
method	"sis", "walk", or "hybrid"
dist	"uniform" or "hypergeometric"
parallel	parallelize the sampling?
includeRejections	should rejections be returned? (note: rejection tables aren't completed in the SIS procedure.)
thin	thinning argument (if using walk method)
format	format of the returned moves, "mat", "vec", or "tab"
dim	the dimensions of the table if "tab" is used, see vec2tab()

Value

a named numeric vector

Author(s)

David Kahle [e-mail], Ruriko Yoshida [e-mail]

Examples

A <- hmat(c(2,2), 1:2)
b <- rep.int(4, 4)

rfiber(10, A = A, b = b)
rfiber(10, A = A, b = b, format = "vec")
rfiber(10, A = A, b = b, format = "tab", dim = c(2, 2))

set.seed(1)
(tab <- rfiber(1, A = A, b = b, format = "tab", dim = c(2, 2))[[1]])
x <- tab2vec(tab)

## Not run:  # requires 4ti2

rfiber(100, A = A, x = x, method = "walk")

library(microbenchmark)
microbenchmark(
  rfiber(100, A = A, b = b),
  rfiber(100, A = A, x = x, method = "walk")
)

# the distribution of the samples is seemingly near-uniform
tabs <- rfiber(1e4, A, b, format = "vec", parallel = TRUE)
table(sapply(tabs, paste, collapse = " "))



tabs <- rfiber(1000, A = A, b = b, format = "vec")
unique(tabs)
lapply(unique(tabs), vec2tab, dim = c(2, 2))

table(sapply(tabs, paste, collapse = " " ))




data(politics); politics
tab2vec(politics)
b <- A %*% tab2vec(politics)
tabs <- rfiber(1000, A = A, b = b)
unique(t(tabs))
table(apply(tabs, 2, paste, collapse = " " )) # roughly uniform





# poisson regression example
J <- 5
A <- rbind(1L, 1:J)
b <- c(5, 15)




system.time(rfiber(1e4, A = A, b = b))
system.time(rfiber(1e4, A = A, b = b, parallel = TRUE))


## End(Not run)

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