In pierrejacob/montecarlodsm: Gibbs sampler for Dempster-Shafer inference in Categorical distributions

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "README-"
)

dempsterpolytope

This package implements a Gibbs sampler for Dempster's inference approach for Categorical distributions using the "simplicial model" or "structure of the second kind". This package is a companion to an article entitled A Gibbs sampler for a class of random convex polytopes by Pierre E. Jacob, Ruobin Gong, Paul T. Edlefsen, Arthur P. Dempster, available at https://arxiv.org/abs/1910.11953v3, and due to appear with a discussion in the Journal of the American Statistical Association, Theory & Methods, in 2021 (see https://www.tandfonline.com/doi/abs/10.1080/01621459.2021.1881523?journalCode=uasa20).

The sampling problem addressed by this package dates back to these articles:

• Dempster, A.P., 1966. New methods for reasoning towards posterior distributions based on sample data. The Annals of Mathematical Statistics, 37(2), pp.355-374.

• Dempster, A.P., 1972. A class of random convex polytopes. The Annals of Mathematical Statistics, pp.260-272.

Important: this is not a general-purpose statistical software. This is just a collection of scripts intended to reproduce figures and tables of a paper. Use at your own risk!

The folder inst/reproduce/ contains the scripts to reproduce the figures of the article. The folder inst/tests/ contains internal checks, and vignettes/ contains tutorials on how to use the package's main functions, which at the moment are just R scripts and not proper "R vignettes".

Installation

The package can be installed from R via:

# install.packages("devtools")
devtools::install_github("pierrejacob/dempsterpolytope")

It depends on the packages Rcpp, RcppEigen, igraph, rcdd, dplyr, lpSolveAPI which can be installed via:

install.packages(c("Rcpp", "RcppEigen", "igraph", "rcdd", "dplyr", "lpSolveAPI"))

Additionally you might want to install other packages, to help with parallel computation:

install.packages(c("doParallel", "doRNG"))

and to help with manipulating results and plotting:

install.packages(c("tidyr", "ggplot2"))

although these packages are not strictly required.

Usage

The "Dempster-Shafer" approach to inference leads, in the case of Categorical distributions, to random convex polytopes. The Gibbs sampler implemented in this package generates such polytopes. The following code shows a random polytope obtained after some iterations of the Gibbs sampler, for the data set (10, 4, 7), meaning 10 observations in category 1, 4 in category 2, 7 in category 3.

library(dempsterpolytope)
set.seed(1)
# count data
counts <- c(10, 4, 7)
# number of MCMC iterations
niterations <- 100
# run Gibbs sampler
gibbs_results <- gibbs_sampler(niterations, counts)
# obtain a K x K matrix representing a convex polytope in the simplex
# by taking the terminal iteration of the Gibbs chain
eta <- gibbs_results$etas[niterations,,]
# There are multiple ways of representing a polytope.
# The H-representation, or "half-plane" representation,
# means the set is represented as points x such that 'constr' times x <= 'rhs' 
Hrep <- eta_linearconstraints(eta)
print(Hrep$constr)
print(Hrep$dir)
print(Hrep$rhs)
# In small dimensions we can also obtain the list of vertices
# corresponding to this polytope, and then the polytope itself
# is the convex hull of these vertices
eta_vertices <- etas_vertices(eta)
print(eta_vertices)
# next we can view the K-simplex as a triangle, with K = 3 here
# and the feasible polytopes as polygons within the triangle
gs <- set_custom_theme()
g <- create_plot_triangle(gs)
cvxpolytope_cartesian.df <- data.frame()
for (iter in 51:100){
  vertices <- etas_vertices(gibbs_results$etas[iter,,])
  cvx_cartesian <- t(apply(vertices, 1, function(row) barycentric2cartesian(row, gs$v_cartesian)))
  average_ <- colMeans(cvx_cartesian)
  o_ <- order(apply(sweep(cvx_cartesian, 2, average_, "-"), 1, function(v) atan2(v[2], v[1])))
  cvx_cartesian <- cvx_cartesian[o_,]
  cvxpolytope_cartesian.df <- rbind(cvxpolytope_cartesian.df, data.frame(cvx_cartesian, iter = iter))
}
gpolytopes <- g + geom_polygon(data = cvxpolytope_cartesian.df, aes(x = X1, y = X2, group = iter), size = 0.25, alpha = .2, fill = 'black', colour = 'black')
gpolytopes

The generated polytopes are black polygons overlaid in the simplex.

pierrejacob/montecarlodsm documentation built on June 16, 2021, 1:06 p.m.