Solution (source code)

= Solution
{c}

Let $m=\mathbb EX$ and $v=\operatorname{Var}(X)$, evaluated from the given <polylogarithm> formulas. Generate a large sample of <independent and identically distributed random variables> with `r_sample_f`, and compute its <sample mean> $\overline X_n$ and <sample variance> $S_n^2$ using `mean` and `var`. The <central limit theorem> gives the first standardized check:
$$
\boxed{Z_m=\frac{\sqrt n(\overline X_n-m)}{\sqrt v}\ \Rightarrow\ N(0,1).}
$$
For the second check, the correct asymptotic <variance> of the <sample variance> involves the fourth <central moment> $\mu_4=\mathbb E[(X-m)^4]$:
$$
\boxed{Z_v=\frac{\sqrt n(S_n^2-v)}{\sqrt{\mu_4-v^2}}\ \Rightarrow\ N(0,1).}
$$
The denominator is not $\sqrt{2v^2}$ unless the <probability distribution> has the requisite fourth <central moment>, as does the <normal distribution>. Here the <probability density function> decays exponentially, so all required <moments> are finite.

To compute $\mu_4$, expand the denominator of the <probability density function> as a <geometric series> and integrate each nonnegative term. The <gamma function> integral gives the <raw moments>
$$
m_j:=\mathbb E[X^j]
=-\frac{j!\operatorname{Li}_{j+1}(1-p)}{\beta^j\log p},\qquad j\geq1.
$$
Thus $v=m_2-m_1^2$ and $\mu_4=m_4-4m_1m_3+6m_1^2m_2-3m_1^4$. The following <R (programming language)> code implements both checks for chosen interior parameters and large `n`:

``
p <- 0.4
beta <- 1.3
n <- 100000
raw_m <- vapply(1:4, function(j) {
  -factorial(j) * polylog(j + 1, 1 - p) / (beta^j * log(p))
}, numeric(1))
m <- raw_m[1]
v <- raw_m[2] - m^2
mu4 <- raw_m[4] - 4 * m * raw_m[3] +
  6 * m^2 * raw_m[2] - 3 * m^4
x <- r_sample_f(n, p, beta)
c(mean_check = sqrt(n) * (mean(x) - m) / sqrt(v),
  variance_check = sqrt(n) * (var(x) - v) / sqrt(mu4 - v^2))
``

For each check separately, an absolute value below $1.96$ is consistent with its approximate $95\%$ <confidence interval>. Repeat with many <independent> samples: the two standardized statistics should each resemble a <standard normal distribution>, with about $5\%$ outside that range. Also vary $p$, $\beta$ and `n`; a systematic discrepancy that grows with $\sqrt n$ suggests a bias. The two checks from the same sample are generally dependent, so passing both is not itself a joint $95\%$ statement, nor can a finite test prove that a generator is correct.

The second <central limit theorem> follows from the <central limit theorem for the sample variance>: $\sqrt n(S_n^2-v)$ differs by a quantity tending to zero in <probability> from $n^{-1/2}\sum_i[(X_i-m)^2-v]$.