1 Introduction
Even though the Edgeworth expansion is often, at the first glance, mostly viewed as a theoretical tool, it has steadily gained interest from the practitioner point of view for several decades. At the time present, there are too many examples of smart applications of these expansions to list them all. For instance, it has successfully been used to improve confidence intervals, either asymptotic or bootstrap ones, giving sharp insights on, among other properties, their coverage or expected length. Hall Hall (1992) and Shao and Tu Shao and Tu (1996) provide a detailed account on this topic. Moreover, the Edgeworth expansion is deeply rooted in one of the more recent outstanding study dealing with confidence interval for a binomial proportion (Brown, Cai, and Das Gupta, Brown et al. (2001, 2002)), in an exponential family (Brown, Cai, and Das Gupta, Brown et al. (2003)) and in discrete distributions (Cai, Cai (2005)).
The major drawback of these expansions are their long and tedious computations. Even when dealing with the most simple statistic the studentized mean, a famous textbook can provide an incorrect formula as the secondorder CornishFisher expansion stated in Theorem 13.5 by Das Gupta in (Das Gupta, 2008, p. 194)
. Make a slight change in the denominator of the usual raw variance estimator to use instead its unbiased counterpart by simply replacing
by to get an exactly distributed statistic in the Gaussian framework and it leads to different secondorder formulas, which accounts for the discrepancies of the Edgeworth expansions for the studentized mean when compared to the Edgeworth expansion of the distribution, for instance as stated in Theorem 11.4.2 by Lehman and Romano in (Lehman and Romano, 2005, p. 460).For a smooth or regular enough model of the mean of an iid sample drawn from an absolutely continuous distribution , first and secondorder Edgeworth expansions are valid. Closed, yet complex, formulas were derived by Withers Withers (1983) and Hall Hall (1992) to compute them from the partial derivatives of
and the central moments of the underlying distribution
. In addition, in such a model that exhibits a valid secondorder Edgeworth expansion, there exists a closed formula for the acceleration constant used in the bootstrap biascorrected and accelerated technique (BCA). Were we able to evaluate this formula, this would help to remedy to “A disadvantage of the bootstrap BCA”, as stated Shao and Tu (Shao and Tu, 1996, p. 140), which is “that the determination of is not easy, especially when the problem under consideration is complex”. Moreover, in some statistical models, as put forward in two cases in Sections 4 and 5, the acceleration constant is free of the unknown parameters of the underlying distribution and thus can be derived exactly.All these reasons account for our use of the Maple computation engine Maple Team (2014), that combines numeric computations with symbolic capabilities, in order to design a tool able to automatically and trustworthy compute Edgeworth expansions as well as CornishFisher ones or the acceleration constant used in BCA bootstrap technique. This is this tool, that we provide as a Maple worksheet.
Maple functions or expressions code can be translated into MATLAB code, which in our setting, is virtually R R Core Team (2014) code. The R language can then be used for display purposes, as we did with Figures 1 and 2, to apply bootstrap techniques of the boot package Canty and Ripley (2014); Davison and Hinkley (1997)
or for any further statistical analysis, for instance to deal with the nonincreasing cumulative distribution functions that may result from the above Edgeworth expansions using the rearrangement operator
Hardy et al. (1952) implemented in the Rearrangement package Chernozhukov et al. (2010); Graybill et al. (2011). Using Withers (1983), it is easy to derive CornishFisher expansions of any higher order and upgrade our worksheet to have Maple compute them.Previous work by Finner and Dickhaus Finner and Dickhaus (2010) only addresses the part on the computation of Edgeworth expansions up to any order. From a (applied) statistical point of view, it is not very useful. On the contrary, we not only carry out the Edgeworth expansions but also the CornishFisher ones, BCA acceleration constant derivation, export to and import in R of the results, then use of the R language to derive increasing rearrangements of the computed expansions, which sums up to a significant additional value of this article for statisticians and practitioners.
2 EdgeworthCornishFisher expansions and BCA bootstrap
2.1 Introducing the Edgeworth expansion
Let
be a sequence of random variables with 0 mean and unit variance. The cumulant generating function
of is defined by: , . The ^{th} cumulant of is given by: , . Consequently, we show that:(1) 
The ^{th}
derivative of the characteristic function of
exists at the point if and only if . Then the ^{th} cumulant of exists if and only if . The entire series expansion and the identification of the coefficients lead to the expression of the first four cumulants, see Hall (1992) for more details:Assume that there exist four real numbers , , , so that the first four cumulants satisfy the following asymptotic expansions:
Then Equation 1 becomes:
A Fourier inversion of this expansion leads to:
uniformly in where
is the probability density function of a standard Gaussian distribution and
(2)  
(3) 
2.2 A model valid for Edgeworth expansions
Let , , , …be independent and identically distributed random column
vectors with mean
denoted by , and . Let be a smooth function satisfying . Denote the th element of a vector by: or ,(4)  
(5) 
We will focus on functions such as:

