5 ESTIMATING THE VARIATIONS OF THE SAMPLE VARIANCE

In many situations, estimating the variance of the population is not only a step for building confidence intervals around the mean, but is significant by itself. For example, an industrial product cannot be used for task requiring some precision when the adequate test shows a large variance. On the other hand, enforcing a useless precision will only induce additional costs. In production, a sudden increase in variability may indicate the appearance of a production fault citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitecrow:stats60.

In such situations, confidence intervals around the variance are to be discussed. When the moments of the population are not known a priori, formulas like (4) cannot be used, and must be replaced by formulas using the moments of the sample. Therefore, a method is needed to compute the expectation of these moments and of their products.

1 An Algorithm for the Expectation of Products

Definition 5.1   The $m$-degree of a monomial $\beta\doteq\prod m_{k}^{\beta_{k}}$ is ${\mathrm{dg}}_{m}\,\beta\doteq\sum\beta_{k}$ i.e. the number of factors $m_{k}$, distinct or not, occurring in $\beta$, while the $x$- degree of the same monomial is ${\mathrm{dg}}_{x}\,\beta\doteq\sum k\,\beta_{k}$ i.e. the number of factors $x_{i}$, distinct or not, occurring in the expansion of $\beta$ into a polynomial in the $x_{i}$. Accordingly, ${\mathrm{dg}}_{\mu }\,\alpha\doteq\sum\alpha_{j}$ and ${\mathrm{dg}}_{x}\,\alpha\doteq\sum j\,\alpha_{j}$ are defined for a monomial $\alpha\doteq\prod\mu _{j}^{\alpha_{j}}$ depending on the population moments.

In order to obtain an unbiased statistic for a monomial $\alpha$, we cannot simply substitute each $\mu _{j}$ by its unbiased statistic. This is obvious for terms like $\mu _{j}^{2}$, but in fact this ever occurs since the moments of the sample are not (fully) independent when considered as random variates. We have to consider all the monomials $\beta$ such that ${\mathrm{dg}}_{x}\,\beta={\mathrm{dg}}_{x}\,\alpha$. After having computed each $E_{\Phi }\left(\beta\right)$, we can eliminate the irrelevant monomials in the $\mu _{j}$ in order to specifically isolate $\alpha$.

Many research papers have been devoted to the computation of the $E_{\Phi }\left(\beta\right)$, citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitetracy65:moments among them. These computations have been completely transformed by the formal computing tools that nowadays are largely accessible. To quote the Knuth's foreword to citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalciteAeqB :

Science is what we understand well enough to explain to a computer. Art is everything else we do. During the past several years an important part of mathematics has been transformed from an Art to a Science. No longer do we need to get a brilliant insight in order to evaluate sums of binomial coefficients, and many similar formulas that arise frequently in practice. We can now follow a mechanical procedure and discover the answers quite systematically.

Algorithm 5.2   In order to obtain the closed form of a given $E_{\Phi }\left(\beta\right)$ :

1) For each $n$ in $\left[2,\, N\right]$ write $\beta$ as an expanded polynomial in the $x_{i}$ ( $1\leq i\leq n$). Then substitute each power $x_{i}^{j}$ ($j>1$) by the corresponding moment $\mu _{j}$ and thereafter any remaining $x_{i}$ by 0 ($\mu =0$ is assumed). For each $n$, the result is a polynomial $P_{n}=\sum_{\alpha}c\left(n,\alpha\right)\times\alpha$ where the summation ranges over all the $\alpha$ such that ${\mathrm{dg}}_{x}\,\alpha={\mathrm{dg}}_{x}\,\beta$ and the $c\left(n,\,\alpha\right)$ are rational numbers.

2) The closed form of each $c\left(n,\,\alpha\right)$ is a quotient of polynomials in $n$, whose degrees cannot exceed ${\mathrm{dg}}_{x}\,\beta$. They can be obtained by the algorithms described in citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalciteAeqB and implemented as gfun in Maple citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitesalvy:gfun94.

