Moments from cumulants

moments_from_cumulants.md

Moments from cumulants

IH moments are difficult to calculate explicitly, but the cumulants are easy. We can use the functional definition of the cumulant generating function to 'invert' back to the moments.

By definition, the moment generating function $M_X(t)$ of a random variable $X$ is: $$M_{X}(t) \equiv \mathop{\mathbb{E}}e^{tX}$$ and the cumulant generating function $C_X(t)$ is:

$$ C_{X}(t) \equiv \log M_{X}(t) = \log \mathop{\mathbb{E}}e^{tX} $$

This implies that:

$$ M'_X = C'_X M_X$$

If $X$ and $Y$ are independent, $M_{X+Y} = M_X M_Y$ so $C_{X+Y} = C_X + C_Y$ by taking $\log$ of both sides. Thus the cumulants of a sum of independent variables is the sum of their cumulants:

$$ \kappa^{X+Y}_n = \kappa^X_n + \kappa^{Y}_n $$

Now expanding to their exponential series forms ($\mu_0=1$ and $\kappa_0=0$)

$$ C_X = \sum_{m=0}^\infty \kappa_m \frac{x^m}{m!} $$

$$ M_X = \sum_{m=0}^\infty \mu_m \frac{x^m}{m!} $$

$$ C'_X = \sum_{m=0}^\infty \kappa_{m+1} \frac{x^m}{m!} $$

$$ M'_X = \sum_{m=0}^\infty \mu_{m+1} \frac{x^m}{m!} $$

Taking the Cauchy product of $C'_X$ and $M_X$ as exponential generating functions and collecting like terms of $x^n$, we get the umbral convolution formula for the moments:

$$\mu_{n+1} = \sum_{k=0}^{n}\mu_{n-k}\kappa_{k+1}\binom{n}{k}$$ $$\mu_{n+1} = \kappa_{n+1}+\sum_{k=0}^{n-1}\mu_{n-k}\kappa_{k+1}\binom{n}{k}$$ letting $m = n+1$ and $j = k+1$ $$\mu_{m} = \kappa_{m}+\sum_{k=0}^{m-2}\mu_{m-1-k}\kappa_{k+1}\binom{m-1}{k}$$ $$\mu_{m} = \kappa_{m}+\sum_{j=1}^{m-1}\mu_{m-j}\kappa_{j}\binom{m-1}{j-1}$$ $$\mu_{m} = \sum_{j=1}^{m}\mu_{m-j}\kappa_{j}\binom{m-1}{j-1}$$

Relabeling:

$$\mu_{n} = \sum_{k=1}^{n}\mu_{n-k}\kappa_{k}\binom{n-1}{k-1}$$

$$$$

For the cumulants of the $\mathrm{Uniform}(0, 1)$ distribution, we have $\kappa^{U}_{n} = \frac{B_{n}}{n}$ where $B_n$ = Bernoulli(n) with Bernoulli(1) := 1/2 and Bernoulli(2j+1) = 0 for j > 0.

Therefore the cumulants of the $\mathrm{IH}(m)$ distribution, as the sum of $m$ standard uniform random variables, are:

$$\kappa^{IH(m)}_n = \frac{mB_n}{n}$$

Python

Note

In this code, n is the moment order and m is parameter of the IH distribution (swapped from above)

Scipy's scipy.special.bernoulli uses the convention that $B_1$ is $-1/2$ so we adjust for that in the code below.

def _munp(self, n, m):
    cumulants = sc.special.bernoulli(n)
    cumulants[1:] /= np.arange(1, n+1)
    cumulants[0] = 0
    cumulants[1] = 1/2
    cumulants *= m
    moments = np.zeros(n+1)
    moments[0] = 1

    for i in range(1, n+1):
        moments[i] = np.dot(cumulants[1:i+1]*moments[i-1::-1], sc.special.binom(i-1,np.arange(0, i)))

    return moments[n]