Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 3 / b

A numerical power series evaluation needs a stopping rule. For

r = ∣ z ∣ < 1

and integer

a \geq 0

, the absolute remainder after

N

terms satisfies the polylogarithm series remainder bound

Li_{a} (z) - k = 1 \sum N \frac{z ^{k}}{k ^{a}} \leq \frac{r ^{N + 1}}{( N + 1 ) ^{a} ( 1 - r )} .

(1)

Indeed

k^{- a} \leq (N + 1)^{- a}

for

k > N

, and the remaining powers sum as a geometric series. This also covers negative

z

. The code sums terms in vectors of at most 4096 entries, checking the bound after each block; there is no loop over individual terms.

polylog <- function(a, z, tol = 1e-10, max_terms = 1000000) {
  scalar <- function(v) {
    is.numeric(v) && !is.complex(v) && length(v) == 1L &&
      is.finite(v)
  }
  if (!scalar(a) || a < 0 || a != floor(a))
    stop("a must be a nonnegative integer")
  if (!scalar(z) || abs(z) >= 1)
    stop("z must be real with abs(z) < 1")
  if (!scalar(tol) || tol <= 0 || tol >= 1)
    stop("tol must lie strictly between 0 and 1")
  if (!scalar(max_terms) || max_terms < 1 ||
      max_terms != floor(max_terms))
    stop("max_terms must be a positive integer")
  if (z == 0) return(0)
  r <- abs(z)
  N <- 0
  ans <- 0
  repeat {
    if (N >= max_terms)
      stop("requested truncation tolerance not reached")
    B <- min(4096, max_terms - N)
    k <- seq.int(N + 1, N + B)
    ans <- ans + sum(z^k * exp(-a * log(k)))
    N <- N + B
    log_bound <- (N + 1) * log(r) - a * log(N + 1) - log1p(-r)
    if (log_bound <= log(tol)) return(ans)
  }
}

This accepts every positive integer order in the question, and also order zero if

0 \in N

is the convention adopted. The bound controls absolute truncation error; floating-point rounding is a separate source of error. Close to

∣ z ∣ = 1

, convergence can be slow. If the resource limit is reached, the function reports failure instead of silently claiming the requested accuracy. Useful exact checks are

Li_{1} (z) = - lo g (1 - z)

and

Li_{0} (z) = z / (1 - z)

ii

 0  0

Solution

 0  0

Let

m = E X

and

v = 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:

Z_{m} = \frac{n ( X _{n} - m )}{v} \Rightarrow N (0, 1) .

(1)

For the second check, the correct asymptotic variance of the sample variance involves the fourth central moment

μ_{4} = E [(X - m)^{4}]

Z_{v} = \frac{n ( S _{n}^{2} - v )}{μ _{4} - v ^{2}} \Rightarrow N (0, 1) .

(2)

The denominator is not

2 v^{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

μ_{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} := E [X^{j}] = - \frac{j ! Li _{j + 1} ( 1 - p )}{β ^{j} l o g p}, j \geq 1.

(3)

Thus

v = m_{2} - m_{1}^{2}

and

μ_{4} = m_{4} - 4 m_{1} m_{3} + 6 m_{1}^{2} m_{2} - 3 m_{1}^{4}

. The following R 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