3) In fact, each denominator is a divisor of $n^{p}\left(n-1\right)^{q}$ where $p+q+2={\mathrm{dg}}_{x}\,\beta$ and $q+1={\mathrm{dg}}_{m}\,\beta$. Thereafter, it remains only to obtain the closed form of a polynomial in $n$ from a list of values and Algorithm 5.3 can be used.

Proof. The rule concerning the degrees is obvious. The rationality of the $c\left(n,\,\alpha\right)$ comes from the binomial coefficients, and the specific value of the denominators in closed form comes from products of powers of $1/n$ (from the definition of the sample mean $m$) by powers of $1/\left(n-1\right)$ (from the definition of the $m_{k}$). $\qedsymbol$

Algorithm 5.3 (Newton)   In what follows $_{k}\delta_{j}$ denotes the $j$-th element of a list named $_{k}\delta$. Let $_{0}\delta$ contains the values taken by a polynomial $p$ at some prescribed abscissas $n_{0},\, n_{1},\,\,\, n_{N}$. In other words, $_{0}\delta_{j}=p\left(n_{j}\right)$. In order to determine $p$, $N\geq{\mathrm{dg}}_{p}\,$ is assumed. For increasing $k$, compute the divided differences $_{k}\delta$ as defined by :

$$_{k}\delta_{j}=\left(_{k-1}\delta_{j+1}-{}_{k-1}\delta_{j}\right)/\left(x_{j+k}-x_{j}\right)$$

Then $p\left(x\right)=\sum_{k=0}^{N}\,_{k}\delta_{0}\prod_{j=0}^{k-1}\left(x-n_{j}\right)$

Even automated, these computations are prone to errors and typos. For example, in page 208 of citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitefisher:moments29, we should have $\cdots\left.1120*n-1120\right)\mu _{3}^{2}\,\mu_{2}$ in the $\mu\left(3^{2}2\right)$ formula instead of $1120*n+1120$. While using his algorithm to compute cumulants by identifications, citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitegood77:kstats has detected a typo in citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalciteud-din54, while this article itself was signaling a typo from another author. An efficient test of correctness is the following :

Proposition 5.4   For a given degree $d={\mathrm{dg}}_{x}\,\beta$ , the determinant of all the $E_{\Phi }\left(\beta\right)$ over the basis of all the $\alpha$ of the same degree splits into linear factors, namely $n$ and $\left(n-1\right)$ in the denominator and $\left(n-j\right)$ where $2\leq j\leq d-1$ in the numerator. For example, when ${\mathrm{dg}}_{x}\,\beta=4,\,5,\,6,\,10,\,11$, we have :

