MultinomFourier

MultinomFourier is an R package that computes p-values of multinomial goodness-of-fit tests,

\[ p = P\{T(Y) \ge T(x)\}, \qquad Y \sim \mathrm{Multinomial}(N, p_1, \ldots, p_K), \]

for additive test statistics \(T\): the log-likelihood ratio \(G^2\), Pearson’s \(\mathcal{X}^2\), the Cressie–Read power divergence, and the probability mass function of the multinomial (exact multinomial test). It is aimed at small p-values, far in the tail, and at instances whose sample space is too large to sum over.

The package has three functions:

Function	Method
`pval_fourier()`	Writes the p-value as a Fourier series whose terms are values of the moment-generating function of the statistic, and sums it until an adaptive stopping rule is met. Its cost grows slowly with \(N\) and \(K\).
`pval_exact()`	Exact p-value by branch pruning with bounds from a dynamic program, without enumerating the sample space. Feasible for a moderate number of categories and trials.
`pval_flexible()`	Chooses between the two according to the size of the sample space.

Installation

# install.packages("remotes")
remotes::install_github("PabRvss/MultinomFourier", build_vignettes = TRUE)

The package contains C code, so it needs the standard tools for building R packages from source: Xcode Command Line Tools on macOS, Rtools on Windows, or a C compiler on Linux.

Example

library(MultinomFourier)

# Note: Reported p-values in comments are illustrative and may vary slightly
# across versions or due to changes in numerical criteria.
p <- c(0.1, 0.2, 0.3, 0.4)   # null probabilities
x <- c(35, 30, 60, 75)       # observed counts, N = 200

fit <- pval_fourier(x, p, stat = "llr")
fit$pval     # 0.00688
fit$status   # 0: the stopping rule was met

pval_exact(x, p, stat = "llr")$pval   # 0.00684

# Power divergence: lambda is required (must be strictly greater than -1)
pval_fourier(x, p, stat = "pd", lambda = 2/3)$pval

All three functions return an object of class "multinomfourier", a list with the p-value, the status, the running time and, for the Fourier method, the details of the series (see ?pval_fourier).

Things to keep in mind

Check status. The Fourier series is truncated by a stopping rule, so pval_fourier() returns an approximation. status = 0 means the stopping rule was met; any other value means the result is not reliable.
The statistic matters. For the log-likelihood ratio ("llr") and the probability mass function ("pmf") the series typically converges in a few hundred to a few thousand terms. For Pearson’s \(\mathcal{X}^2\) ("chi2"), and often for the power divergence unless \(\lambda\) is close to 0, it may not converge at all. In that case use pval_exact() when the instance is small enough.
rel_eps is not an error bound. It sets how stable the partial sums must be before stopping (default 1e-3); the actual relative error can be larger.

Documentation

Reference manual: MultinomFourier.pdf in this repository, or help(package = "MultinomFourier").
Worked example with election data: vignette("polling_tables_multfourier", package = "MultinomFourier").

Reference

Subiabre, F. and Thraves, C. Efficient computation of p-values for multinomial tests. Manuscript.

Authors

Charles Thraves, Pablo Rivas (maintainer), Felipe Subiabre.