Asymptotic Expansion of Central Limit Theorem

1  Introduction

1.1  CLT problem

1.2  Strategy & Overview

• Fourier transform of probability distribution, followed by brief statement of the CLT and the questions it poses.
• A study of moments: normalization, propagation/cumulation, semi-invariants.
• Concise derivation of CLT and limiting values of cumulative moments
• Asymptotic expansion: motivation, Tchebychev-Hermite series, domain of convergence, distribution function.

  • Examples: uniform and photon distributions.

• How large does n have to be for a low-order truncated expansion or the gaussian limit to apply? Answer from the Fourier transform.

  • Examples: uniform and photon distributions.

2  Heuristic introduction to CLT

2.1  Characteristic function

If we expand the exp(i s t) term as a Taylor series, then it is clear that the characteristic function f(s) contains all the moments, [`(t^k)], of F(t).

2.1.1  gaussian example

2.2  Brief statement of CLT

1. How large does n have to be for the CLT to apply?
2. How accurately do you wish to represent the distribution of the cumulative sum?

3. What moments of the distribution are important, and how faithfully should they be represented?

The second question is an indispensable complement to the first and means ``what relative fractional accuracy maintained over what range of the argument?'' Typically we shall consider a range of ±3 standard deviations because 99% of events lie within this range for the gaussian distribution. The third question is particularly profitable when only certain moments contribute to the evolution of some characteristic of a physical system. We shall delay a deeper discussion of the CLT until after introduction of further useful mathematical tools.

3  Moments

3.1  Normalized moments

Let s² º [`(t²)] be the variance, then s is the r.m.s. width about the mean; and for the gaussian distribution s is equal to one standard deviation. If in place of F_t(t) we deal with F_x(x) = sF_t(t = xs) where x º t/s, then we shall have a set of moments which are length-scale invariant. Thus we define the normalized moments as:

3.2  Cumulative moments

Now we have the relation

3.3  Exponential form and semi-invariants

The coefficients b_k are called semi-invariants[2]. The obvious method to find the general coefficient b_k is to form the k^th derivative of A and of B, then set s = 0 and compare the two expressions. This procedure has to be followed with k ascending from 0 to l and back-substitution of all previous solutions b_j(a0,a1,a2,...) with 0 £ j < k to find closed form expressions. The first few terms are as follows:

3.3.1  cumulative moments

3.3.2  remainder form

3.4  Reduced variable

4  Central Limit Theorem

4.0.1  qualitative argument

4.1  Gaussian distributed moments

Hence a possible statement of the CLT is that the probability density function F[T_n/(sÖn)] tends in the limit n®¥ to a distribution whose moments are identical with those of a gaussian Equations (7) or (31).

Another way to use these ideas is to estimate how large must n be that the gaussian distribution gives some moment to some desired absolute accuracy. For example, for the gaussian (which has zero skewness) to give the skewness [`(X<sub>n</sub><sup>3</sup>)] to an accuracy of 1% requires Ön ³ 100×|b<sub>3</sub>|3! = 100[`(x³)]. Likewise, for the gaussian approximation of the true function F(X_n) to give the curtosis [`(X<sup>4</sup><sub>n</sub>)] to a relative fractional accuracy of 1% requires n ³ 100 |b<sub>4</sub>|2!/b<sub>2</sub><sup>2</sup> = 100|[`(x⁴)]/3-1|.

4.2  Gaussian distribution

5  Asymptotic expansions

5.1  Importance of higher order semi-invariants

5.1.1  convolution

When G (or an approximation to it) is known analytically, one can quite profitably use the convolution form (44); and this is used in sections 7.1.2 and 7.2.2. A variant of this idea, which can be numerically more efficient, is to represent H by its truncated Taylor expansion:

5.2  The Tchebychev-Hermite expansion

5.3  Truncation of H and domain of convergence

5.3.1  sufficiently large n

5.4  Symmetric distribution

5.4.1  sufficiently large n

The conditions (53) and (57) are much more easily derived (section 7) in the transform space s by direct consideration of the characteristic function.

5.5  Convolution

For a given domain [^z], there is a minimum n for which the truncated Taylor series converges, and its value depends on the form of F_t, the single event probability density. For example, for the top hat function n ³ 5 whereas for the single photon spectrum n ³ 270 when [^z] = 3.

5.6  Distribution function