, where is the (scalar) parameter estimated by and is the asymptotic variance of ;

, where is an estimate of .
We can assume that is a known function; as stated in (Hall, 1992), the following function fulfills the previous requirements: . If used, stands for the nonstudentized case and the studentized one. This “smooth function model” (Hall, 1992) allows to study problems where is a mean, or a variance, or a ratio of means or variances, or a difference of means or variances, or a correlation coefficient, etc.
Example 2.1.
Let be a random sample from an univariate population with mean and variance , and if we wished to estimate , then we would take , , , . This would ensure that , (the sample mean), , and
Example 2.2.
Let be a random sample from an univariate population with mean and variance , and if we wished to estimate , then we would take , , , , . In this case, we have:
The cases where is a correlation coefficient (a function of five means), or a variance ratio (a function of four means), among others, may be treated similarly.
Since and , by a Taylor expansion: , , where means that and
Using Hall (1992) and Withers (1983), it is easy to derive Edgeworth expansions of any higher order and complete our worksheet to have Maple compute them. One can also try to use Finner and Dickhaus work Finner and Dickhaus (2010).
Remark 2.4.
Edgeworth expansions can also be used to derive asymptotic expansions for cumulative density functions of distribution. A case in point is the distribution, see Lehman and Romano in Formulas 11.73 and 11.74 of Theorem 11.4.2, (Lehman and Romano, 2005, p. 460) for the secondorder Edgeworth expansion of the distribution. An easy way to retrieve an exactly distributed statistic in a iid Gaussian framework, is to consider with the square root of the unbiased variance estimator (). Unfortunately, when applying Theorem 2.3 to the sample mean, we get the expansion of with the square root of the raw variance estimator (). This subtle change in the normalizing constant leads to discrepancies between the two secondorder terms of the Edgeworth expansion of the distribution and of the studentized mean. Using asymptotic expansions, we show in the provided Maple worksheet, and report in Section 3.4, how to switch from the expansion of the studentized mean to the one of the distribution.
2.3 CornishFisher expansion
The CornishFisher expansion is the expansion of the
level quantile of an estimator from its Edgeworth expansion:
Denote by the level quantile of a standard Gaussian distribution. Equation 6 leads to an expansion of as a polynomial in :
Even though the true values of the parameter , whose value may be required to compute the polynomials and in the CornishFisher expansion when working with non pivotal statistics, are unknown, Hall (1992) and Shao and Tu (1996) showed that this issue can be solved using bootstrap techniques to ensure the interval still features the secondorder probability matching property.
The same remark as above applies, using Withers (1983), it is easy to derive CornishFisher expansions of any higher order and complete our worksheet to have Maple compute them.
2.4 BCA bootstrap
In a framework valid for secondorder Edgeworth expansions, bootstrap bias corrected and accelerated (BCA) are secondorder probability matching intervals (Hall, 1992; Shao and Tu, 1996). Yet, even though the simpler bootstrap intervals also feature secondorder probability matching, BCA ones are often preferred (Hall, 1992, p. 93 and pp. 133134) and their definition will be now shortly recalled. Note that the BCA intervals are transformation respecting and thus invariant under reparametrization.
Let be a random sample from the population studied and a parameter of interest. Define for the bootstrap distribution of and, for all in , . In this notation, . The sequences of quantiles would be “centered” empirically if were equal to . In general, this identity does not hold, and we might wish to correct for centering error. Let be the estimate of the centering correction or bias correction. An onesided BCA interval for , with nominal coverage , is given by where:
The quantity
amounts for a skewness correction. An explicit formula for the acceleration constant
, valid for nonparametric as well as parametric bootstrap, is worked out by Hall (Hall, 1992, p. 132):An estimate is also proposed by Hall (Hall, 1992, p. 133):
A twosided BCA interval is given by .
3 Expanding the mean
Let be the skewness and
the excess kurtosis. All the explicit results in this section were derived from the theoretical ones stated in Section
2 using Maple. For this first application we will follow closely the the steps detailed in the file SI_BertrandMaumy_Edgeworth_mean.mw.3.1 Setting the framework
We begin with a restart statement to be sure to start with a fresh session.
> restart;
We first need to set within Maple the framework of the Exemple 2.1 in Section 2.2. We want to highlight that the user only needs to define the function, see Section 2.2, all the other quantities required for the computation, for both the nonstudentized or studentized cases, been automatically derived from this very definition as soon as the previous function remained called .
> g:=proc(x1) return x1; end proc;
Automatically retrieve the value of the dimension , see Section 2.2, for the nonstudentized case (Dim:=1) and for the studentized one (Dims:=2). Enlarge the vector to take into account higher moments required for the computations related to the studentized statistic.
> Dim:=nops([op(1,op(1,g))]);Dims:=2*Dim;NtoL:=[evaln(x(1 .. Dims))];
Even the function is automatically defined by this code chunk.
> h := proc (x::(seq(name))) option remember; local NNN, derg; Ψfor NNN in seq(1..(_npassed/2)) do ΨΨderg[NNN]:=diff(g(op([1..Dim],[_passed])), ΨΨ[op(NNN,[_passed])]); Ψend do; return sqrt(sum(sum(derg[iii]*derg[jjj]*(op(iii+jjj,[_passed])op(iii, [_passed])*op(jjj,[_passed])),jjj=1..(_npassed/2)),iii=1..(_npassed/2))); end proc;
We begin with the link between raw moments, from first to fourth order, and centered ones. We assume is the skewness, the excess kurtosis and , , standardized central moments. Can be tuned with a specific distribution. The variable names and had to be used to avoid errors due to confusion with two builtin R functions Gamma and kappa if the result is exported to the R statistical language.
> EspX:=[]:for ii in [seq(1..nops(NtoL))] do EspX:=[op(EspX),cat(x,ii)=expand(sum(binomial(ii, j)*’muj’* sigma^j*mu^(iij), ’j’ = 0 .. ii))]: end do: EspX;
We reflect some known facts about the raw moments so that Maple can use them when simplifying complex expressions.
> for AA in [seq(1..nops(NtoL))] do if ‘mod‘(AA,2)=1 then assume(NtoL[AA],real) else assume(NtoL[AA]>0) end if end do; NtoL;
Several code chunks, not shown, now define various quantities of interests either for the computation itself of for speeding it up, namely: , , its partial derivatives , , its partial derivatives , see Section 2.2. The last preparation step, is to compute all the moments that appear in Section 2.2’s formulas. We only report the first and the last of the five code chunks required to carry out this work.
> for ii from 1 to 4*Dims do XXii:=t^iiexpand(subs([mu4 = kappa1+3, mu3 = Gamma1, mu2 = 1, mu1 = 0, mu0 = 1], sum(binomial(ii, j)* ’muj’*sigma^j*mu^(iij), ’j’ = 0 .. ii))): end do:
> for ii from 1 to Dims do for jj from 1 to Dims do for kk from 1 to Dims do for ll from 1 to Dims do Muii"."jj"."kk"."ll:=simplify( subs([seq(t^(4*Dimsi)=t^(4*Dimsi)XX(4*Dimsi),i=0..(4*Dims1))], expand(XXii*XXjj*XXkk*XXll))); end do; end do; end do; end do;
> ExpGaussian:=[Gamma1=0,kappa1=0]: for ii from 5 to 4*Dims do ExpGaussian:= [op(ExpGaussian),muii=(‘mod‘(ii+1, 2))*doublefactorial(ii1)]: end do: ExpGaussian;
3.2 Acceleration constant
The purpose of the following code chunk is to compute involved in the derivation of the value of the acceleration constant , see Section 2.4.
> A:=simplify(add(add(add(PartAs([NtoL[i]])*PartAs([NtoL[j]])* PartAs([NtoL[k]])*Mui"."j"."k,i=1..Dims),j=1..Dims),k=1..Dims)): A:=simplify(subs(EspX,A),assume=positive);
As a consequence, we can state that the value of the acceleration constant is :
3.3 Nonstudentized statistic
We now compute the values , , and defined in Theorem 2.3.
> k12:=simplify(1/2*(add(add(PartA0([NtoL[i],NtoL[j]])*Mui"."j,i= 1..Dim),j=1..Dim))): > k12:=simplify(subs(EspX,k12));
> k31:=k31:=simplify(add(add(add(PartA0([NtoL[i]])*PartA0([NtoL[j]])*PartA0( [NtoL[k]])*Mui"."j"."k,i=1..Dim),j=1..Dim),k=1..Dim)+3*add(add( add(add(PartA0([NtoL[i]])*PartA0([NtoL[j]])*PartA0([NtoL[k],NtoL[l]])* Mui"."k*Muj"."l,i=1..Dim),j=1..Dim),k=1..Dim),l=1..Dim)): > k31:=simplify(subs(EspX,k31));
It is easy to derive results for a given distribution by replacing moments with their actual values. For instance, in the Gaussian case and . In the following, we will report only the results for the general case.
> eval(k31,[Gamma1=0,kappa1=0]);
> k22:=simplify(add(add(add(PartA0([NtoL[i]])*PartA0([NtoL[j],NtoL[k]])* Mui"."j"."k,i=1..Dim),j=1..Dim),k=1..Dim)+1/2*add(add(add(add( PartA0([NtoL[i],NtoL[j]])*PartA0([NtoL[k],NtoL[l]])*Mui"."k* Muj"."l,i=1..Dim),j=1..Dim),k=1..Dim),l=1..Dim)+add(add(add(add( PartA0([NtoL[i]])*PartA0([NtoL[j],NtoL[k],NtoL[l]])*Mui"."j* Muk"."l,i=1..Dim),j=1..Dim),k=1..Dim),l=1..Dim)): > k22:=subs(EspX,k22);
> k41:=simplify(add(add(add(add(PartA0([NtoL[i]])*PartA0([NtoL[j]])*PartA0( [NtoL[k]])*PartA0([NtoL[l]])*(Mui"."j"."k"."l3* Mui"."j*Muk"."l),i=1..Dim),j=1..Dim),k=1..Dim),l=1..Dim)+12* add(add(add(add(add(PartA0([NtoL[i]])*PartA0([NtoL[j]])*PartA0([NtoL[k]])* PartA0([NtoL[l],NtoL[m]])*Mui"."l*Muj"."k"."m,i=1..Dim), j=1..Dim),k=1..Dim),l=1..Dim),m=1..Dim)+12*add(add(add(add(add(add(PartA0( [NtoL[i]])*PartA0([NtoL[j]])*PartA0([NtoL[k],NtoL[l]])*PartA0([NtoL[m], NtoL[o]])*(Mui"."k*Muj"."m*Mul"."o),i=1..Dim),j=1..Dim), k=1..Dim),l=1..Dim),m=1..Dim),o=1..Dim)+2/3*add(add(add(add(add(add(PartA0( [NtoL[i]])*PartA0([NtoL[j]])*PartA0([NtoL[k]])*PartA0([NtoL[l],NtoL[m], NtoL[o]])*(6*Mui"."m*Muk"."l*Muj"."o),i=1..Dim), j=1..Dim),k=1..Dim),l=1..Dim),m=1..Dim),o=1..Dim)): > k41:=subs(EspX,k41);
In order to compute and , we now evaluate Equations 2 and 3.
> P1:=(k12+1/6*k31*(x^21)): > P1:=simplify(P1);
> P2:=x*(1/2*(k22+k12^2)+1/24*(k41+4*k12*k31)*(x^23)+1/72*k31^2* (x^410*x^2+15)): > P2:=collect(collect(collect(P2,x),Gamma1),kappa1);
> P11:=P1: > P11:=collect(collect(collect(simplify(subs(EspX,P11)),x),Gamma1),kappa1);
> P21:=P1*diff(P1,x)P2x*P1*P1/2: > P21:=collect(collect(collect(simplify(subs(EspX,P21)),x),Gamma1),kappa1);
3.4 Studentized statistic
Same value as Shao and Tu (Shao and Tu, 1996, p. 145).
We now provide the expansion of the distribution, cf. remark 2.4. See the worksheet on how we used the polynomials and and an asymptotic expansion of the standard Gaussian CDF, density and polynomials and to get the polynomials and of the secondorder Edgeworth expansion of the distribution.
These values are the one given by Lehman and Romano in Formulas 11.73 and 11.74 of Theorem 11.4.2, (Lehman and Romano, 2005, p. 460) for the secondorder Edgeworth expansion of the distribution.
Corrects the value of given in Theorem 13.5 in (Das Gupta, 2008, p. 194). Same value as the one found by Shao and Tu in (Shao and Tu, 1996, p. 145) and Finner and Dickhaus in Finner and Dickhaus (2010).
3.5 Export results from Maple to R
> fd := fopen("A_mean.R", WRITE, TEXT);
> cg1 := CodeGeneration[’Matlab’](A, resultname = ’A’, output = string);
> fprintf(fd, cg1);
> fclose(fd);
4 Expanding the variance
All the explicit results in this section were derived from the theoretical ones stated in Section 2 using Maple. The steps and details of the computations are available in the file SI_BertrandMaumy_Edgeworth_variance.mw. To demonstrate how is to simple to switch from the mean expansion in Section 3 setting to another one, namely the variance expansion one, the only code chunk from that must be modified is the following one
> g:=proc(x1,x2) return x2x1^2; end proc;
4.1 Acceleration constant
In the Gaussian case, we find the same value as Shao and Tu in (Shao and Tu, 1996, p. 138) by substituting with the well known values .
4.2 Nonstudentized statistic
Same as the value provided by Hall in (Hall, 1992, p. 75 and 76).
The values of , , and can be found in the dedicated Maple worksheet.
4.3 Studentized statistic
Same value as the one found by Hall in (Hall, 1992, p. 76).
The values of , , and can be found in the dedicated Maple worksheet.
5 A more complex example
All the explicit results in this section were derived from the theoretical ones stated in Section 2 using Maple. The steps and details of the computations are available in the file SI_BertrandMaumy_Edgeworth_ML_normal_std_sym.mw.
5.1 Introduction
One burning issue that commonly arises in official quality laboratories or in the industries of various sectors, namely chemistry, pharmacy, food processing, etc, is the estimation of the proportion of experimental results that lie within two limits. For instance, in official quality laboratories, is the the range of the acceptable values and is the proportion of valid units.
Let , where , be a sequence of independent experimental results that follow a Gaussian distribution with mean and variance . We denote the sample mean and the sample variance respectively by:
We introduce the proportion defined by:
(7) 
where stands for the cumulative distribution function of the standard Gaussian distribution, and , with , are the “limits of acceptance”, whose values depend on the active regulatory norm.
Using Equation 7, one can put forward as the maximum likelihood (ML) point estimator of :
Comments
There are no comments yet.