Dr. John Kolassa

Revision as of 19:55, 4 September 2008

It is convenient notationally to adopt Einstein's summation convention, so \(\xi_r X^r\) denotes the linear combination \(\xi_1 X^1 + \cdots + \xi_k X^k\), the square of the linear combination is \((\xi_r X^r)^2 = \xi_r\xi_s X^r X^s\) a sum of \(k^2\) terms, and so on for higher powers. The Taylor expansion of the moment generating function \(M(\xi) = E(\exp(\xi_r X^r)\) is \( M(\xi) = 1 + \xi_r \kappa^r + \textstyle{\frac1{2!}} \xi_r\xi_s \kappa^{rs} + \textstyle{\frac1{3!}} \xi_r\xi_s \xi_t \kappa^{r s t} +\cdots. \) The cumulants are defined as the coefficients \(\kappa^{r,s}, \kappa^{r,s,t},\ldots\) in the Taylor expansion \( \log M(\xi) = \xi_r \kappa^r + \textstyle{\frac1{2!}} \xi_r\xi_s \kappa^{r,s} + \textstyle{\frac1{3!}} \xi_r\xi_s \xi_t \kappa^{r,s,t} +\cdots. \) This notation does not distinguish first-order moments from first-order cumulants, but commas separating the superscripts serve to distinguish higher-order cumulants from moments.

Comparison of coefficients reveals that the each moment \(\kappa^{rs}, \kappa^{r s t},\ldots\) is a sum over partitions of the superscripts, each term in the sum being a product of cumulants: \begin{eqnarray*} \kappa^{rs}&=&\kappa^{r,s} + \kappa^r\kappa^s\\ \kappa^{r s t}&=&\kappa^{r,s,t} + \kappa^{r,s}\kappa^t + \kappa^{r,t}\kappa^s + \kappa^{s,t}\kappa^r + \kappa^r\kappa^s\kappa^t\\ &=&\kappa^{r,s,t} + \kappa^{r,s}\kappa^t[3] + \kappa^r\kappa^s\kappa^t\\ \kappa^{r s t u}&=&\kappa^{r,s,t,u} + \kappa^{r,s,t}\kappa^u[4] + \kappa^{r,s}\kappa^{t,u}[3] + \kappa^{r,s}\kappa^t\kappa^u[6] + \kappa^r\kappa^s\kappa^t\kappa^u. \end{eqnarray*} Each parenthetical number indicates a sum over distinct partitions having the same block sizes, so the fourth-order moment is a sum of 15 distinct cumulant products. In the reverse direction, each cumulant is also a sum over partitions of the indices. Each term in the sum is a product of moments, but with coefficient \((-1)^{\nu-1} (\nu-1)!\) where \(\nu\) is the number of blocks: \begin{eqnarray*} \kappa^{r,s} &=& \kappa^{rs} - \kappa^r\kappa^s\\ \kappa^{r,s,t} &=& \kappa^{r s t} - \kappa^{rs}\kappa^t[3] + 2 \kappa^r\kappa^s\kappa^t\\ \kappa^{r,s,t,u} &=& \kappa^{r s t u} - \kappa^{r s t}\kappa^u[4] - \kappa^{rs}\kappa^{t u}[3] + 2 \kappa^{rs}\kappa^t\kappa^u[6] - 6 \kappa^r\kappa^s\kappa^t\kappa^u \end{eqnarray*}

These relationships are an instance of M\"obius inversion on the partition lattice.

Partition notation serves one additional purpose. It establishes moments and cumulants as special cases of generalized cumulants, which includes objects of the type \(\kappa^{r,st} = {\rm cov}(X^r, X^s X^t)\), \(\kappa^{rs, t u} = {\rm cov}(X^r X^s, X^t X^u)\), and \(\kappa^{rs, t, u}\) with incompletely partitioned indices. These objects arise very naturally in statistical work involving asymptotic approximation of distributions. They are intermediate between moments and cumulants, and have characteristics of both.

