p-Values of Multinomial Goodness-of-Fit Tests by Fourier Series Inversion

Computes p-values of multinomial goodness-of-fit tests based on an additive statistic: the log-likelihood ratio, Pearson's chi-squared, the Cressie-Read power divergence, and the exact multinomial test. The main function, pval_fourier(), writes the p-value as a Fourier series in the moment-generating function of the statistic and sums it until an adaptive stopping rule is met, which makes small p-values of large instances computationally feasible. The function pval_exact() computes the exact p-value by branch pruning with bounds from a dynamic program, and pval_flexible() chooses between the two according to the size of the instance.

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:

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.