Exact Distribution of the Sample Variance from a
Gamma Parent Distribution
T. Royen
Fachhochschule Bingen, University of Applied Sciences,
Berlinstrasse 109, D–55411 Bingen, Germany
e-mail: royen@fh-bingen.de
Abstract
Several representations of the exact cdf of the sum of squares of n independent
identically gamma–distributed random variables Xi are given, in particular by a
series of gamma distribution functions. Using a characterization of the gamma
distribution by Laha, an expansion of the exact distribution of the sample variance
is derived by a Taylor series approach with the former distribution as its leading
term. In particular for integer orders α some further series are provided, including
a convex combination of gamma distributions for α = 1 and nearly of this type
for α > 1. Furthermore, some representations of the distribution of the angle Φ
between (X1, ...,Xn) and (1, ...,1) are given by orthogonal series. All these series
are based on the same sequence of easily computed moments of cos(Φ).
AMS 2000 Subject Classifications: 62E15, 62H10
Keywords: Exact distribution of quadratic forms in non–normal random vari-
ables, Exact distribution of the sample variance, Gamma distribution, Exponential
distribution
1 Introduction
The distribution of the sample variance s2 in non–normal cases has attracted sporadic
attention during the last eight decades. Early investigations, concerning a gamma par-
ent distribution can be traced at least to Craig (1929) and Pearson (1929), who used
moment approximations, later investigated more thoroughly by Bowman and Shen-
ton(1983). For general distributions various approximations were proposed by Box
(1953), Roy and Tiku (1962), Tan and Wong (1977) and by Mudholkar and Trivedi
(1981). The latter authors recommend transformations of the Wilson–Hilferty type.
Exact results seem to be available only for a mixture of two normal distributions (Hyre-
nius 1950, Mudholkar and Trivedi 1981). See also chapter 4.7 in Mathai and Provost
(1992). The interest in the distribution of s2 or 1/s was also stimulated by investiga-
tions on the estimation of process capability indices, see e.g. Pearn, Kotz and Johnson
(1992) and Pearn and Kotz (2006).
http://arxiv.org/abs/0704.1415v1
Apart of its intrinsic theoretical value, an analytical representation of the cdf of s2
enables more rapid and more accurate evaluations than Monte Carlo methods. Besides,
it provides a base to investigate the accuracy of the different proposed approximations.
Now let X1, . . . ,Xn be a random sample with mean X̄ and sample variance
S2 = (n− 1)−1 ∑ni=1(Xi − X̄)2 from a gamma distribution with pdf
gα(x) = x
α−1 exp(−x)/Γ(α), (1)
and cdf
Gα(x) =
gα(ξ )dξ .
Xi, Z :=
X2i , U
2 := ZY−2, (2)
we have
(n− 1)S2 = (U2 − 1/n)Y2, 1/
n ≤U ≤ 1. (3)
By Laha (1954) — see also Lukacs (1955) — the mutual independence of X and the
sample coefficient of variation S/X =
n((U2 − 1/n)/(1− 1/n))1/2 was shown to
characterize the family of gamma distributions indexed by α . This is the key to derive
the distribution of S2 essentially from the distribution of Z. At first, it follows with the
angle Φ between (X1, . . . ,Xn) and (1, . . . ,1) that
Hα ,n(r) := Pr{Z ≤ r2} = EU(Gαn(r/U))
= EΦ(Gαn(r
ncosΦ)),
n− 1S ≤ r} = EU(Gαn(r(U2 − 1/n)−1/2)
= EΦ(Gαn(r
ncotΦ)).
With the cdf Fα ,n of U it follows also
Hα ,n(r) = EY (Fα ,n(r/Y ))
= rαn(Γ(αn))−1
xαn−1Fα ,n(x
−1)e−rxdx.
h∗α ,n(r) := r
−αnHα ,n(r) (7)
is the Laplace transform (L.t.) of
f ∗α ,n(x) := (Γ(αn))
−1xαn−1Fα ,n(x
−1), 0 < x ≤
n. (8)
A power series for h∗α ,n – and consequently the cdf Hα ,n – is derived in (12) . . . (16)
of the follwing section. Due to (7), (8) the coefficients of this series determine the
moments of f ∗α ,n and therefore Fα ,n and the cdf in (5).
For a cdf a probability mixture representation is more appealing than a power series.
Such a mixture of gamma distributions for Hα ,n is found in (20) in section 2.
Then the reader might go directly to theorem 4.1, which is the main result con-
cerning the actual computation of the cdf of S2. By means of a few lines of code of a
computer algebra system some tables of this cdf with at least eight correct digits have
been computed for some values of α and n.
Representations for f ∗α ,n by an orthogonal series with Legendre polynomials and
by a Fourier sinus series are derived too in section 2. Because of the relations (7), (8)
these two series provide two further representations of h∗α ,n or Hα ,n. The corresponding
series for Fα ,n(u) in (31) and the orthogonal series derived from (32) can be used within
(5). However, the application of theorem 4.1 is simpler and more accurate.
In spite of the numerical use of the double series in theorem 4.1 it would be theo-
retically more satisfying to have a single alternating power series or even a probability
mixture representation for the cdf of S2. The latter is accomplished exactly for the
exponential case (α = 1) and nearly for integer α > 1 by theorem 4.2, which contains
also an alternating power series. However, for numerical purposes, theorem 4.2 is (at
present) not considered as a competitor to theorem 4.1 since the computation of the
required coefficients is more cumbersome. The method is explained in section 3. It
is based on a representation of the cdf of tanΦ by a polynomial on a certain section
[0,ϕn] of the domain of Φ. This formula is found in theorem 3.1. The proofs of all the
theorems are given in the appendix.
Throughout the paper formulas from the handbook of mathematical functions by
Abramowitz and Stegun are cited by A.S. and their number. The symbol ∑(k) stands
for summation over all possible decompositions k = k1+ . . .+kn with non–negative in-
tegers k j. Moments of positive random variables occur also with non-integer exponents
and for defective distributions.
2 The distribution of the sum of squares
The L.t. of Z1 = X
1 is given by
ψα(t) = (2Γ(α))−1
zα/2−1 exp(−
z− tz)dz
= (2Γ(α))−1
(−1)k Γ((α + k)/2)
t−(α+k)/2 (9)
= (2Γ(α))−1t−α/2 exp
(α + 1
with Kummer’s confluent hypergeometric function M. In particular
ψ1(t) =
with the cdf Φ(x) = (1+ erf(x/
2))/2 of the standard normal distribution. For in-
teger values (1−α)/2 or (2−α)/2 Kummer’s M is given by Hermite polynomials
(A.S. 13.6.17/18).
The cdf Hα ,n(r) = Pr{Z ≤ r2} can be obtained by the Fourier inversion formula,
but some further representations are useful. If ψα(t) is written as t−α/2βα(t−1/2), then
(ψα(t))n = t−αn/2
βα ,n,k t−k/2 (11)
βα ,n,k = (−1)k(2Γ(α))−n ∑
Γ((α + ki)/2)
(βα(y−1))nyk
dϕ , y = ρeiϕ . (12)
ρ > 0, −π < ϕ ≤ π , k ∈ N0.
The parameter ρ was inserted here only for numerical considerations. The βα ,n,k are
more quickly computed recursively for n = 1,2,4,8, . . . than by the integrals.
Laplace inversion in (11) implies
Hα ,n(r) = r
βα ,n,k
Γ(1+(αn+ k)/2)
rk. (13)
By the series expansion of Gαn in (4) we obtain with the moments
µα ,n,k := E
n cosΦ)αn+k
= n(αn+k)/2γα ,n,k (14)
Hα ,n(r) =
Γ(αn)
µα ,n,k
αn+ k
(−r)k
and by comparison with (13) the relation
(−1)k
µα ,n,k = 2Γ(αn)
βα ,n,k
Γ((αn+ k)/2)
. (16)
Also the following representations of the distribution functions of ∑ni=1 X
2,U and
tanΦ are based essentially on the sequence (µα ,n,k) or equivalently (γα ,n,,k).
A probability mixture representation for Hα ,n by gamma distribution functions is
obtained as follows: With any scale factor λ let be
µα ,n,k,λ = λ
−(αn+k)µα ,n,k . (17)
Then, multiplying
Hα ,n(r) = λ
(λ r)αn−1
Γ(αn)
µα ,n,k,λ (−λ r)k/k!
by e−λ r ∑∞ℓ=0(λ r)
ℓ/ℓ! we obtain
λ gαn(λ r)
δα ,n,k,λ (λ r)
k/k! (18)
with the alternating kth order differences
δα ,n,k,λ = (−∆)kµα ,n,0,λ =
(−1) j
µα ,n, j,λ (19)
and consequently
Hα ,n(r) =
δα ,n,k,λ
αn+ k− 1
Gαn+k(λ r). (20)
In particular with λ =
n this is a probability mixture since
δα ,n,k,√n = E ((cosΦ)
αn(1− cosΦ)k)> 0 and
αn+ k− 1
δα ,n,k,√n = 1 (21)
by virtue of the binomial series with 0≤ 1−cosΦ≤ 1−1/
n. However, for numerical
approximations different values of λ should be more suitable, e.g.
))1/(αn)
= µ1/(αn)α ,n,0 =
2Γ(αn)
(2Γ(α))n
(Γ(α/2))n
Γ(αn/2)
)1/(αn)
with µα ,n,0,λ = 1.
Besides, two further representations of Hα ,n are obtained from the following or-
thogonal series for f ∗αn from (7), (8). With
xk f ∗αn(x)dx =
Γ(αn)
µα ,n,k
αn+ k
, (23)
implied by (14), (15), and the shifted Legendre polynomials
P∗k (y) =
p∗k, jy
j, 0 ≤ y ≤ 1,
(P∗k (y))
2dy =
2k+ 1
, (24)
(A.S. 22.2.11), it follows
f ∗α ,n(x) =
cα ,n,kP
k (x/
cα ,n,k = (2k+ 1)
nαn/2
Γ(αn)
p∗k j
γα ,n, j
αn+ j
Thus, with (A.S. 11.4.26) and the modified spherical Bessel functions (A.S. 10.2.2) we
Hα ,n(r) = r
αn exp
(−1)kcα ,n,k
Ik+1/2
. (26)
Consequently, the Fourier transform of f ∗α ,n is representable by
h∗α ,n(−it) = exp
ikcα ,n,k jk
with the spherical Bessel functions jk from (A.S. 10.1). For absolutely large real t this
series is numerically more suitable than the power series from (15). Such values are
required as coefficients in the following Fourier sinus expansion for f ∗α ,n.
f ∗α ,n(x) =
bα ,n,m sin
, 0 ≤ x ≤
n , (28)
bα ,n,m =
h∗α ,n
n ∑k≥0
m+k odd
(−1)(m+k−1)/2cα ,n,k jk(mπ/2) .
Because of lemma A.2 the function f ∗αn has a square integrable derivative at least for
min(αn,n)> 2 which entails the absolute uniform convergence of the orthogonal series
in (25) and (28). Inserting (28) into (6) leads to
Hα ,n(r) = r
nr2 +m2π2
bα ,n,m
1− (−1)me−r
. (30)
Finally, with the functions
Gβ (ζ ) :=
Γ(β + k/2)
= (Γ(β ))−1M(1,β ;ζ 2)+ (Γ(β + 1/2))−1ζM(1,β + 1/2;ζ 2), β > 0,
= ζ 2(1−β )eζ
Gβ−1(ζ
2)+Gβ−1/2(ζ
, if β > 1,
the following integral representation can be derived from (11), (12), (13):
Hα ,n(r) =
βα(y−1)
G1+αn/2(ry)
y = ρeiϕ , ρ > 0.
3 The distribution of the angle Φ
Two representations of the cdf Fα ,n of U = (
ncosΦ)−1 follow directly from (25) and
(28).
Fα ,n(u) = Γ(αn)uαn−1
cα ,n,kP
= Γ(αn)uαn−1
bα ,n,m sin
Besides, an orthogonal expansion with Legendre polynomials for the density of U2 is
obtained by the moments
E (U2k) =
Γ(αn)k!
Γ(αn+ 2k)∑
Γ(α + 2ki)
Γ(α)ki!
, (32)
derived from (2).
Direct use of these orthogonal series within (5) is possible to get the cdf of S,
but the double series in theorem 4.1 is numerically more favourable. As mentioned
in the introduction it would be theoretically more satisfying to have a single alternat-
ing series or a probability mixture for this cdf. A single power series is not directly
available by a power series expansion of Gαn in (5) since the moments of cotΦ do
not all exist. However, at least for integer α , this obstacle is avoided by splitting
the domain [0,arctan(
n− 1)] of Φ by ϕn = arctan(1/
n− 1), which is the max-
imal angle ϕ for which the whole cone {Φ ≤ ϕ} is contained completely within
{x1, . . . ,xn|x1, . . . ,xn ≥ 0}. For integer α the theorem 3.1 below provides a polynomial
representation of the cdf of tanΦ restricted to 0≤ tanϕ ≤ (n−1)−1/2. This simple rep-
resentation is used within (5) to integrate over tanϕ ≤ (n−1)−1/2. The integration over
tanϕ > (n− 1)−1/2 needs only truncated moments of cotΦ = (1− cos2 Φ)−1/2 cosΦ,
obtained by the binomial expansion with the truncated moments of cosΦ, given below
in (36).
Theorem 3.1 Let α be a positive integer, then
Wα ,n(t) := Pr{tanΦ ≤ t}
Γ(αn)
(Γ(α))nnαn/2
[(α−1)n/2]
aα ,n,2 j
tn−1+2 j
n− 1+ 2 j
, 0 ≤ t ≤ (n− 1)−1/2,
where the coefficients aα ,n,2 j are given as the unique solutions of the linear equations
[(α−1)n/2]
+β n+m− j
aα ,n,2 j =
(2m)!
Γ(1/2+β +m j)/(2m j)! , α = 2β + 1
(2m+ n)!
nm+n/2
Γ(1/2+β +m j)/(2m j + 1)! , α = 2β
m, j = 0, . . . , [(α − 1)n/2].
With any ρ > 0 the sums ∑(m) within the right hand sides of (34) are also given by
e−imϕ
dϕ , α = 2β + 1,
α + 1
α + 1
e−imϕ
dϕ , α = 2β ,
After Kummer’s transformation (A.S. 13.1.27) the above confluent hypergeometric
functions M can also be expressed by Hermite polynomials due to (A.S. 13.6.17/18).
The truncated moments γ̄α ,n,k = Ē ((cosΦ)αn+k), defined by integration over
ϕ ≥ ϕn, follow from theorem 3.1 by straightforward calculation.
γ̄α ,n,k = γα ,n,k −
Γ(αn)
(Γ(α))nnαn/2
[(α−1)n/2]
aα ,n,2 jB
2 + j,
(α−1)n+k+1
2 − j;
with the incomplete beta function B(a,b;x) =
a−1(1− t)b−1dt. In particular with
α = 1 we have
γ̄1,n,k = γ1,n,k −
(n− 1)!
π (n−1)/2
Γ( n−12 )
. (37)
Two further remarks: The equations (34) can also be solved by explicit matrix in-
version. The matrix
2 +β n+m− j
, m, j = 0, . . . ,N = [(α − 1)n/2], has the
structure Γ(x)Diag
. . . ,(x)m, . . .
(x+m)N− j
with factorials (x)m = x(x+ 1) . . .(x+
m−1). The inverse (pm j(x)) of ((x+m)N− j) contains only polynomials pm j of degree
m in its mth row. These polynomials can be obtained by interpolation from the m+ 1
values
pm j(−(N − k)) =
(k = 0, . . . ,m)
(−1) j+k+1
( j+ k−m)!(N− ( j+ k))!
, m− k ≤ j ≤ N − k
0 , otherwise
The cdf F1,n of U1,n can be applied to a goodness of fit test of the hypothesis H0 : F = F0
against H1 : F 6=F0, where F0 is any specified continuous cdf. Under H0 the test statistic
∑ni=1(F0(Xi:k)−F0(Xi−1:k))2, obtained from an ordered sample Xi:k, i= 1, . . . ,k = n−1,
(X0:k =−∞, Xn:k = ∞), has the same distribution as U21,n. The asymptotic normality of
U21,n can be derived from the asymptotic bivariate normal distribution of
(Yn − n)/
(Zn − 2n)/
with correlation ρ = 2/
5. The convergence to a normal distribution
is slow and might be accelerated by a suitable transformation.
4 The distribution of the sample variance
With the functions
Hα ,n,m(r) := E (Gαn+m(r
ncosΦ))
rαn+m
Γ(αn+m)
µα ,n,m+k
αn+m+ k
(−r)k
m ∈ N0, µα ,n,k = E (
ncosφ)αn+k = n(αn+k)/2γα ,n,k from (14) we obtain
Theorem 4.1 Let X1, . . . ,Xn be independent random variables with den-
sity (Γ(α))−1e−xxα−1, α > 0, then the cdf of the sample variance S2 =
(n− 1)−1 ∑ni=1(Xi − X̄)2 is given by
(n− 1)S2 ≤ z = r2
Γ(αn+ 2 j)
Γ(αn) j!n j
Hα ,n,2 j(
z) (39)
Γ(αn)
(αn+ k)/2+ j− 1
µα ,n,2 j+k
αn+ k
(−r)k
δα ,n, j,k,λ
αn+ k− 1
Gαn+k(λ r)
with the alternating kth order differences
δα ,n, j,k,λ =
(−1)ℓ
(αn+ ℓ)/2+ j− 1
µα ,n,2 j+ℓ
λ αn+2 j+ℓ
and any λ > 0. The series in (39) is absolutely convergent with terms of the order
j−(n+1)/2
, where the O–constant depends only on α and n.
As a simple numerical example we obtain with α = 1 and n = 10 the value
Pr{S ≤ 2}= 0.98530379 . . . with at least eight correct digits.
Finally, the method, explained before theorem 3.1, is applied to obtain the alterna-
tive representations in theorem 4.2 for integer α , including in particular a probability
mixture for the exponential case α = 1. The degree of the polynomials in theorem 3.1
and the rate of convergence of the binomial series in (43) limit the numerical use of
theorem 4.2 to moderate values of αn.
The truncated moments
M̄α ,n,k = Ē
(cotΦ)αn+k
(αn+ k)/2+ j− 1
γ̄α ,n,k+2 j (43)
are used with γ̄α ,n,k from (36), obtained by integration over ϕ ≥ ϕn = arccot(
n− 1).
Theorem 4.2 With the truncated moments from (43) and (36), the pdf wα ,n and the cdf
Wα ,n of tanΦ on 0 ≤ ϕ ≤ ϕn from theorem 3.1, the distribution of S is given by each of
the following three formulas (44), (45), (46):
n− 1S ≤ r
Γ(αn)
(−1)k
γ̄α ,n,k
αn+ k
Pαn,k(r
Gαn(r
n/t)wα ,n(t)dt
with the polynomials
Pαn,k(x) =
[k/2]
(αn+ k)/2
xk−2 j
(k− 2 j)!
Γ(αn)
M̄α ,n,k
αn+ k
+ Wα ,n(1/
n− 1)Gαn
n(n− 1)
n(n−1)
gαn(x)Wα ,n
The infinite series in (45) is also given by
∆̄α ,n,k,λ
αn+ k− 1
Gαn+k(λ r) (46)
with any λ > 0 and the alternating kth order differences
∆̄α ,n,k,λ = (−∆)kM̄α ,n,0,λ =
(−1)ℓ
M̄α ,n,ℓ,λ ,
M̄α ,n,ℓ,λ = Ē
ncotΦ
)αn+ℓ
Remarks: In particular, with α = 1 the only coefficient in w1,n(t) is a1,n,0 =
(n− 1)bn−1 with the volume bn−1 = π (n−1)/2/Γ
of the (n− 1)–unit ball. Then
the last two terms in (45) are reduced to
(n− 1)!
n(n(n− 1))(n−1)/2
bn−1Gn−1(r
n(n− 1)). (47)
The truncated moments γ̄1,n,k were given in (37). Only positive coefficients ∆̄α ,n,k,λ
arise in (46) with λ =
n(n− 1) since 1 > 1− cotΦ√
n−1 ≥ 0. The last term in (45) is
bounded by Wα ,n(1/
n− 1)
1−Gαn(r
n(n− 1))
, which is often neglectible. For
actual computations the choice
Ē ((cotΦ)αn)
1−Wα ,n(1/
n− 1)
)1/(αn)
might be more favourable. With this λ the leading term in (46) becomes
1−Wα ,n
Gαn(λ r).
It would be desirable to find a similar probability mixture representation for the cdf
of S2 for all α > 0 but presumedly, also in case of its existence, the computation of the
required coefficients would not be easy.
Acknowledgement: The author would like to thank S. Kotz for his strong interest
in this work and some hints at the literature within a personal communication.
Appendix
Proof of theorem 3.1. Let α be any positive integer and ω the Lebesgue measure on
the (n− 1)–unit sphere Sn−1. Then
(x1 + . . .+ xk)
(A.1)
= 2Γ((αn+ k)/2)−1k!∑
Γ(α + k j)/2)/k j! , α + k j odd.
After an orthogonal transformation x = Ty with y1 = n
−1/2 ∑ni=1 xi = cosϑ1, followed
by transformation to polar coordinates, we obtain with c j = cosϑ j, s j = sinϑ j the same
integral value as
n+ s1τ j(ϑ2, . . . ,ϑn−1))α−1dω
= nk/2
(α−1)n
α − 1
n(ℓ−(α−1)n)/2
n−1− j
j dϑ j
(α−1)n+k−ℓ
n−2+ℓ
1 dϑ1
= nk/2
[(α−1)n/2]
aα ,n,2 jB
(α − 1)n+ k+ 1
. (A.2)
By comparison with (A.1) it follows
[(α−1)n/2]
(α − 1)n+ k+ 1
aα ,n,2 j
(A.3)
= 2n−k/2k!∑
Γ((α + k j)/2)/k j! .
With α = 2β + 1 only even k j = 2m j occur with sum k = 2m, and with α = 2β only
k j = 2m j + 1 with sum k = 2m+ n, which leads to the linear equations in (34).
If ℓ is odd then the coefficients aα ,n,ℓ vanish since
sn−1− jj dϑ j
is a linear combination of the integrals
smnn−1
M j+n−1− j
j dϑ j , (Mn−1 = 0) ,
which are different from zero only for even m2, . . . ,mn, ∑nj=2 m j = ∑
j=1 ℓ j = ℓ.
If t = tanϕ ≤ (n− 1)−1 we obtain with Sn−1,ϕ = Sn−1 ∩{ϑ1 ≤ ϕ} that
Pr{tanΦ ≤ t} = (Γ(α))−n
{ϑ1≤ϕ}
e−x j xα−1j dx j
= (Γ(α))−n
Sn−1,ϕ
ραn−1 exp(−ρ
n+ s1τ j)α−1dω
Γ(αn)
(Γ(α))nnαn/2
Sn−1,ϕ
c−αn1
n+ s1τ j)α−1dω
Γ(αn)
(Γ(α))nnαn/2
[(α−1)n/2]
aα ,n,2 j
sn−2+2 j1 c
−n−2 j
1 dϑ1 ,
where the integrals are given by (tanϕ)n−1+2 j/(n−1+2 j, which provides the asserted
result in (33). ✷
Comparing the series
(−1)k
(z/y)(αn+k)/2
αn+ k
= Γ(αn)Gαn(
with the corresponding one for Γ(αn+ 2 j)Gαn+2 j(
z/y) we obtain the identity
Γ(αn+ 2 j)
Γ(αn)
Gαn+2 j
(A.4)
by differentiation.
Lemma A.1
2αn/2
, y,z > 0. (A.5)
Proof. Using (A.4) and Cauchy’s integral formula for the derivatives of
Γ(αn+ 2 j)Gαn+2 j(
zαn/2+ j−1 exp(−
the bound in (A.5) follows from
( j− 1)!
Γ(αn+ 2 j)Gαn+2 j(
= lim
ζ αn/2+ j−1 exp(−
ζ )(ζ − z)− jdζ
yαn/2−1 exp(−
y/2)dy = 2αn/2Γ(αn)/π ,
where the way of integration is
≤ ϕ ≤
∣R ≥ y ≥−R
Lemma A.2 Let Φ denote the angle between (X1, . . . ,Xn) and (1, . . . ,1) from (4),
Wα ,n(t) the cdf of tanΦ and bn−1 = π (n−1)/2/Γ( n+12 ) the volume of the (n− 1)–unit
ball, then
Wα ,n(t)∼
Γ(αn)bn−1
(Γ(α))nn(α−1/2)n
tn−1, t ↓ 0, (A.6)
and in particular
W1,n(t) =
(n− 1)!
bn−1t
n−1, 0 ≤ t ≤ (n− 1)−1/2. (A.7)
Furthermore
E ((cosΦ)s)∼
Γ(αn)(2π)(n−1)/2
(Γ(α))nn(α−1/2)n
s−(n−1)/2, s ↑ ∞. (A.8)
Proof. Let x,xg denote the arithmetic and the geometric mean of (x1, . . . ,xn ≥ 0 and
ϕ the angle between (x1, . . . ,xn) and (1, . . . ,1). Setting xi = (1+ εi)x with any fixed x
we have the relation
(xi − x)2 =
ε2i = n tan
and for ϕ ↓ 0
xng = x
(1+ εi)∼ xn exp
∼ xn.
It follows
Wα ,n(t) = Pr{Φ ≤ ϕ}=
gα(xi)dxi
= (Γ(α))−n
n(α−1)
g exp(−nx)dx1 . . .dxn
∼ (Γ(α))−nbn−1
x(α−1)n exp(−nx)(x
n tanϕ)n−1dx
which provides (A.6).
For α = 1 this asymptotic relation can be replaced by an equation if ϕ ≤ ϕn =
arctan(1/
n− 1), which is the largest angle ϕ for which the whole cone {Φ ≤ ϕ} is
contained within {x1, . . . ,xn|x1, . . . ,xn ≥ 0}.
To prove (A.8) we obtain from (A.6)
Pr{| ln(cosΦ)| ≤ ε} = Pr{ln(1+ tan2 Φ)≤ 2ε}
∼ Wαn
∼ cα ,nε(n−1)/2, ε ↓ 0,
with a factor cα ,n determined by (A.6). For s → ∞ the relation (A.8) follows from
the Hardy–Littlewood–Karamata Tauber theorem (cf. e.g. chapter 13 in Feller 1971),
applied to
E ((cosΦ)s) = E (exp(−s| ln(cosΦ)|)) . ✷
Proof of theorem 4.1. With any ϑ ∈ (0,1) EU
(z/(U2 −ϑ/n))1/2
can be
represented by the Taylor series
With U−1 =
ncosΦ and the bounds from Lemma A.1 the terms of this series are
absolutely bounded by
2αn/2
(cosΦ)2 j
ϑ j/ j = O(ϑ j j−(n+1)/2)
due to (A.8) with s = 2 j and on O–constant depending only on α and n.
For ϑ ↑ 1 it follows from (6) and (A.4) that
Pr{(n− 1)S2 ≤ z = r2} = EU
U2 − 1/n
Γ(αn+ 2 j)
Γ(αn) j!n j
Gαn+2 j
Γ(αn+ 2 j)
Γ(αn) j!n j
Hα ,n,2 j(
with the functions Hα ,n,2 j from (38).
With any λ > 0 we obtain
Γ(αn+ 2 j)
Γ(αn) j!n j
Hα ,n,2 j(
Γ(αn)
(αn+ k)/2+ j− 1
µα ,n,2 j+k
λ αn+2 j+k
(−1)k
(λ r)αn−1+k
λ 2/n
After multiplication by exp(−λ r)∑∞ℓ=0(λ r)ℓ/ℓ! the series (41) is obtained by integra-
tion over r. ✷
Proof of theorem 4.2. With the pdf wα ,n and the cdf Wα ,n of tanΦ from theorem 3.1
and the truncated moments M̄α ,n,k = Ē
(cotΦ)αn+k
from (43) it follows with (5) that
n− 1S ≤ r} = E
Gαn(r
ncotΦ)
Gαn(r
n/t)wα ,n(t)dt +
Γ(αn)
M̄α ,n,k
αn+ k
After integration by parts and the substitution t → r
n/x the above integral becomes
n(n− 1)
n(n−1)
gαn(x)Wαn
Inserting the series for M̄α ,n,k from (43) formula (44) follows by rearranging the
resulting absolutely convergent double series according to the truncated moments
γ̄α ,n,k = Ē
(cosΦ)αn+k
from (36).
The series in (46) is obtained by the same way as (21). ✷
References
[1] Abramowitz, M. and Stegun, I. (1968) Handbook of Mathematical Functions.
Dover Publications, Inc., New York
[2] Bowman, K.O. and Shenton, L.R. (1983) The Distribution of the Standard Devi-
ation and Skewness in Random Samples from a Gamma Density — A New Look
at a Craig–Pearson Study. Oak Ridge National Laboratory/ CSD–109
[3] Box, G.E.P. (1953) Nonnormality and Tests on Variances, Biometrika 40,
318–335
[4] Craig, C.C. (1929) Sampling when the Parent Population is of Pearson’s Type III,
Biometrika 21, 287–293
[5] Feller, W. (1971) An Introduction of Probability Theory and Its Applications,
Vol. II, 2nd ed. John Wiley & Sons, New York
[6] Hyrenius, J. (1950) Distribution of Student–Fisher’s t in Samples from Com-
pound Normal Functions, Biometrika 37, 429–442
[7] Laha, R.G. (1954) On a Characterization of the Gamma Distribution, Annals of
Mathematical Statistics 25, 784–787
[8] Lukacs, E. (1955) A Characterization of the Gamma Distribution, Annals of
Mathematical Statistics 26, 319–324
[9] Mathai, A.M., and Provost, S.B. (1992) Quadratic Forms in Random Variables.
Marcel Dekker, Inc., New York
[10] Mudholkar, G.S. and Trivedi, M.C. (1981) A Gaussian Approxiamtion to the
Distribution of the Sample Variance for Nonnormal Populations, Journal of the
Americal Statistical Association 76, 479–485
[11] Pearn, W.L. and Kotz, S. (2006) Encyclopedia and Handbook of Process Capa-
bility Indices: A Comprehensive Exposition of Quality Control Measures. (Ser.
on Quality, Reliability and Engineering Statistics, Vol. 12). World Scientific Pub-
lishing Company, Singapore
[12] Pearn, W.L., Kotz, S. and Johnson, N.L. (1992) Distributional and Inferential
Properties of Process Capability Indices, Journal of Quality Technology, 24,
216–230
[13] Pearson, E.S. (1929) Note on Dr. Craig’s paper, Biometrika 21, 294–302
[14] Roy, J. and Tiku, M.L. (1962) A Laguerre Series Approximation to the Sampling
Distribution of the Variance, Sanky ā, Series A, 24, 181–184
[15] Tan, W.Y. and Wong, S.P. (1977) On the Roy–Tiku Approximation of Sample
Variances from Nonnormal Universes, Journal of the American Statistical Asso-
ciation 72, 875–880
	Introduction
	The distribution of the sum of squares
	The distribution of the angle bold0mu mumu 
	The distribution of the sample variance

  Several representations of the exact cdf of the sum of squares of n
independent gamma-distributed random variables Xi are given, in particular by a
series of gamma distribution functions. Using a characterization of the gamma
distribution by Laha, an expansion of the exact distribution of the sample
variance is derived by a Taylor series approach with the former distribution as
its leading term. In particular for integer orders alpha some further series
are provided, including a convex combination of gamma distributions for alpha =
1 and nearly of this type for alpha > 1. Furthermore, some representations of
the distribution of the angle Phi between (X1,...,Xn) and (1,...,1) are given
by orthogonal series. All these series are based on the same sequence of easily
computed moments of cos(Phi).

Introduction
The distribution of the sample variance s2 in non–normal cases has attracted sporadic
attention during the last eight decades. Early investigations, concerning a gamma par-
ent distribution can be traced at least to Craig (1929) and Pearson (1929), who used
moment approximations, later investigated more thoroughly by Bowman and Shen-
ton(1983). For general distributions various approximations were proposed by Box
(1953), Roy and Tiku (1962), Tan and Wong (1977) and by Mudholkar and Trivedi
(1981). The latter authors recommend transformations of the Wilson–Hilferty type.
Exact results seem to be available only for a mixture of two normal distributions (Hyre-
nius 1950, Mudholkar and Trivedi 1981). See also chapter 4.7 in Mathai and Provost
(1992). The interest in the distribution of s2 or 1/s was also stimulated by investiga-
tions on the estimation of process capability indices, see e.g. Pearn, Kotz and Johnson
(1992) and Pearn and Kotz (2006).
http://arxiv.org/abs/0704.1415v1
Apart of its intrinsic theoretical value, an analytical representation of the cdf of s2
enables more rapid and more accurate evaluations than Monte Carlo methods. Besides,
it provides a base to investigate the accuracy of the different proposed approximations.
Now let X1, . . . ,Xn be a random sample with mean X̄ and sample variance
S2 = (n− 1)−1 ∑ni=1(Xi − X̄)2 from a gamma distribution with pdf
gα(x) = x
α−1 exp(−x)/Γ(α), (1)
and cdf
Gα(x) =
gα(ξ )dξ .
Xi, Z :=
X2i , U
2 := ZY−2, (2)
we have
(n− 1)S2 = (U2 − 1/n)Y2, 1/
n ≤U ≤ 1. (3)
By Laha (1954) — see also Lukacs (1955) — the mutual independence of X and the
sample coefficient of variation S/X =
n((U2 − 1/n)/(1− 1/n))1/2 was shown to
characterize the family of gamma distributions indexed by α . This is the key to derive
the distribution of S2 essentially from the distribution of Z. At first, it follows with the
angle Φ between (X1, . . . ,Xn) and (1, . . . ,1) that
Hα ,n(r) := Pr{Z ≤ r2} = EU(Gαn(r/U))
= EΦ(Gαn(r
ncosΦ)),
n− 1S ≤ r} = EU(Gαn(r(U2 − 1/n)−1/2)
= EΦ(Gαn(r
ncotΦ)).
With the cdf Fα ,n of U it follows also
Hα ,n(r) = EY (Fα ,n(r/Y ))
= rαn(Γ(αn))−1
xαn−1Fα ,n(x
−1)e−rxdx.
h∗α ,n(r) := r
−αnHα ,n(r) (7)
is the Laplace transform (L.t.) of
f ∗α ,n(x) := (Γ(αn))
−1xαn−1Fα ,n(x
−1), 0 < x ≤
n. (8)
A power series for h∗α ,n – and consequently the cdf Hα ,n – is derived in (12) . . . (16)
of the follwing section. Due to (7), (8) the coefficients of this series determine the
moments of f ∗α ,n and therefore Fα ,n and the cdf in (5).
For a cdf a probability mixture representation is more appealing than a power series.
Such a mixture of gamma distributions for Hα ,n is found in (20) in section 2.
Then the reader might go directly to theorem 4.1, which is the main result con-
cerning the actual computation of the cdf of S2. By means of a few lines of code of a
computer algebra system some tables of this cdf with at least eight correct digits have
been computed for some values of α and n.
Representations for f ∗α ,n by an orthogonal series with Legendre polynomials and
by a Fourier sinus series are derived too in section 2. Because of the relations (7), (8)
these two series provide two further representations of h∗α ,n or Hα ,n. The corresponding
series for Fα ,n(u) in (31) and the orthogonal series derived from (32) can be used within
(5). However, the application of theorem 4.1 is simpler and more accurate.
In spite of the numerical use of the double series in theorem 4.1 it would be theo-
retically more satisfying to have a single alternating power series or even a probability
mixture representation for the cdf of S2. The latter is accomplished exactly for the
exponential case (α = 1) and nearly for integer α > 1 by theorem 4.2, which contains
also an alternating power series. However, for numerical purposes, theorem 4.2 is (at
present) not considered as a competitor to theorem 4.1 since the computation of the
required coefficients is more cumbersome. The method is explained in section 3. It
is based on a representation of the cdf of tanΦ by a polynomial on a certain section
[0,ϕn] of the domain of Φ. This formula is found in theorem 3.1. The proofs of all the
theorems are given in the appendix.
Throughout the paper formulas from the handbook of mathematical functions by
Abramowitz and Stegun are cited by A.S. and their number. The symbol ∑(k) stands
for summation over all possible decompositions k = k1+ . . .+kn with non–negative in-
tegers k j. Moments of positive random variables occur also with non-integer exponents
and for defective distributions.
2 The distribution of the sum of squares
The L.t. of Z1 = X
1 is given by
ψα(t) = (2Γ(α))−1
zα/2−1 exp(−
z− tz)dz
= (2Γ(α))−1
(−1)k Γ((α + k)/2)
t−(α+k)/2 (9)
= (2Γ(α))−1t−α/2 exp
(α + 1
with Kummer’s confluent hypergeometric function M. In particular
ψ1(t) =
with the cdf Φ(x) = (1+ erf(x/
2))/2 of the standard normal distribution. For in-
teger values (1−α)/2 or (2−α)/2 Kummer’s M is given by Hermite polynomials
(A.S. 13.6.17/18).
The cdf Hα ,n(r) = Pr{Z ≤ r2} can be obtained by the Fourier inversion formula,
but some further representations are useful. If ψα(t) is written as t−α/2βα(t−1/2), then
(ψα(t))n = t−αn/2
βα ,n,k t−k/2 (11)
βα ,n,k = (−1)k(2Γ(α))−n ∑
Γ((α + ki)/2)
(βα(y−1))nyk
dϕ , y = ρeiϕ . (12)
ρ > 0, −π < ϕ ≤ π , k ∈ N0.
The parameter ρ was inserted here only for numerical considerations. The βα ,n,k are
more quickly computed recursively for n = 1,2,4,8, . . . than by the integrals.
Laplace inversion in (11) implies
Hα ,n(r) = r
βα ,n,k
Γ(1+(αn+ k)/2)
rk. (13)
By the series expansion of Gαn in (4) we obtain with the moments
µα ,n,k := E
n cosΦ)αn+k
= n(αn+k)/2γα ,n,k (14)
Hα ,n(r) =
Γ(αn)
µα ,n,k
αn+ k
(−r)k
and by comparison with (13) the relation
(−1)k
µα ,n,k = 2Γ(αn)
βα ,n,k
Γ((αn+ k)/2)
. (16)
Also the following representations of the distribution functions of ∑ni=1 X
2,U and
tanΦ are based essentially on the sequence (µα ,n,k) or equivalently (γα ,n,,k).
A probability mixture representation for Hα ,n by gamma distribution functions is
obtained as follows: With any scale factor λ let be
µα ,n,k,λ = λ
−(αn+k)µα ,n,k . (17)
Then, multiplying
Hα ,n(r) = λ
(λ r)αn−1
Γ(αn)
µα ,n,k,λ (−λ r)k/k!
by e−λ r ∑∞ℓ=0(λ r)
ℓ/ℓ! we obtain
λ gαn(λ r)
δα ,n,k,λ (λ r)
k/k! (18)
with the alternating kth order differences
δα ,n,k,λ = (−∆)kµα ,n,0,λ =
(−1) j
µα ,n, j,λ (19)
and consequently
Hα ,n(r) =
δα ,n,k,λ
αn+ k− 1
Gαn+k(λ r). (20)
In particular with λ =
n this is a probability mixture since
δα ,n,k,√n = E ((cosΦ)
αn(1− cosΦ)k)> 0 and
αn+ k− 1
δα ,n,k,√n = 1 (21)
by virtue of the binomial series with 0≤ 1−cosΦ≤ 1−1/
n. However, for numerical
approximations different values of λ should be more suitable, e.g.
))1/(αn)
= µ1/(αn)α ,n,0 =
2Γ(αn)
(2Γ(α))n
(Γ(α/2))n
Γ(αn/2)
)1/(αn)
with µα ,n,0,λ = 1.
Besides, two further representations of Hα ,n are obtained from the following or-
thogonal series for f ∗αn from (7), (8). With
xk f ∗αn(x)dx =
Γ(αn)
µα ,n,k
αn+ k
, (23)
implied by (14), (15), and the shifted Legendre polynomials
P∗k (y) =
p∗k, jy
j, 0 ≤ y ≤ 1,
(P∗k (y))
2dy =
2k+ 1
, (24)
(A.S. 22.2.11), it follows
f ∗α ,n(x) =
cα ,n,kP
k (x/
cα ,n,k = (2k+ 1)
nαn/2
Γ(αn)
p∗k j
γα ,n, j
αn+ j
Thus, with (A.S. 11.4.26) and the modified spherical Bessel functions (A.S. 10.2.2) we
Hα ,n(r) = r
αn exp
(−1)kcα ,n,k
Ik+1/2
. (26)
Consequently, the Fourier transform of f ∗α ,n is representable by
h∗α ,n(−it) = exp
ikcα ,n,k jk
with the spherical Bessel functions jk from (A.S. 10.1). For absolutely large real t this
series is numerically more suitable than the power series from (15). Such values are
required as coefficients in the following Fourier sinus expansion for f ∗α ,n.
f ∗α ,n(x) =
bα ,n,m sin
, 0 ≤ x ≤
n , (28)
bα ,n,m =
h∗α ,n
n ∑k≥0
m+k odd
(−1)(m+k−1)/2cα ,n,k jk(mπ/2) .
Because of lemma A.2 the function f ∗αn has a square integrable derivative at least for
min(αn,n)> 2 which entails the absolute uniform convergence of the orthogonal series
in (25) and (28). Inserting (28) into (6) leads to
Hα ,n(r) = r
nr2 +m2π2
bα ,n,m
1− (−1)me−r
. (30)
Finally, with the functions
Gβ (ζ ) :=
Γ(β + k/2)
= (Γ(β ))−1M(1,β ;ζ 2)+ (Γ(β + 1/2))−1ζM(1,β + 1/2;ζ 2), β > 0,
= ζ 2(1−β )eζ
Gβ−1(ζ
2)+Gβ−1/2(ζ
, if β > 1,
the following integral representation can be derived from (11), (12), (13):
Hα ,n(r) =
βα(y−1)
G1+αn/2(ry)
y = ρeiϕ , ρ > 0.
3 The distribution of the angle Φ
Two representations of the cdf Fα ,n of U = (
ncosΦ)−1 follow directly from (25) and
(28).
Fα ,n(u) = Γ(αn)uαn−1
cα ,n,kP
= Γ(αn)uαn−1
bα ,n,m sin
Besides, an orthogonal expansion with Legendre polynomials for the density of U2 is
obtained by the moments
E (U2k) =
Γ(αn)k!
Γ(αn+ 2k)∑
Γ(α + 2ki)
Γ(α)ki!
, (32)
derived from (2).
Direct use of these orthogonal series within (5) is possible to get the cdf of S,
but the double series in theorem 4.1 is numerically more favourable. As mentioned
in the introduction it would be theoretically more satisfying to have a single alternat-
ing series or a probability mixture for this cdf. A single power series is not directly
available by a power series expansion of Gαn in (5) since the moments of cotΦ do
not all exist. However, at least for integer α , this obstacle is avoided by splitting
the domain [0,arctan(
n− 1)] of Φ by ϕn = arctan(1/
n− 1), which is the max-
imal angle ϕ for which the whole cone {Φ ≤ ϕ} is contained completely within
{x1, . . . ,xn|x1, . . . ,xn ≥ 0}. For integer α the theorem 3.1 below provides a polynomial
representation of the cdf of tanΦ restricted to 0≤ tanϕ ≤ (n−1)−1/2. This simple rep-
resentation is used within (5) to integrate over tanϕ ≤ (n−1)−1/2. The integration over
tanϕ > (n− 1)−1/2 needs only truncated moments of cotΦ = (1− cos2 Φ)−1/2 cosΦ,
obtained by the binomial expansion with the truncated moments of cosΦ, given below
in (36).
Theorem 3.1 Let α be a positive integer, then
Wα ,n(t) := Pr{tanΦ ≤ t}
Γ(αn)
(Γ(α))nnαn/2
[(α−1)n/2]
aα ,n,2 j
tn−1+2 j
n− 1+ 2 j
, 0 ≤ t ≤ (n− 1)−1/2,
where the coefficients aα ,n,2 j are given as the unique solutions of the linear equations
[(α−1)n/2]
+β n+m− j
aα ,n,2 j =
(2m)!
Γ(1/2+β +m j)/(2m j)! , α = 2β + 1
(2m+ n)!
nm+n/2
Γ(1/2+β +m j)/(2m j + 1)! , α = 2β
m, j = 0, . . . , [(α − 1)n/2].
With any ρ > 0 the sums ∑(m) within the right hand sides of (34) are also given by
e−imϕ
dϕ , α = 2β + 1,
α + 1
α + 1
e−imϕ
dϕ , α = 2β ,
After Kummer’s transformation (A.S. 13.1.27) the above confluent hypergeometric
functions M can also be expressed by Hermite polynomials due to (A.S. 13.6.17/18).
The truncated moments γ̄α ,n,k = Ē ((cosΦ)αn+k), defined by integration over
ϕ ≥ ϕn, follow from theorem 3.1 by straightforward calculation.
γ̄α ,n,k = γα ,n,k −
Γ(αn)
(Γ(α))nnαn/2
[(α−1)n/2]
aα ,n,2 jB
2 + j,
(α−1)n+k+1
2 − j;
with the incomplete beta function B(a,b;x) =
a−1(1− t)b−1dt. In particular with
α = 1 we have
γ̄1,n,k = γ1,n,k −
(n− 1)!
π (n−1)/2
Γ( n−12 )
. (37)
Two further remarks: The equations (34) can also be solved by explicit matrix in-
version. The matrix
2 +β n+m− j
, m, j = 0, . . . ,N = [(α − 1)n/2], has the
structure Γ(x)Diag
. . . ,(x)m, . . .
(x+m)N− j
with factorials (x)m = x(x+ 1) . . .(x+
m−1). The inverse (pm j(x)) of ((x+m)N− j) contains only polynomials pm j of degree
m in its mth row. These polynomials can be obtained by interpolation from the m+ 1
values
pm j(−(N − k)) =
(k = 0, . . . ,m)
(−1) j+k+1
( j+ k−m)!(N− ( j+ k))!
, m− k ≤ j ≤ N − k
0 , otherwise
The cdf F1,n of U1,n can be applied to a goodness of fit test of the hypothesis H0 : F = F0
against H1 : F 6=F0, where F0 is any specified continuous cdf. Under H0 the test statistic
∑ni=1(F0(Xi:k)−F0(Xi−1:k))2, obtained from an ordered sample Xi:k, i= 1, . . . ,k = n−1,
(X0:k =−∞, Xn:k = ∞), has the same distribution as U21,n. The asymptotic normality of
U21,n can be derived from the asymptotic bivariate normal distribution of
(Yn − n)/
(Zn − 2n)/
with correlation ρ = 2/
5. The convergence to a normal distribution
is slow and might be accelerated by a suitable transformation.
4 The distribution of the sample variance
With the functions
Hα ,n,m(r) := E (Gαn+m(r
ncosΦ))
rαn+m
Γ(αn+m)
µα ,n,m+k
αn+m+ k
(−r)k
m ∈ N0, µα ,n,k = E (
ncosφ)αn+k = n(αn+k)/2γα ,n,k from (14) we obtain
Theorem 4.1 Let X1, . . . ,Xn be independent random variables with den-
sity (Γ(α))−1e−xxα−1, α > 0, then the cdf of the sample variance S2 =
(n− 1)−1 ∑ni=1(Xi − X̄)2 is given by
(n− 1)S2 ≤ z = r2
Γ(αn+ 2 j)
Γ(αn) j!n j
Hα ,n,2 j(
z) (39)
Γ(αn)
(αn+ k)/2+ j− 1
µα ,n,2 j+k
αn+ k
(−r)k
δα ,n, j,k,λ
αn+ k− 1
Gαn+k(λ r)
with the alternating kth order differences
δα ,n, j,k,λ =
(−1)ℓ
(αn+ ℓ)/2+ j− 1
µα ,n,2 j+ℓ
λ αn+2 j+ℓ
and any λ > 0. The series in (39) is absolutely convergent with terms of the order
j−(n+1)/2
, where the O–constant depends only on α and n.
As a simple numerical example we obtain with α = 1 and n = 10 the value
Pr{S ≤ 2}= 0.98530379 . . . with at least eight correct digits.
Finally, the method, explained before theorem 3.1, is applied to obtain the alterna-
tive representations in theorem 4.2 for integer α , including in particular a probability
mixture for the exponential case α = 1. The degree of the polynomials in theorem 3.1
and the rate of convergence of the binomial series in (43) limit the numerical use of
theorem 4.2 to moderate values of αn.
The truncated moments
M̄α ,n,k = Ē
(cotΦ)αn+k
(αn+ k)/2+ j− 1
γ̄α ,n,k+2 j (43)
are used with γ̄α ,n,k from (36), obtained by integration over ϕ ≥ ϕn = arccot(
n− 1).
Theorem 4.2 With the truncated moments from (43) and (36), the pdf wα ,n and the cdf
Wα ,n of tanΦ on 0 ≤ ϕ ≤ ϕn from theorem 3.1, the distribution of S is given by each of
the following three formulas (44), (45), (46):
n− 1S ≤ r
Γ(αn)
(−1)k
γ̄α ,n,k
αn+ k
Pαn,k(r
Gαn(r
n/t)wα ,n(t)dt
with the polynomials
Pαn,k(x) =
[k/2]
(αn+ k)/2
xk−2 j
(k− 2 j)!
Γ(αn)
M̄α ,n,k
αn+ k
+ Wα ,n(1/
n− 1)Gαn
n(n− 1)
n(n−1)
gαn(x)Wα ,n
The infinite series in (45) is also given by
∆̄α ,n,k,λ
αn+ k− 1
Gαn+k(λ r) (46)
with any λ > 0 and the alternating kth order differences
∆̄α ,n,k,λ = (−∆)kM̄α ,n,0,λ =
(−1)ℓ
M̄α ,n,ℓ,λ ,
M̄α ,n,ℓ,λ = Ē
ncotΦ
)αn+ℓ
Remarks: In particular, with α = 1 the only coefficient in w1,n(t) is a1,n,0 =
(n− 1)bn−1 with the volume bn−1 = π (n−1)/2/Γ
of the (n− 1)–unit ball. Then
the last two terms in (45) are reduced to
(n− 1)!
n(n(n− 1))(n−1)/2
bn−1Gn−1(r
n(n− 1)). (47)
The truncated moments γ̄1,n,k were given in (37). Only positive coefficients ∆̄α ,n,k,λ
arise in (46) with λ =
n(n− 1) since 1 > 1− cotΦ√
n−1 ≥ 0. The last term in (45) is
bounded by Wα ,n(1/
n− 1)
1−Gαn(r
n(n− 1))
, which is often neglectible. For
actual computations the choice
Ē ((cotΦ)αn)
1−Wα ,n(1/
n− 1)
)1/(αn)
might be more favourable. With this λ the leading term in (46) becomes
1−Wα ,n
Gαn(λ r).
It would be desirable to find a similar probability mixture representation for the cdf
of S2 for all α > 0 but presumedly, also in case of its existence, the computation of the
required coefficients would not be easy.
Acknowledgement: The author would like to thank S. Kotz for his strong interest
in this work and some hints at the literature within a personal communication.
Appendix
Proof of theorem 3.1. Let α be any positive integer and ω the Lebesgue measure on
the (n− 1)–unit sphere Sn−1. Then
(x1 + . . .+ xk)
(A.1)
= 2Γ((αn+ k)/2)−1k!∑
Γ(α + k j)/2)/k j! , α + k j odd.
After an orthogonal transformation x = Ty with y1 = n
−1/2 ∑ni=1 xi = cosϑ1, followed
by transformation to polar coordinates, we obtain with c j = cosϑ j, s j = sinϑ j the same
integral value as
n+ s1τ j(ϑ2, . . . ,ϑn−1))α−1dω
= nk/2
(α−1)n
α − 1
n(ℓ−(α−1)n)/2
n−1− j
j dϑ j
(α−1)n+k−ℓ
n−2+ℓ
1 dϑ1
= nk/2
[(α−1)n/2]
aα ,n,2 jB
(α − 1)n+ k+ 1
. (A.2)
By comparison with (A.1) it follows
[(α−1)n/2]
(α − 1)n+ k+ 1
aα ,n,2 j
(A.3)
= 2n−k/2k!∑
Γ((α + k j)/2)/k j! .
With α = 2β + 1 only even k j = 2m j occur with sum k = 2m, and with α = 2β only
k j = 2m j + 1 with sum k = 2m+ n, which leads to the linear equations in (34).
If ℓ is odd then the coefficients aα ,n,ℓ vanish since
sn−1− jj dϑ j
is a linear combination of the integrals
smnn−1
M j+n−1− j
j dϑ j , (Mn−1 = 0) ,
which are different from zero only for even m2, . . . ,mn, ∑nj=2 m j = ∑
j=1 ℓ j = ℓ.
If t = tanϕ ≤ (n− 1)−1 we obtain with Sn−1,ϕ = Sn−1 ∩{ϑ1 ≤ ϕ} that
Pr{tanΦ ≤ t} = (Γ(α))−n
{ϑ1≤ϕ}
e−x j xα−1j dx j
= (Γ(α))−n
Sn−1,ϕ
ραn−1 exp(−ρ
n+ s1τ j)α−1dω
Γ(αn)
(Γ(α))nnαn/2
Sn−1,ϕ
c−αn1
n+ s1τ j)α−1dω
Γ(αn)
(Γ(α))nnαn/2
[(α−1)n/2]
aα ,n,2 j
sn−2+2 j1 c
−n−2 j
1 dϑ1 ,
where the integrals are given by (tanϕ)n−1+2 j/(n−1+2 j, which provides the asserted
result in (33). ✷
Comparing the series
(−1)k
(z/y)(αn+k)/2
αn+ k
= Γ(αn)Gαn(
with the corresponding one for Γ(αn+ 2 j)Gαn+2 j(
z/y) we obtain the identity
Γ(αn+ 2 j)
Γ(αn)
Gαn+2 j
(A.4)
by differentiation.
Lemma A.1
2αn/2
, y,z > 0. (A.5)
Proof. Using (A.4) and Cauchy’s integral formula for the derivatives of
Γ(αn+ 2 j)Gαn+2 j(
zαn/2+ j−1 exp(−
the bound in (A.5) follows from
( j− 1)!
Γ(αn+ 2 j)Gαn+2 j(
= lim
ζ αn/2+ j−1 exp(−
ζ )(ζ − z)− jdζ
yαn/2−1 exp(−
y/2)dy = 2αn/2Γ(αn)/π ,
where the way of integration is
≤ ϕ ≤
∣R ≥ y ≥−R
Lemma A.2 Let Φ denote the angle between (X1, . . . ,Xn) and (1, . . . ,1) from (4),
Wα ,n(t) the cdf of tanΦ and bn−1 = π (n−1)/2/Γ( n+12 ) the volume of the (n− 1)–unit
ball, then
Wα ,n(t)∼
Γ(αn)bn−1
(Γ(α))nn(α−1/2)n
tn−1, t ↓ 0, (A.6)
and in particular
W1,n(t) =
(n− 1)!
bn−1t
n−1, 0 ≤ t ≤ (n− 1)−1/2. (A.7)
Furthermore
E ((cosΦ)s)∼
Γ(αn)(2π)(n−1)/2
(Γ(α))nn(α−1/2)n
s−(n−1)/2, s ↑ ∞. (A.8)
Proof. Let x,xg denote the arithmetic and the geometric mean of (x1, . . . ,xn ≥ 0 and
ϕ the angle between (x1, . . . ,xn) and (1, . . . ,1). Setting xi = (1+ εi)x with any fixed x
we have the relation
(xi − x)2 =
ε2i = n tan
and for ϕ ↓ 0
xng = x
(1+ εi)∼ xn exp
∼ xn.
It follows
Wα ,n(t) = Pr{Φ ≤ ϕ}=
gα(xi)dxi
= (Γ(α))−n
n(α−1)
g exp(−nx)dx1 . . .dxn
∼ (Γ(α))−nbn−1
x(α−1)n exp(−nx)(x
n tanϕ)n−1dx
which provides (A.6).
For α = 1 this asymptotic relation can be replaced by an equation if ϕ ≤ ϕn =
arctan(1/
n− 1), which is the largest angle ϕ for which the whole cone {Φ ≤ ϕ} is
contained within {x1, . . . ,xn|x1, . . . ,xn ≥ 0}.
To prove (A.8) we obtain from (A.6)
Pr{| ln(cosΦ)| ≤ ε} = Pr{ln(1+ tan2 Φ)≤ 2ε}
∼ Wαn
∼ cα ,nε(n−1)/2, ε ↓ 0,
with a factor cα ,n determined by (A.6). For s → ∞ the relation (A.8) follows from
the Hardy–Littlewood–Karamata Tauber theorem (cf. e.g. chapter 13 in Feller 1971),
applied to
E ((cosΦ)s) = E (exp(−s| ln(cosΦ)|)) . ✷
Proof of theorem 4.1. With any ϑ ∈ (0,1) EU
(z/(U2 −ϑ/n))1/2
can be
represented by the Taylor series
With U−1 =
ncosΦ and the bounds from Lemma A.1 the terms of this series are
absolutely bounded by
2αn/2
(cosΦ)2 j
ϑ j/ j = O(ϑ j j−(n+1)/2)
due to (A.8) with s = 2 j and on O–constant depending only on α and n.
For ϑ ↑ 1 it follows from (6) and (A.4) that
Pr{(n− 1)S2 ≤ z = r2} = EU
U2 − 1/n
Γ(αn+ 2 j)
Γ(αn) j!n j
Gαn+2 j
Γ(αn+ 2 j)
Γ(αn) j!n j
Hα ,n,2 j(
with the functions Hα ,n,2 j from (38).
With any λ > 0 we obtain
Γ(αn+ 2 j)
Γ(αn) j!n j
Hα ,n,2 j(
Γ(αn)
(αn+ k)/2+ j− 1
µα ,n,2 j+k
λ αn+2 j+k
(−1)k
(λ r)αn−1+k
λ 2/n
After multiplication by exp(−λ r)∑∞ℓ=0(λ r)ℓ/ℓ! the series (41) is obtained by integra-
tion over r. ✷
Proof of theorem 4.2. With the pdf wα ,n and the cdf Wα ,n of tanΦ from theorem 3.1
and the truncated moments M̄α ,n,k = Ē
(cotΦ)αn+k
from (43) it follows with (5) that
n− 1S ≤ r} = E
Gαn(r
ncotΦ)
Gαn(r
n/t)wα ,n(t)dt +
Γ(αn)
M̄α ,n,k
αn+ k
After integration by parts and the substitution t → r
n/x the above integral becomes
n(n− 1)
n(n−1)
gαn(x)Wαn
Inserting the series for M̄α ,n,k from (43) formula (44) follows by rearranging the
resulting absolutely convergent double series according to the truncated moments
γ̄α ,n,k = Ē
(cosΦ)αn+k
from (36).
The series in (46) is obtained by the same way as (21). ✷
References
[1] Abramowitz, M. and Stegun, I. (1968) Handbook of Mathematical Functions.
Dover Publications, Inc., New York
[2] Bowman, K.O. and Shenton, L.R. (1983) The Distribution of the Standard Devi-
ation and Skewness in Random Samples from a Gamma Density — A New Look
at a Craig–Pearson Study. Oak Ridge National Laboratory/ CSD–109
[3] Box, G.E.P. (1953) Nonnormality and Tests on Variances, Biometrika 40,
318–335
[4] Craig, C.C. (1929) Sampling when the Parent Population is of Pearson’s Type III,
Biometrika 21, 287–293
[5] Feller, W. (1971) An Introduction of Probability Theory and Its Applications,
Vol. II, 2nd ed. John Wiley & Sons, New York
[6] Hyrenius, J. (1950) Distribution of Student–Fisher’s t in Samples from Com-
pound Normal Functions, Biometrika 37, 429–442
[7] Laha, R.G. (1954) On a Characterization of the Gamma Distribution, Annals of
Mathematical Statistics 25, 784–787
[8] Lukacs, E. (1955) A Characterization of the Gamma Distribution, Annals of
Mathematical Statistics 26, 319–324
[9] Mathai, A.M., and Provost, S.B. (1992) Quadratic Forms in Random Variables.
Marcel Dekker, Inc., New York
[10] Mudholkar, G.S. and Trivedi, M.C. (1981) A Gaussian Approxiamtion to the
Distribution of the Sample Variance for Nonnormal Populations, Journal of the
Americal Statistical Association 76, 479–485
[11] Pearn, W.L. and Kotz, S. (2006) Encyclopedia and Handbook of Process Capa-
bility Indices: A Comprehensive Exposition of Quality Control Measures. (Ser.
on Quality, Reliability and Engineering Statistics, Vol. 12). World Scientific Pub-
lishing Company, Singapore
[12] Pearn, W.L., Kotz, S. and Johnson, N.L. (1992) Distributional and Inferential
Properties of Process Capability Indices, Journal of Quality Technology, 24,
216–230
[13] Pearson, E.S. (1929) Note on Dr. Craig’s paper, Biometrika 21, 294–302
[14] Roy, J. and Tiku, M.L. (1962) A Laguerre Series Approximation to the Sampling
Distribution of the Variance, Sanky ā, Series A, 24, 181–184
[15] Tan, W.Y. and Wong, S.P. (1977) On the Roy–Tiku Approximation of Sample
Variances from Nonnormal Universes, Journal of the American Statistical Asso-
ciation 72, 875–880
	Introduction
	The distribution of the sum of squares
	The distribution of the angle bold0mu mumu 
	The distribution of the sample variance