Every generalized cumulant can be expressed as a sum of certain products of ordinary cumulants. Some examples are as follows: \begin{eqnarray*} \kappa^{rs, t} &=& \kappa^{r,s,t} + \kappa^r\kappa^{s,t} + \kappa^s \kappa^{r,t}\\ &=& \kappa^{r,s,t} + \kappa^r\kappa^{s,t}[2]\\ \kappa^{rs,t u} &=& \kappa^{r,s,t,u} + \kappa^{r,s,t}\kappa^u[4] + \kappa^{r,t}\kappa^{s,u}[2] + \kappa^{r,t}\kappa^s\kappa^u[4]\\ \kappa^{rs,t,u} &=& \kappa^{r,s,t,u} + \kappa^{r,t,u}\kappa^s[2] + \kappa^{r,t}\kappa^{s,u}[2] \end{eqnarray*} Each generalized cumulant is associated with a partition \(\tau\) of the given set of indices. For example, \(\kappa^{rs,t,u}\) is associated with the partition \(\tau=rs|t|u\) of four indices into three blocks. Each term on the right is a cumulant product associated with a partition \(\sigma\) of the same indices. The coefficient is one if the least upper bound \(\sigma\vee\tau\) has a single block, otherwise zero. Thus, with \(\tau=rs|t|u\), the product \(\kappa^{r,s}\kappa^{t,u}\) does not appear on the right because \(\sigma\vee\tau = rs|t u\) has two blocks.

As an example of the way these formulae may be used, let \(X\) be a scalar random variable with cumulants \(\kappa_1,\kappa_2,\kappa_3,\ldots\). By translating the second formula in the preceding list, we find that the variance of the squared variable is \( {\rm var}(X^2) = \kappa_4 + 4\kappa_3\kappa_1 + 2\kappa_2^2 + 4\kappa_2\kappa_1^2, \) reducing to \(\kappa_4 + 2\kappa_2^2\) if the mean is zero.

Exponential families

Let \( f\) be a probability distribution on an arbitrary measurable space \( ({\mathcal X},\nu)\) , and let \( t\colon{\mathcal X}\to{\cal R}\) be a real-valued random variable with cumulant generating function \( K(\cdot)\) , finite in a set \( \Theta\) containing zero in the interior. The family of distributions on \( {\mathcal X}\) with density \( f_\theta(x) = e^{\theta t(x)} f(x) / M(\theta) = e^{\theta t(x) - K(\theta)} f(x) \) indexed by \( \theta\in\Theta\) is called the exponential family associated with \( f\) and the canonical statistic~\( t\) . In statistical physics, the normalizing constant \( M(\theta)\) is called the partition function.

Two examples suffice to illustrate the idea. In the first example, \( {\mathcal X} = \{1,2,\ldots\}\) is the set of natural numbers, \( f(x) \propto 1/x^2\) and \( t(x) = -\log(x)\) . The associated exponential family is \( f_\theta(x) = x^{-\theta}/\zeta(\theta)\) , where \( \zeta(\theta)\) is the Riemann zeta function with real argument \( \theta > 1\) .

In the second example, \( {\mathcal X}={\cal X}_n\) is the symmetric group or the set of permutations of \( n\) letters, \( x\in{\mathcal X}_n\) is a permutation, \( t(x)\) is the number of cycles, \( f(x) = 1/n!\) is the uniform distribution, and \( M_n(\xi) = \Gamma(n+e^\xi)/(n!\, \Gamma(e^\xi))\) for all real~\( \xi\) . The exponential family of distributions on permutations of \( [n]\) is \( f_{n,\theta}(x) = \frac{\Gamma(\lambda)\, \lambda^{t(x)}} {\Gamma(n+\lambda)}, \) the same as the the distribution generated by the Chinese restaurant process with parameter \( \lambda = e^\theta\) . The associated marginal distribution on partitions, the Ewens distribution on partitions of \( [n]\) , is also of the exponential-family form with canonical statistic equal to the number of blocks or cycles. This number \( t(x)\) is a random variable whose cumulants are the derivatives of \( \log M(\cdot)\) evaluated at the parameter~\( \theta\) .