Using the computer program Mathematica we have created an expansion of E up to order 1/Ö{n¹³} and up to order z³⁸ with the 12 moments as free variables. Once n and the moments are substituted, the expansion collapses to a Taylor series in powers of z. Fortran-77 subroutines have been written to evaluate the expansion and generate coefficients for the Taylor series. In addition, we have written a further subroutine to evaluate the inverse D^-1. Typically the correction term E(z)e^-z2/2/Ö{2p} is much smaller than F(z). Hence we can adopt a two step strategy. In the first step we invert a good approximation to Erf. In the second step we use Newton-Raphson iteration to solve D_t(z) = q for z using an exact form for Erf. A useful and highly accurate approximation is

5.7  Acknowledgement

6  Examples

6.1  Uniform distribution

6.1.1  moments & semi-invariants

6.1.2  distribution of cumulative sum

The Figure 3 below shows shows the result of inverting D_Zn for n = 5 and 10⁵ random trials.

Figure 1: H(z,5) constructed by asymptotic expansion versus Z_n

Figure 2: Exact H(z,5) by numerical integration (of FT^-1) versus Z_n

In section 7.1.2, an approximation to H based on retaining the curtosis is given which is claimed to be valid for n ³ 20. We invite the reader to study Figures 4, 5 and 14; the agreement between them is excellent.

When the asymptotic expansion of H(z,n) is truncated at 1/Ö{n¹³} and n = 20, the domain of convergence is [^z] » 4.3. Consequently, we can graph and compare the expansion with a numerical evaluation out to Z_n = 4, as below Figures 6 and 7.

Figure 4: H(z,20) by asymptotic expansion versus Z_n

Figure 5: Exact H(z,20) by numerical integration versus Z_n

Figure 6: H(z,20) constructed by asymptotic expansion versus Z_n

Figure 7: Exact H(z,20) by numerical integration (of FT^-1) versus Z_n

6.2  Single photon distribution

6.2.1  moments & semi-invariants

Now we expand the characteristic function as a Taylor series: FT{P_x} = A(s). We find:

6.2.2  distribution of cumulative sum

Figure 9: H(z,270) versus Z_n

Figure 10: H(z,6000) versus Z_n

The Figure 11 below shows shows the result of inverting D_Zn for n = 270 and 10⁵ random trials; the departure from a gaussian distribution is obvious.

In section 7.2.2, an approximation to H based on retaining the skewness is given which is claimed to be valid for n ³ 6000. We invite the reader to compare Figures 10 and 20; the agreement between them is excellent.

