= Solution
Assuming independent digits under <Benford law>, put $N=\sum_i y_i$ and $p_i=\log_{10}(1+1/i)$. The count vector has a <multinomial distribution>, with expected counts $E_i=Np_i$. Suitable <predictive discrepancy statistics> include the <Pearson chi-squared statistic> and <multinomial deviance>:
$$
\boxed{T_P=\sum_{i=1}^9\frac{(y_i-E_i)^2}{E_i},\qquad
T_G=2\sum_{i:y_i>0}y_i\log(y_i/E_i).}
$$
The zero-count terms of $T_G$ have limiting value zero. A <Monte Carlo method> gives a direct null comparison even when expected counts are small.
For example, an original R implementation is:
``
benford_check <- function(y, B = 9999L) {
p <- log10(1 + 1/(1:9))
n <- sum(y)
expected <- n*p
observed <- sum((y - expected)^2/expected)
replicas <- rmultinom(B, size = n, prob = p)
simulated <- colSums((replicas - expected)^2/expected)
(1 + sum(simulated >= observed))/(B + 1)
}
``
Each column returned by `rmultinom` is a replicated count vector. R recycles the nine expected counts down each column. The add-one ratio is a <Monte Carlo test> estimate. In WinBUGS, alternatively generate a replicated `dmulti` vector using fixed $p_i,N$, compute its discrepancy and monitor exceedance of the observed value. The null does not estimate unknown digit probabilities. Dependence or selection in the accounts would require an appropriate simulation model.
Back to article page