In the multi-parameter case, \( t\colon{\mathcal X}\to{\cal R}^p\) is a random vector and \( \xi\colon{\mathcal R}^p\to{\cal R}\) is a linear functional, \( M(\xi) = E(e^{\xi(t)})\) is the joint moment generating function. It is sometimes convenient to employ Einstein's implicit summation convention in the form \( \theta(t) = \theta_i t^i\) where \( t^1,\ldots, t^p\) are the components of \( t(x)\) , and \( \theta_1,\ldots, \theta_p\) are the coefficients of the linear functional. For simplicity of notation in what follows, \( {\mathcal X}={\cal R}^p\) and \( t(x) = x\) is the identity function. An exponential-family distribution in \( {\mathcal R}^p\) has the form \( f_\theta(x)=\exp(x^j\theta_j-g(x)-\varphi(\theta)) \) for given functions \( g\) and \( \varphi\) . Integration shows that the distribution \( f_\theta\) has cumulant generating function \( K_\theta(\xi)=\varphi(\theta+\xi)-\varphi(\theta)\) . The cumulants of \( X\sim f_\theta\) are equal to the derivatives of \( \varphi\) at the parameter~\( \theta\) .

Calculus of cumulants

The umbral calculus is a syntax or formal system consisting of certain operations on objects called umbrae, mimicking addition and multiplication of independent real-valued random variables \citep{1994}. To each real-valued sequence \( 1, a_1, a_2,\ldots\) there corresponds an umbra \( \alpha\) such that \( E(\alpha^r) = a_r\) . This freedom gives rise to special umbrae, the singleton and Bell umbra, corresponding to no real-valued random variable. Using these special umbrae, one develops the formal notion of an \( \alpha\)-cumulant umbra \( \chi\cdot\alpha\) by formal product operations in the syntax. Properties of cumulants, \( k\) -statistics and other polynomial functions are then derived by purely formal combinatorial operations. \citet{2008} present details.

Approximation of distributions

Edgeworth approximation

Suppose that \(Y\) is a random variable that arises as the sum of \(n\) independent and identically-distributed summands, each of which has mean \(0\), unit variance, and cumulants \(\kappa_r\), and \(X=Y/\sqrt{n}\). For ease of exposition, assume that cumulants of all orders exist. Then, using (\ref{ndep}), the cumulant generating function of \(X\) is given by \(K(\xi)=\xi^2/2 +\kappa_3\xi^3/(6\sqrt{n}) +\kappa_4\xi^4/(24 n) +\cdots\), and the moment generating function of \(X\) is given by \( K(\xi)=\exp(\xi^2/2)\exp(\kappa_3\xi^3/(6\sqrt{n})+\kappa_4\xi^4/(24 n)+\cdots) \) Exponentiating the second factor gives \( K(\xi)=\exp(\xi^2/2)\left(1\!+\!{{\kappa_3\xi^3}\over{6\sqrt{n}}}\!+\!{{\kappa_4\xi^4}\over{24 n}}\!+\!\cdots\!+\! {\textstyle{\frac12}} \left[ {{\kappa_3\xi^3}\over{6\sqrt{n}}}\!+\!{{\kappa_4\xi^4}\over{24 n}}\!+\!\cdots\right]^2\!+\!\!\cdots\right). \) Reordering terms in powers of sample size, <math kseries> =\exp(\xi^2/2)\left(1+{{\kappa_3\xi^3}\over{6\sqrt{n}}}+{{\kappa_4\xi^4}\over{24 n}}+ {{\kappa_3^2\xi^6}\over{72 n}}+\cdots\right). </math> Repeated application of integration by parts to (\ref{mgfdef}) shows that <math mgfderiv> \xi^r M(\xi) =\int_{-\infty}^\infty\exp(\xi x)(-1)^r f^{(r)}(x)~d x, </math> where \(f^{(r)}\) denotes the derivative of \(f\) of order \(r\). Relation (\ref{mgfderiv}) holds if \(f\) and its derivatives go to zero quickly as \(\vert x\vert\to\infty\). Applying (\ref{mgfderiv}) to the normal density \(\phi(x)=\exp(-x^2/2)/\sqrt{2\pi}\), and applying the result to (\ref{kseries}), gives \( M(\xi)\approx\int_{-\infty}^\infty\exp(\xi x)\phi(x)\left[1+{{\kappa_3 h^3(x)}\over{6\sqrt{n}}}+{{\kappa_4h^4(x)}\over{24 n}}+ {{\kappa_3^2h^6(x)}\over{72 n}}\right]~d x \) for \(h^r(x)=(-1)^r\phi^{(r)}(x)/\phi(x)\), and, since the relationship giving the moment generating function in terms of the density is invertible, and that the inversion process is properly smooth, \citet{1907} approximates the density of \(X\) by <math edser> e_4(x)=\phi(x)\left[1+{{\kappa_3 h^3(x)}\over{6\sqrt{n}}}+{{\kappa_4h^4(x)}\over{24 n}}+ {{\kappa_3^2h^6(x)}\over{72 n}}\right]. </math> In fact, when the summands contributing to \(S\) have a density and cumulants of order at least 5, the error in the approximation, multiplied by \(n^{3/2}\), remains bounded. The functions \(h^r\) defined above are the Hermite polynomials. The approximation (\ref{edser}) is known as the Edgeworth series. The subscript refers to the number of cumulants used in its definition. This series can be used to approximate either the cumulative distribution function or survival function through term-wise integration.