$$\Delta_{\,4}=\frac{\left(n-2\right)\left(n-3\right)}{n\,\left(n-1... ...ft(n-3\right)^{2}\left(n-4\right)\left(n-5\right)}{n^{3}\,\left(n-1\right)^{4}}$$

$$\Delta_{10}=\frac{\left(n-2\right)^{11}\left(n-3\right)^{10}\left... ...(n-7\right)^{2}\left(n-8\right)\left(n-9\right)}{n^{18}\,\left(n-1\right)^{22}}$$

$$\Delta_{11}=\frac{\left(n-2\right)^{14}\left(n-3\right)^{12}\left... ...t(n-8\right)^{2}\left(n-9\right)\left(n-10\right)}{n^{28}\left(n-1\right)^{27}}$$ (9)

Proof. Clearly, the denominator splits into powers of $n$ and $n-1$. To prove the form of the numerator, let us consider $1/\Delta$. When expressing a given $\alpha$ in the vector space spanned by the $E_{\Phi }\left(\beta\right)$, the l.c.m. of the involved denominators is necessarily $\prod\left(n-j\right)$ since (1) the degree of this polynomial has the maximal value and (2) each $n-j$ is required since the $E_{\Phi }\left(\beta\right)$ cannot form a basis when $j<{\mathrm{dg}}_{x}\,\alpha$. When computing the determinant of the $\alpha$ over the $E_{\Phi }\left(\beta\right)$, this denominator is elevated to power $\char93 \beta$ : any exponent in the numerator cannot exceed this value (in fact, they form a decreasing sequence). $\qedsymbol$

Remark 5.5   It appears that no factors $n-1$ are canceling. When ${\mathrm{dg}}_{x}\,\beta=11,$ this leads to power :

$$27=\sum_{{\mathrm{dg}}_{x}\,\beta=11}\left({\mathrm{dg}}_{m}\,\beta-1\right)$$

At the same time, many powers of $n$ are canceling, reducing the total degree of denominator from $126$ ($14$ denominators, each of ninth degree) to only $27+28=55$.

2 Some Results

As a direct application of this algorithm, we have the following results.

Well Known Result 5.6   When $n>$3, let us define $V$ by (10). Then $E_{\Phi }\left(V\right)=var_{\Phi }\left(m_{2}\right)$ .

$$V\doteq\frac{1}{\left(n-2\right)\left(n-3\right)}\left(\sum\left(x-m\right)^{4}-\frac{n^{2}-3}{n}\,m_{2}^{2}\right)$$ (10)

Proof. The value of $var_{\Phi }\left(m_{2}\right)$ is given by (4). It's total degree is $4$. It can be seen citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitefisher:moments29 that :

$$E_{\Phi }\left(m_{2}^{2}\right)=\frac{1}{n}\,\mu_{4}+\frac{n^{2}-... ...(m_{4}\right)=\frac{n^{2}-3n+3}{n^{2}}\,\mu_{4}+\frac{6n-9}{n^{2}}\,\mu_{2}^{2}$$ (11)

The result follows by elimination. The $\left(n-2\right)\left(n-3\right)$ denominator "comes" from the fact that an undetermined expression is needed when $n=2$ and $n=3$ since, for these values, we have either $E_{\Phi }\left(m_{4}\right)=2\,E_{\Phi }\left(m_{2}^{2}\right)$ or $E_{\Phi }\left(m_{4}\right)=E_{\Phi }\left(m_{2}^{2}\right)$, reducing the dimension of the vector space. $\qedsymbol$

Proposition 5.7   The skewness of statistic $m_{2}$ is given by :

$$\gamma_{1}\left(m_{2}\right) = \frac{1}{\sqrt{n-1}}\,\frac{\left(n-1\right)^{2}\kappa_{6}+12\, n... ...,\left(\left(n-1\right){\it\kappa_{4}}+2\, n\,{\it\kappa_{2}}^{2}\right)^{3/2}}$$

$$= {\displaystyle \frac{1}{\sqrt{n}}\,\left(\frac{{\it\kappa_{6}}+12... ...2\,{\it\kappa_{2}}^{2}\right)^{3/2}}+\mathrm{O}\left(\frac{1}{n}\right)\right)}$$ (12)

Remark   When $\varphi$ is Gaussian, all cumulants are 0 except from $\kappa_{2}$ and this formula reduces to $\sqrt{8/\left(n-1\right)}$, as it should be, since $m_{2}$ is $\chi_{n-1}^{2}$ distributed in this special case. In the general case, the distribution of $m_{2}$ is necessarily skew, and (12) shows that the asymptotic skewness is ever at most $\mathrm{O}\left(1/\sqrt{n}\right)$.

Proposition 5.8   The variance of estimator $V$ is :

$$\begin{multline*} var(n\,V)= \,\,{\frac {{1}}{n}}\,\mu_8 -{\frac {4 \left( n-7 \... ...-4\,n+7 \right) }{ \left( n-1 \right) ^{2}n}}\,{\mu_3}\,{\mu_5} \end{multline*}$$

$$\begin{multline*} var(n\,V)= \,\, {\frac {1}{n}}\,{\kappa_8}+{\frac {24}{n-1}}\,... ...5 \right) }{ \left( n-1 \right) ^{2}}} \,{\kappa_3}\,{\kappa_5} \end{multline*}$$

Remark   According to citeseppen-1000cite##1##2(##1@tempswa , ##2)##1##2##3##1 ##3internalcitefisher:moments29, a better looking expression is obtained by transcoding moments into cumulants. But simplification is only on rational, exactly known coefficients. The underlying complexity, due to such a number of quite canceling terms remains the same : in the cumulants formula, all signs are positive, but the cumulants themselves aren't necessarily positive (even the cumulants of even index).