7  How large `n' for the CLT to apply

Because the range of the gaussian exp(-s²/2) is intimately related to its transform partner exp(-z²/2) so it follows that, for n large enough to satisfy the above conditions, the true F_Zn will be well approximated (i.e. to 1% or better) by a gaussian over the Z_n-range ± three standard deviations, or the T_n-range of ±3Öns. If an approximation of F_Zn is required that is accurate over a smaller range of Z_n, then one can relax the 3 standard deviations requirement appearing in equations (71); and the ``sufficiently large value of n" will decrease.

Let us denote by N_d the n for which a single term dominates nR(s/Ön) for |s| < 3. Further, let us call the value of n for which this term becomes negligible N_¥. Once a desired accuracy for F_Xn and the range over which that accuracy must be maintained have been chosen, then the sufficiently large n for the CLT to be applicable, N_¥, can be computed solely from a knowledge of the semi-invariants b_k. The above arguments are not a proof. However, condition (71) is identical with that obtained in the Z_n-space and, further, the two detailed numerical examples we shall present confirm all of these expectations. These examples are based on the notion of analytically performing the inverse Fourier transform of the remainder e^nR, but truncating R(s) at either b₃s³ or b₄s⁴. If we retain b₃s³ then we reliably model the evolution of skewness; and if we keep b₄s⁴, then we accurately model the curtosis. In fact, because of the exponentiation, the terms b₃ and b₄ ``feed up" and so also contribute to improving the estimates of the 5<sup>th</sup>, 6<sup>th</sup> and 7<sup>th</sup> moments, etc, even if we drop the terms b<sub>5</sub>,b<sub>6</sub>,b<sub>7</sub> from r(s). For the range of values N<sub>d</sub> < n < N<sub>¥</sub>, these numerical examples may be used as a ``bench mark'' for the asymptotic expansions derived earlier. In the following sections we apply these ideas to our two ``working examples''.

7.1  Uniform distribution

7.1.1  N_d and N_¥

We may draw the following conclusions, valid for |s| < 3.

• For n < 26, the term b₄s⁴/n is not an adequate approximation for the remainder function n×R(s/sqrtn).
• The term in s⁴ becomes the dominant remainder term for n ³ N_d = 26.
• The dominant remainder term becomes negligible for n ³ N_¥ = 405.
• For n > 405, the approximation of F_Zn by a gaussian distribution will be accurate to 1% or better over the range Z_n = [-3,+3].

7.1.2  Retain curtosis in characteristic function

From the convolution integral (74) we may write

Figure 14 shows H(z,20) plotted versus Z_n. For n = 20, the b₄s⁴ terms starts to become a good approximation to R(s). However, as can be seen in figure 15 using the gaussian term alone to represent F_Zn would yield an error of up to 7% at Z_n = 3.

Figure 14: H(z,20) versus Z_n.

Figure 15: H(z,20) versus Z_n.

Figure 16 shows H(z,400) plotted versus Z_n. We predicted above that H would differ from unity by only 1%. The numerical values appearing in the plot confirm accuracy to 0.4% out to Z_n = 3. Hence our prediction was slightly conservative. When n = 400, figure 17 indicates that if the gaussian approximation is used in place of the true F_Zn, then the error will be less than 2% out to Z_n = 4.

Figure 16: H(z,400) versus Z_n.

Figure 17: H(z,400) versus Z_n.

7.2  Single photon distribution

7.2.1  N_d and N_¥

We may draw the following conclusions, valid for |s| < 3.

• For n < 6783, the term b₃s³/Ön is not an adequate approximation for the remainder function nR(s/Ön).
• The term in s³ becomes the dominant remainder term for n ³ N_d = 6783.
• The dominant remainder term becomes negligible for n ³ N_¥ = 2.7156×10⁶.
• For n > 2.716×10⁶, the approximation of F_Zn(Z_n) by a gaussian distribution will be accurate to 1% or better over the range Z_n = [-3,+3].

7.2.2  retain the skewness in characteristic function

Figure 18: G(z,1)

As before, from we may write the convolution integral in the form (76) but with G given by Equation (80). Again the task is to show that H(z,n) is approximately constant over the range Z_n = [-3,+3]. Figure 19 shows that H(z,1) attempts to restore the shape of the true single-photon spectrum. However, because for n < 6000 the term b₃s³/Ön is a poor approximation to nR(s/Ön), so it follows that e^-z2/2H(z,1)/Ö{2p} is a poor approximation to the true P_Zn(z).

Figure 20 shows H for large enough n that b₃s³/Ön has become the dominant term and is a sufficient approximation to n×R(s/Ön) for frequencies |s| £ 3. However, n is not yet large enough that F_Zn can be well approximated by the gaussian term alone. H(z,6000) differs from unity by about 9% at the extremities of the range Z_n = [-3,+3].

For very large values of n, ``Mathematica" has problems with numerical accuracy when performing the numerical integrations that give H(z,n). These difficulties arise not because the integrand is oscillatory, but rather because one term of the integrand is essentially 0/0 and taking the limit numerically is prone to errors that can only be avoided by truncating the integral. The problems are anticipated to be much more acute for negative values of z. Figure 21 shows H for n = 10⁵. H(z) is given reliably for z > 0 and indicates that H deviates from unity by about 5% at Z_n = 3.

Figure 20: H(z,6000) versus Z_n.

Figure 21: Numerical approximation to H(z,10⁵) versus Z_n.

Figure 22 shows H for n = 10⁶ large enough that the CLT starts to apply with some confidence. Unfortunately, n is also large enough that the numerical integrals which approximate to H are wildly inaccurate for z < 0. However, the numerical approximation to H is still acceptable for z > 0. Figure 23 shows that H differs from unity by less than 1% at Z_n = 3. The reader will recall that we predicted above (section 7.2.1) that this level of accuracy would not occur until n = 2.7×10⁶, and so our prediction was somewhat conservative.

Figure 22: Numerical approximation to H(z,10⁶).

Figure 23: H(z,10⁶) versus Z_n mostly positive.

8  Conclusion

References

[1]

A.I. Khinchin: Mathematical Foundations of Statistical Mechanics, Dover Publications, New York (1949).

[2]

B.V. Gnedenko & A.N. Kolmorogrov: Limit Distributions for Sums of Independent Random Variables, Addison-Wesley Publishing Co.

[3]

Cramer: Mathematical Methods of Statistics, Princeton Univ. Press (1946).

[4]

M. Loeve: Probability Theory, Vol.I, Springer-Verlag (1977).

[5]

J.M. Jowett: Introductory Statistical Mechanics for Electron Storage Rings, SLAC-PUB-4033 (1986).

[6]

I. Gradshteyn & I. Rhyzik: Integrals, Series and Products, Academic Press (1980).

Footnotes:

¹ This must be so, else the transform f(s) does not exist.