The preceding discussion is intended to be heuristic; \citet{2006} presents a rigorous derivation, along with the natural extension to random vectors.

Saddlepoint approximation

The approximation (\ref{edser}) to the density \(f(x)\) has the property that \(|f(x)-e_r(x)|\leq C n^{-(r-1)/2}\), for some constant \(C\), when the cumulant of order \(r+1\) exists; \(C\) does not depend on \(x\). A similar bound holds for the relative error \((f(x)-e_r(x))/f(x)\), only when \(x\) is restricted to a finite interval. Because of the polynomial factor multiplying the first omitted term in (\ref{edser}), the relative error can be expected to behave poorly. One might prefer an approximation that maintains good behavior for values of \(X\) in a range that increases as \(n\) increases; specifically, one might prefer an approximation that performs well for values of \(\bar Y=X/\sqrt{n}\) in a fixed interval.

Assume again that random variables \(Y_j\) are independent and identically distributed, each with a cumulant generating function \(K(\xi)\) finite for \(\xi\) in a neighborhood of \(0\). As above, define the exponential family \( f_{\bar Y}(\bar y;\theta)=\exp(\theta\bar y-K(\theta))f_{\bar Y}(\bar y). \) One can then choose a value of \(\theta\) depending on \(\bar y\) that makes \(f_{\bar Y}(\bar y;\theta)\) easy to approximate, and the exponential family relationship to derive an approximation for \(f_{\bar Y}(\bar y)\). Conventionally we choose \(\hat\theta\) to satisfy <math speqn> K'(\hat\theta)=\bar y; </math> this makes the expectation of the distribution with density \(f_{\bar Y}(\cdot;\hat\theta)\) equal to the observed value. One then applies (\ref{edser}), with the scale of the ordinate changed to reflect the fact that we are approximating the distribution of \(X/\sqrt{n}\), to obtain \( f_{\bar Y}(\bar y)\approx\exp(-\hat\theta\bar y+K(\hat\theta)) n\phi(0)\left[1+{{\kappa_3 h^3(0)}\over{6\sqrt{n}}}+{{\kappa_4h^4(0)}\over{24 n}}+ {{\kappa_3^2h^6(0)}\over{72 n}}\right]. \) Using the fact that \(h^3(0)=0\), \(h^4(0)=3\), and \(h^6(0)=-15\), we obtain <math spser> f_{\bar Y}(\bar y)\approx{{n}\over{\sqrt{2\pi}}} \exp(K(\hat\theta)-\hat\theta\bar y) \left[1+{{\hat\kappa_4}\over{8 n}}- {{5\hat\kappa_3^2}\over{24 n}}\right]. </math> Here \(\hat\kappa_j\) are calculated from the derivatives of \(K\) in the preceding manner, but in this case evaluated at \(\hat\theta\). This approximation may only be applied to values of \(\bar y\) for which (\ref{speqn}) has solutions in an open neighborhood of 0. Expression (\ref{spser}) represents the saddlepoint approximation to the density of the mean \(\bar Y\); since \(f_{\bar Y}(\bar y;\theta)\) has a cumulant generating function defined on an open set containing \(0\), cumulants of all orders exist, the Edgeworth series including \(\kappa_6\) may be applied to \(f_{\bar Y}(\bar y;\theta)\), and so the error in the Edgeworth series is of order \(O(1/n^2)\). Hence the error in (\ref{spser}) is of the same order, and in this case, is relative and uniform for values of \(\bar y\) in a bounded subset of an open subset on which (\ref{speqn}) has a solution. This approximation was introduced to the statistics literature by \citet{1954}.

