Data
The example focuses on the Escuela Centenario
polling location in the city of Iquique. This location comprises 15
polling stations, numbered 276 through 290. Most stations contain 400
registered voters, while the final two stations are smaller.
The vote matrix below records votes for four mayoral candidates, null
or blank votes, and abstention. Columns correspond to individual polling
stations.
library(MultinomFourier)
votes <- matrix(
c(
# Candidate A
37, 34, 31, 34, 37, 37, 35, 32, 40, 37, 37, 38, 36, 22, 25,
# Candidate B
39, 37, 36, 36, 26, 29, 28, 39, 39, 34, 28, 32, 35, 15, 27,
# Candidate C
99, 107, 110, 108, 127, 107, 108, 125, 89, 100, 109, 103, 98, 74, 70,
# Candidate D
83, 84, 8, 87, 70, 102, 88, 74, 82, 86, 69, 86, 86, 49, 44,
# Null and blank votes
46, 37, 42, 41, 34, 37, 45, 26, 62, 45, 48, 49, 41, 21, 21,
# Abstention
96, 101, 173, 94, 106, 88, 96, 104, 88, 98, 109, 92, 104, 57, 50
),
nrow = 6,
byrow = TRUE
)
rownames(votes) <- c(
"Candidate A",
"Candidate B",
"Candidate C",
"Candidate D",
"Null/blank",
"Abstention"
)
colnames(votes) <- as.character(276:290)
polling_station_sizes <- colSums(votes)
polling_station_sizes
#> 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290
#> 400 400 400 400 400 400 400 400 400 400 400 400 400 238 237
The first 13 stations have exactly 400 registered voters each (\(N_m = 400\)), while stations 289 and 290
have 238 and 237 voters, respectively.
Visual Inspection
Because station sizes differ slightly, proportions provide a fairer
visual comparison than raw counts. Each cell shows the percentage of
votes.
prop <- sweep(votes, 2, polling_station_sizes, "/")
heat_data <- as.data.frame(as.table(prop))
names(heat_data) <- c("Option", "Station", "Proportion")
# Cell labels: percentage with one decimal, in white on dark cells
heat_data$Label <- sprintf("%.1f", 100 * heat_data$Proportion)
dark_cell <- heat_data$Proportion > min(prop) + 0.55 * diff(range(prop))
heat_data$TextColor <- ifelse(dark_cell, "white", "black")
if (requireNamespace("ggplot2", quietly = TRUE)) {
ggplot2::ggplot(
heat_data,
ggplot2::aes(x = Station, y = Option, fill = Proportion)
) +
ggplot2::geom_tile(color = "white", linewidth = 0.25) +
ggplot2::geom_text(
ggplot2::aes(label = Label, color = TextColor),
size = 2.6
) +
ggplot2::scale_color_identity() +
ggplot2::scale_fill_distiller(
palette = "YlOrRd",
direction = 1,
labels = function(x) paste0(100 * x, "%")
) +
ggplot2::scale_y_discrete(limits = rev(rownames(votes))) +
ggplot2::labs(
x = "Polling station",
y = NULL,
fill = "Vote proportion"
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(
panel.grid = ggplot2::element_blank(),
axis.text.x = ggplot2::element_text(angle = 0, hjust = 0.5)
)
} else {
# Rows reversed so that Candidate A appears at the top
z <- t(prop[rev(rownames(prop)), , drop = FALSE])
op <- par(mar = c(5, 8, 2, 2))
image(
seq_len(nrow(z)), seq_len(ncol(z)), z,
col = hcl.colors(100, "YlOrRd", rev = TRUE),
axes = FALSE,
xlab = "Polling station",
ylab = ""
)
axis(1, at = seq_len(nrow(z)), labels = rownames(z), las = 1)
axis(2, at = seq_len(ncol(z)), labels = colnames(z), las = 1)
lab <- matrix(sprintf("%.1f", 100 * z), nrow = nrow(z))
dark <- z > min(z) + 0.55 * diff(range(z))
text(
rep(seq_len(nrow(z)), ncol(z)), rep(seq_len(ncol(z)), each = nrow(z)),
lab, cex = 0.6, col = ifelse(dark, "white", "black")
)
par(op)
}
Inspection of the heatmap highlights Station 278:
Candidate D received only 8 votes (2%), compared with 70–100 votes
(17%–25%) in neighboring stations, while abstention was markedly higher
(173 votes).
While visually salient, graphical inspection alone does not quantify
how statistically improbable this outcome is under the multinomial
model.
Leave-One-Out Reference Distribution
Let \(X_m = (X_{m1}, \ldots,
X_{mK})\) denote the vector of observed counts at station \(m\), with total sample size \(N_m = \sum_{k=1}^K X_{mk}\).
To evaluate whether station \(m\) is
compatible with the location as a whole without allowing an anomalous
station to contaminate its own reference distribution, we estimate the
null probability vector using a leave-one-out
scheme:
\[\widehat p_k^{(-m)}
=
\frac{\sum_{j \neq m} X_{jk}}{\sum_{j \neq m} N_j},
\qquad k = 1, \ldots, K.\]
Under this formulation, station \(m\) contributes zero weight to its own
reference vector \(\widehat{\mathbf{p}}^{(-m)}\).
loo_p <- sapply(
seq_len(ncol(votes)),
function(m) {
other_counts <- rowSums(votes[, -m, drop = FALSE])
other_counts / sum(other_counts)
}
)
rownames(loo_p) <- rownames(votes)
colnames(loo_p) <- colnames(votes)
# Reference distribution for station 278
round(loo_p[, "278"], 4)
#> Candidate A Candidate B Candidate C Candidate D Null/blank Abstention
#> 0.0912 0.0842 0.2700 0.2066 0.1048 0.2432
Detailed Analysis of Station 278
We model the counts at station 278 as
\[X_{278} \sim
\mathrm{Multinomial}\left(N_{278} = 400, \;
\boldsymbol{\pi}\right),\]
and test whether the vote shares \(\boldsymbol{\pi}\) at this station match
those of the rest of the location:
\[\begin{aligned}
H_0 &: \; \boldsymbol{\pi} = \widehat{\mathbf{p}}^{(-278)}, \\
H_1 &: \; \boldsymbol{\pi} \neq \widehat{\mathbf{p}}^{(-278)}.
\end{aligned}\]
We use the log-likelihood ratio statistic, which measures how far the
observed vote shares are from the reference shares:
\[G^2(x) = 2 \sum_{k=1}^{K} x_k
\ln\left(\frac{x_k}{N_{278}\, \widehat
p_k^{(-278)}}\right),\]
where \(x = (x_1, \ldots, x_K)\) are
the observed counts and terms with \(x_k =
0\) are taken as zero. The \(p\)-value is the probability, under \(H_0\), of a value of the statistic at least
as large as the observed one:
\[p = \mathbb{P}_{H_0}\left(G^2(X_{278})
\ge G^2(x)\right).\]
Other statistics are compared in a later section.
We first compare observed counts with expected counts under \(\widehat{\mathbf{p}}^{(-278)}\):
comparison_278 <- data.frame(
Option = rownames(votes),
Observed = votes[, "278"],
Expected = round(polling_station_sizes["278"] * loo_p[, "278"], 1),
Difference = round(votes[, "278"] - polling_station_sizes["278"] * loo_p[, "278"], 1),
row.names = NULL
)
knitr::kable(
comparison_278,
format = "html",
escape = FALSE,
align = "l",
caption = "Observed vs. expected counts for Polling Station 278."
)
Observed vs. expected counts for Polling Station 278.
|
Option
|
Observed
|
Expected
|
Difference
|
|
Candidate A
|
31
|
36.5
|
-5.5
|
|
Candidate B
|
36
|
33.7
|
2.3
|
|
Candidate C
|
110
|
108.0
|
2.0
|
|
Candidate D
|
8
|
82.7
|
-74.7
|
|
Null/blank
|
42
|
41.9
|
0.1
|
|
Abstention
|
173
|
97.3
|
75.7
|
Fourier inversion: pval_fourier()
The main method of MultinomFourier is
pval_fourier(). It writes the exact \(p\)-value as a Fourier series and sums it
until a stopping rule is met. The partial sums are accelerated by
extrapolation, and the summation stops as soon as either B
consecutive terms are each smaller than rel_eps / B
relative to the partial sum, or the extrapolated values vary by at most
a relative rel_eps over the last half of the terms.
Throughout this vignette we use B = 50 (the default) and
rel_eps = 1e-6, tighter than the default 1e-3:
on these data it lowers the relative error by about two orders of
magnitude, at the cost of a few thousand terms instead of a few hundred.
The result is still an approximation whose accuracy depends on how fast
the series converges; rel_eps sets how stable the sum must
be, not a bound on its error.
fit_278 <- pval_fourier(
x = votes[, "278"],
p = loo_p[, "278"],
rel_eps = 1e-6
)
fit_278
#> Multinomial Test (MultinomFourier)
#> ----------------------------------
#> Method : fourier
#> Statistic : Log-Likelihood Ratio
#> p-value : 9.2015753e-33
#> Status : 0 (Converged)
#> Gamma : 0.963835
#> Time gamma : 0.007766 seconds
#> W : 909.592
#> Terms : 9480
#> Time : 0.308082 seconds
The returned multinomfourier object is a list; its main
fields are:
fit_278$pval: computed \(p\)-value (\(\approx\) 9.2 × 10−33).fit_278$stat: test statistic used
("Log-Likelihood Ratio").fit_278$n_terms: number of Fourier series terms
evaluated (9480 here).fit_278$time: total execution time in seconds (0.31 s
here).fit_278$method: computational method used
("fourier").fit_278$status: completion status ID (0
indicates that the stopping rule was met).
The remaining fields (the exponential tilt gamma, the
period half-width W and four convergence diagnostics such
as err_rel_ext) are described in
?pval_fourier.
A \(p\)-value of order \(10^{-33}\) indicates that the pattern
observed at station 278 would be extremely unlikely under the null
model.
Exact computation: pval_exact()
The multinomial support for station 278 contains
\[\binom{N + K - 1}{K - 1} = \binom{400 +
6 - 1}{6 - 1} = \binom{405}{5} = 88,978,137,003 \approx 8.9 \times
10^{10}\]
outcomes. pval_exact() does not enumerate them one by
one. It first computes, by dynamic programming, the smallest and largest
value that the remaining part of the statistic can take for every number
of categories and voters still to be assigned. It then assigns the
counts one category at a time. When these bounds show that every
completion of a partial assignment is at least as extreme as the
observed outcome, the probability of that whole branch is added at once;
when they show that none is, the branch is discarded. Only branches that
can fall on either side of the observed value are explored further. The
result is the exact \(p\)-value, and
the computation remains feasible here even though the support is
large:
exact_fit_278 <- pval_exact(
x = votes[, "278"],
p = loo_p[, "278"]
)
exact_fit_278$pval
#> [1] 9.201581e-33
On two threads (the default) it takes 4.8 s, against 0.31 s for
pval_fourier(), and the relative difference between the two
\(p\)-values is 6.0 × 10−7.
The number of branches that must be explored, and hence the cost of
pval_exact(), still grows quickly with the number of
categories \(K\) and the sample size
\(N\), which is where
pval_fourier() is needed.
Automatic choice: pval_flexible()
pval_flexible() chooses between the two methods based on
the size of the support: below a cutoff of \(\log_{10}(\text{support}) = 10\) it
computes the exact \(p\)-value, and
above it uses Fourier inversion. The cutoff is deliberately
conservative, because the cost of pval_exact() depends on
the observed counts and is highest for extreme outcomes like this one.
For station 278, \(\log_{10}(\text{support})
\approx 10.95\), so it uses pval_fourier():
pval_flexible(
x = votes[, "278"],
p = loo_p[, "278"],
rel_eps = 1e-6
)$method
#> [1] "fourier"
Screening All Polling Stations
We now apply pval_fourier() to all 15 stations in
Escuela Centenario, each against its own leave-one-out reference vector,
and check its results against the exact \(p\)-values from pval_exact().
The data are used as they are, so the comparison covers the instance
actually being analyzed, including the extreme tail at station 278.
Because the \(p\)-values span more
than 30 orders of magnitude, we compare the two methods through the
relative error \(|\hat p_{\text{Fourier}} /
p_{\text{Exact}} - 1|\) rather than the absolute difference,
which would be uninformative for small \(p\)-values.
p_fourier_all <- sapply(colnames(votes), function(m) {
pval_fourier(votes[, m], loo_p[, m], rel_eps = 1e-6)$pval
})
p_exact_all <- sapply(colnames(votes), function(m) {
# Station 278 was already computed above
if (m == "278") return(exact_fit_278$pval)
pval_exact(votes[, m], loo_p[, m])$pval
})
rel_err_all <- abs(p_fourier_all / p_exact_all - 1)
screening <- data.frame(
Station = colnames(votes),
Voters = polling_station_sizes,
p_Fourier = fmt_p(p_fourier_all),
p_Exact = fmt_p(p_exact_all),
Rel_Error = fmt_p(rel_err_all),
row.names = NULL
)
knitr::kable(
screening,
format = "html",
escape = FALSE,
align = "l",
caption = "Leave-one-out p-values for all stations (LLR statistic): Fourier inversion vs. exact computation."
)
Leave-one-out p-values for all stations (LLR statistic): Fourier
inversion vs. exact computation.
|
Station
|
Voters
|
p_Fourier
|
p_Exact
|
Rel_Error
|
|
276
|
400
|
0.667
|
0.667
|
5.7 × 10−5
|
|
277
|
400
|
0.870
|
0.870
|
7.4 × 10−5
|
|
278
|
400
|
9.2 × 10−33
|
9.2 × 10−33
|
6.0 × 10−7
|
|
279
|
400
|
0.775
|
0.775
|
6.0 × 10−5
|
|
280
|
400
|
0.122
|
0.122
|
2.9 × 10−5
|
|
281
|
400
|
0.030
|
0.030
|
5.2 × 10−6
|
|
282
|
400
|
0.604
|
0.604
|
5.2 × 10−6
|
|
283
|
400
|
0.029
|
0.029
|
1.5 × 10−5
|
|
284
|
400
|
0.003
|
0.003
|
5.8 × 10−5
|
|
285
|
400
|
0.804
|
0.804
|
1.3 × 10−4
|
|
286
|
400
|
0.584
|
0.584
|
4.6 × 10−5
|
|
287
|
400
|
0.520
|
0.520
|
6.3 × 10−5
|
|
288
|
400
|
0.808
|
0.808
|
3.5 × 10−5
|
|
289
|
238
|
0.511
|
0.511
|
3.2 × 10−5
|
|
290
|
237
|
0.257
|
0.257
|
6.2 × 10−6
|
The largest relative error across the 15 stations is 1.3 ×
10−4, at station 285. At station 278, where \(p \approx\) 9.2 × 10−33, the
relative error is 6.0 × 10−7. Note that rel_eps
controls the stopping rule, not the error: the observed relative errors
can exceed it.
On a \(-\log_{10}(p)\) scale, where
higher values indicate stronger evidence against the reference
distribution, the two methods give the same profile:
screening_plot_data <- data.frame(
Station = rep(as.integer(colnames(votes)), 2),
Method = rep(c("Fourier (pval_fourier)", "Exact (pval_exact)"), each = ncol(votes)),
neg_log10_p = -log10(c(p_fourier_all, p_exact_all))
)
if (requireNamespace("ggplot2", quietly = TRUE)) {
ggplot2::ggplot(
screening_plot_data,
ggplot2::aes(x = Station, y = neg_log10_p, color = Method, shape = Method, linetype = Method)
) +
ggplot2::geom_line(linewidth = 0.7) +
ggplot2::geom_point(size = 2.5) +
ggplot2::scale_x_continuous(breaks = 276:290) +
ggplot2::scale_color_manual(values = c("Fourier (pval_fourier)" = "#1f77b4", "Exact (pval_exact)" = "#ff7f0e")) +
ggplot2::labs(
title = "LLR p-values by polling station at Escuela Centenario",
x = "Polling station",
y = expression(-log[10](p)),
color = "Method",
shape = "Method",
linetype = "Method"
) +
ggplot2::theme_minimal(base_size = 11) +
ggplot2::theme(legend.position = "bottom")
} else {
matplot(
276:290,
-log10(cbind(p_fourier_all, p_exact_all)),
type = "b",
pch = c(19, 17),
col = c("blue", "orange"),
lty = c(1, 2),
xlab = "Polling station",
ylab = expression(-log[10](p)),
xaxt = "n",
main = "LLR p-values by polling station at Escuela Centenario"
)
axis(1, at = 276:290, labels = 276:290)
legend("topright", legend = c("Fourier (pval_fourier)", "Exact (pval_exact)"),
col = c("blue", "orange"), pch = c(19, 17), lty = c(1, 2), bty = "n")
}
Station 278 stands out: its \(p\)-value is about 29 orders of magnitude
smaller than that of any other station (the next smallest is station
284, with \(p \approx\) 0.003).
Comparing Multinomial Test Statistics
MultinomFourier supports several classic statistics. All
of them test the same hypotheses \(H_0\) and \(H_1\) stated above for station 278; they
differ in how they order the possible outcomes from least to most
extreme, and therefore in which outcomes count as “at least as extreme”
as the observed one.
For a station with observed counts \(x =
(x_1, \ldots, x_K)\), total \(N =
\sum_k x_k\) and null probabilities \(p
= (p_1, \ldots, p_K)\) (here \(\widehat{\mathbf{p}}^{(-m)}\)), the
multinomial sample space is
\[\Omega_N = \left\{ y \in \{0, 1, \ldots,
N\}^K : \sum_{k=1}^{K} y_k = N \right\},\]
and the probability of an outcome \(y \in
\Omega_N\) under \(H_0\) is
\[P(y) = \frac{N!}{y_1! \cdots y_K!}
\prod_{k=1}^{K} p_k^{y_k}.\]
For a statistic \(T\) in which
larger values are more extreme, the \(p\)-value is
\[p = \sum_{y \in \Omega_N} P(y) \,
\mathbf{1}\left\{ T(y) \ge T(x) \right\}.\]
Probability Mass Function ("pmf").
Outcomes are ordered by their own probability under \(H_0\): an outcome is at least as extreme as
\(x\) when it is at most as probable.
The \(p\)-value is
\[p_{\text{PMF}} = \sum_{y \in \Omega_N}
P(y) \, \mathbf{1}\left\{ P(y) \le P(x) \right\}.\]
Pearson’s Chi-Square ("chi2"). Outcomes
are ordered by the squared distance between observed and expected
counts, scaled by the expected counts:
\[\mathcal{X}^2(y) = \sum_{k=1}^{K}
\frac{(y_k - N p_k)^2}{N p_k},
\qquad
p_{\mathcal{X}^2} = \sum_{y \in \Omega_N} P(y) \, \mathbf{1}\left\{
\mathcal{X}^2(y) \ge \mathcal{X}^2(x) \right\}.\]
Log-Likelihood Ratio ("llr"). The
statistic \(G^2\) defined above for
station 278, with \(p_{G^2}\) given by
the same formula with \(T = G^2\).
Power Divergence ("pd"). The family of
Cressie and Read (1984), indexed by a real parameter \(\lambda\):
\[2 I^{\lambda}(y) =
\frac{2}{\lambda(\lambda + 1)} \sum_{k=1}^{K} y_k \left[ \left(
\frac{y_k}{N p_k} \right)^{\lambda} - 1 \right],
\qquad \lambda \neq 0, -1,\]
with the values at \(\lambda = 0\)
and \(\lambda = -1\) defined as limits.
In this package, the power-divergence statistic requires \(\lambda > -1\) (values \(\lambda \le -1\) produce an error). The
family contains Pearson’s statistic (\(\lambda
= 1\)) and the log-likelihood ratio (\(\lambda \to 0\), available directly as
stat = "llr"). We use \(\lambda =
2/3\), which lies between the two.
test_specs <- list(
"PMF" = list(stat = "pmf"),
"Pearson (Chi2)" = list(stat = "chi2"),
"Log-Likelihood Ratio" = list(stat = "llr"),
"Power Divergence (2/3)"= list(stat = "pd", lambda = 2/3)
)
# Run a p-value function on station m with the statistic described by spec
run_test <- function(fun, m, spec, ...) {
args <- list(x = votes[, m], p = loo_p[, m], stat = spec$stat, ...)
if (!is.null(spec$lambda)) args$lambda <- spec$lambda
do.call(fun, args)
}
fits_by_test <- lapply(
test_specs,
function(spec) lapply(colnames(votes), function(m) run_test(pval_fourier, m, spec, rel_eps = 1e-6))
)
p_by_test <- sapply(fits_by_test, function(fits) sapply(fits, `[[`, "pval"))
status_by_test <- sapply(fits_by_test, function(fits) sapply(fits, `[[`, "status"))
rownames(p_by_test) <- rownames(status_by_test) <- colnames(votes)
# A dagger marks results where the series did not converge
p_display <- matrix(fmt_p(p_by_test), nrow = nrow(p_by_test), dimnames = dimnames(p_by_test))
p_display[status_by_test != 0] <- paste0(p_display[status_by_test != 0], " †")
comparison <- data.frame(
Station = colnames(votes),
p_display,
check.names = FALSE,
row.names = NULL
)
knitr::kable(
comparison,
format = "html",
escape = FALSE,
align = "l",
caption = "Leave-one-out p-values across different test statistics (Fourier inversion). † marks results where the series did not converge (status != 0)."
)
Leave-one-out p-values across different test statistics (Fourier
inversion). † marks results where the series did not converge (status !=
0).
|
Station
|
PMF
|
Pearson (Chi2)
|
Log-Likelihood Ratio
|
Power Divergence (2/3)
|
|
276
|
0.633
|
0.660
|
0.667
|
0.663
|
|
277
|
0.872
|
0.869
|
0.870
|
0.870
|
|
278
|
9.2 × 10−33
|
6.6 × 10−18 †
|
9.2 × 10−33
|
1.7 × 10−23 †
|
|
279
|
0.769
|
0.769
|
0.775
|
0.771
|
|
280
|
0.135
|
0.125
|
0.122
|
0.124
|
|
281
|
0.030
|
0.022 †
|
0.030
|
0.024
|
|
282
|
0.612
|
0.603
|
0.604
|
0.603
|
|
283
|
0.032
|
0.041 †
|
0.029
|
0.037
|
|
284
|
0.002
|
0.002 †
|
0.003
|
0.002
|
|
285
|
0.787
|
0.800
|
0.804
|
0.801
|
|
286
|
0.590
|
0.592
|
0.584
|
0.589
|
|
287
|
0.499
|
0.509
|
0.520
|
0.513
|
|
288
|
0.800
|
0.807
|
0.808
|
0.807
|
|
289
|
0.545
|
0.522
|
0.511
|
0.518
|
|
290
|
0.236
|
0.240
|
0.257
|
0.245
|
For typical stations (e.g. 276, 277, 282, 285), all statistics yield
comparable \(p\)-values. Cells marked †
are those where pval_fourier() reached its limit of
max_terms terms without meeting its convergence criterion
(status != 0), so the value shown there is not reliable.
This happens for Pearson’s chi-square at stations 278, 281, 283, 284,
and for the power divergence with \(\lambda =
2/3\) at station 278; PMF and LLR converge at every station. The
table below compares, for station 278, the Fourier result of each
statistic with the exact \(p\)-value
computed by pval_exact():
exact_278 <- sapply(
test_specs,
function(spec) run_test(pval_exact, "278", spec)$pval
)
fourier_278 <- lapply(fits_by_test, function(fits) fits[[which(colnames(votes) == "278")]])
comparison_278 <- data.frame(
Statistic = names(test_specs),
p_Fourier = fmt_p(sapply(fourier_278, `[[`, "pval")),
Terms = sapply(fourier_278, `[[`, "n_terms"),
Status = sapply(fourier_278, `[[`, "status"),
p_Exact = fmt_p(exact_278),
row.names = NULL
)
knitr::kable(
comparison_278,
format = "html",
escape = FALSE,
align = "l",
caption = "Station 278: Fourier inversion vs. exact p-values for each statistic."
)
Station 278: Fourier inversion vs. exact p-values for each statistic.
|
Statistic
|
p_Fourier
|
Terms
|
Status
|
p_Exact
|
|
PMF
|
9.2 × 10−33
|
7518
|
0
|
9.2 × 10−33
|
|
Pearson (Chi2)
|
6.6 × 10−18
|
10000
|
2
|
1.6 × 10−20
|
|
Log-Likelihood Ratio
|
9.2 × 10−33
|
9480
|
0
|
9.2 × 10−33
|
|
Power Divergence (2/3)
|
1.7 × 10−23
|
10000
|
2
|
8.9 × 10−24
|
At this station, for PMF and LLR the Fourier series converges (7518
and 9480 terms) and agrees with the exact value. For Pearson’s
chi-square and the power divergence it does not converge, and raising
max_terms does not fix this: with 100,000 terms the
chi-square result even becomes negative. The reason is structural.
Pearson’s statistic is \(\mathcal{X}^2 =
\sum_k x_k^2 / (N p_k) - N\) with integer counts \(x_k\), so its values lie close to a
lattice; the terms of the Fourier series then recur instead of
cancelling, and the truncated sum does not settle. The power divergence
with \(\lambda = 2/3\) behaves
similarly at this station. In practice, check the status of
every pval_fourier() result. When it is not 0, use
pval_exact() if the instance is tractable, or a statistic
such as LLR or PMF, for which the series converges.
Using the exact values for the two non-converged cells, all four
statistics flag station 278, with \(p\)-values between 9.2 × 10−33
and 1.6 × 10−20. At this station, the detection of the
anomaly does not depend on the statistic, although the size of the \(p\)-value varies by about 12 orders of
magnitude.
