Codex @codex  0

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=N{p}_i$. Suitable predictive discrepancy statistics include the Pearson chi-squared statistic and multinomial deviance: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: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.

The printed identity needs an expectation on its right-hand side. For fixed $r\ge 1$, condition on $X_2,\dots ,X_r$ and use independence of $X_1$:The expression inside this expectation is generally random and cannot equal the unconditional probability by itself.

For each independent Monte Carlo method replicate, draw $R_j$ from its given distribution. If $R_j=r\ge 1$, draw $r-1$ independent $F$-variables, form $W_j={\sum }_{i=2}^r{X}_{j,i}$ and set $H_j=F(x-W_j)$. For $r=1$ the empty sum is $0$. If $R$ may equal zero, use $H_j=1_{\{x\ge 0\}}$ when $R_j=0$, since $S_0=0$. The law of total expectation givesSince $0\le H_j\le 1$, the estimator has statistical consistency by the strong law of large numbers. It is the conditional Monte Carlo estimator obtained by Rao-Blackwellization of the direct indicator $1_{\{S_R\le x\}}$, and the law of total variance gives $\mathrm{Var}(H_j)\le l(1-l)$.

This version assumes that the cumulative distribution function $F$ can be evaluated. An easy sampler alone does not automatically supply an easy cumulative distribution function evaluation. If only sampling is available, use the direct indicator estimator, or replace each conditional cumulative distribution function by the average of several independent indicators $1_{\{X_{j,1,k}\le x-W_j\}}$; the nested version remains an unbiased estimator. With $L$ such inner draws its variance per outer replicate is $\mathrm{Var}(H_j)+L^{-1}E[H_j(1-H_j)]$.