The Edgeworth series for the density was trivially integrated to obtain an approximation to tail probabilities. Integration of the saddlepoint approximation is more delicate. Two main approaches have been investigated. \citet{1987} expresses \(f_{\bar Y}(\bar y)\) exactly as a complex integral involving \(K(\xi)\), integrates with respect to \(\bar y\) to obtain another complex integral, and reviews techniques for approximating the resulting integrals. \citet{1982} and \citet{1980} derive tail probability approximations based on approximately integrating (\ref{spser}) with respect to \(\bar y\) directly.

These saddlepoint and Edgeworth approximations have multivariate and conditional extensions. \citet{1988} exploits the conditional saddlepoint tail probability approximation to perform inference in canonical exponential families.

Samples and sub-samples

A function \(f\colon{\mathcal R}^n\to{\cal R}\) is symmetric if \(f(x_1 ,\ldots, x_n) = f(x_{\pi(1)} ,\ldots, x_{\pi(n)})\) for each permutation \(\pi\) of the arguments. For example, the total \(T_n = x_1 + \cdots + x_n\), the average \(T_n/n\), the min, max and median are symmetric functions, as are the sum of squares \(S_n = \sum x_i^2\), the sample variance \(s_n^2 = (S_n - T_n^2/n)/(n-1)\) and the mean absolute deviation \(\sum |x_i - x_j|/(n(n-1))\).

A vector \(x\) in \({\mathcal R}^n\) is an ordered list of \(n\) real numbers \((x_1 ,\ldots, x_n)\) or a function \(x\colon[n]\to{\mathcal R}\) where \([n]=\{1 ,\ldots, n\}\). For \(m \le n\), a 1--1 function \(\varphi\colon[m]\to[n]\) is a sample of size~ \(m\), the sampled values being \(x\varphi = (x_{\varphi(1)} ,\ldots, x_{\varphi(m)})\). All told, there are \(n(n-1)\cdots(n-m+1)\) distinct samples of size~ \(m\) that can be taken from a list of length~ \(n\). A \emph{sequence} of functions \(f_n\colon{\mathcal R}^n\to{\cal R}\) is consistent under sub-sampling if, for each \(f_m, f_n\), \( f_n(x) = {\rm ave} _\varphi f_m(x\varphi), \) where \({\rm ave} _\varphi\) denotes the average over samples of size~ \(m\). For \(m=n\), this condition implies only that \(f_n\) is a symmetric function.

Although the total and the median are both symmetric functions, neither is consistent under sub-sampling. For example, the median of the numbers \((0,1,3)\) is one, but the average of the medians of samples of size two is 4/3. However, the average \(\bar x_n = T_n/n\) is sampling consistent. Likewise the sample variance \(s_n^2 = \sum(x_i - \bar x)^2/(n-1)\) with divisor \(n-1\) is sampling consistent, but the mean squared deviation \(\sum(x_i - \bar x_n)^2/n\) with divisor \(n\) is not. Other sampling consistent functions include Fisher's \(k\)-statistics, the first few of which are \(k_{1,n} = \bar x_n\), \(k_{2,n} = s_n^2\) for \(n\ge 2\), \begin{eqnarray*} k_{3,n} &=& n\sum(x_i - \bar x_n)^3/((n-1)(n-2))\\ k_{4,n} &=& \end{eqnarray*} defined for \(n\ge 3\) and \(n\ge 4\) respectively.

For a sequence of independent and identically distributed random variables, the \(k\)-statistic of order~ \(r\le n\) is the unique symmetric function such that \(E(k_{r,n}) = \kappa_r\). \citet{1929} derived the variances and covariances. The connection with finite-population sub-sampling was developed by \citet{1950}.

\bibliography{cumulants}