Mpmath: Asymptotic expansions

Mpmath: Asymptotic expansions#

Edgeworth expansion: general approximation of the pdf, cdf and sf#

ctx.edgeworth(x, cumulants)#

where ctx is ipm, dec, mpm, or gmp.

Calculates the Edgeworth expansion approximations of the pdf, cdf and sf by using the cumulants of the distribution.

See also Wikipedia [1182], Blinnikov and Moessner [37], Wikipedia [1176], Wikipedia [1413].

Using the first 4 cumulants, we have

\[\text{Pr}\left[\frac{u-\mu}{\sigma} \le x \right] = \Phi(x) + \phi(x)\left[\frac{\gamma_1}{6}(1-x^2) + \frac{\gamma_2}{24}(3x-x^3) + \frac{\gamma_1^2}{72}(15x-10x^3+x^5)\right] + O(N^{-3/2})\]

Withers and Nadarajah [1607] and Withers and Nadarajah [1608] have shown that the Edgeworth expansion can also be expressed in terms of standardized cumulants \(\alpha_r = \kappa_r/\kappa_2^{r/2}\) (with \(\alpha_1 = \alpha_2 = 0)\), partial exponential Bell polynomials \(B_{r,k}(\boldsymbol{\alpha}), \boldsymbol{\alpha} = (\alpha_1, \cdots, \alpha_r)\), and Hermite polynomials \(H_r(x)\):

\[f(x) = \phi(x) + \phi(x) \sum_{j=1}^{\infty} h_{1j}(x) n^{-j/2}, \quad \text{where}\]
\[h_{1j}(x) = \sum_{k=1}^j B_{r,k}(\boldsymbol{\alpha}) H_{r}(x) / r! , \quad \text{with } r = j+2k.\]
\[F(x) = \Phi(x) - \phi(x) \sum_{j=1}^{\infty} h_{0j}(x) n^{-j/2},\]
\[1 - F(x) = \Phi(-x) + \phi(x) \sum_{j=1}^{\infty} h_{0j}(x) n^{-j/2}, \quad \text{where}\]
\[h_{0j}(x) = \sum_{k=1}^j B_{r,k}(\boldsymbol{\alpha}) H_{r-1}(x) / r! , \quad \text{with } r = j+2k.\]

The partial Bell polynomials can be computed efficiently by a recurrence relation:

\[B_{r,k}(\boldsymbol{\alpha}) =\sum _{i=1}^{r-k+1}{\binom {r-1}{i-1}}\alpha_{i}B_{r-i,k-1}(\boldsymbol{\alpha}), \quad \text{where }\]
\[B_{0,0}=1; \quad B_{r,0}=0{\text{ (for }}r\geq 1); \quad B_{0,k}=0{\text{ (for }}k\geq 1); \quad B_{r,1}(\boldsymbol{\alpha})=\alpha_r.\]

The Hermite polynomials are defined by

\[H_0(x) = 1, \quad H_1(x) = x, \quad H_{r+1}(x) = x H_r(x)- r H_{r-1}(x)\]

The following code shows how to call the Edgeworth expansion for the central cdf.

def EdgeworthWhithersDemo(self):
    Order = 30
    F = 100
    X = mp.mpf(110)

    kappa = mp.matrix(Order+1, 1)
    kappa[1] = F
    for i in range(2, Order+1):
        kappa[i] = kappa[i - 1] * 2 * (i - 1)

    LeftTail1, RightTail1 = self.EdgeworthWhithers(X, Order, kappa)
    print("EdgeworthWhithers LeftTail1: ", LeftTail1, " RightTail1: ",  RightTail1)

This gives the following results:

Hello!
X:  110.0 , Z:  0.7071067811865475244008443621048490392848 ,
mean:  100.0 , sigma:  14.1421356237309504880168872420969807857
s3:  0.7602499389065232688413733269459822643682 ,
s4:  0.3106965603769277448709088912288091945547

j:  1 , s2:  0.007323188157795373284114366290203815819688
j:  2 , s2:  0.000244106271926512442803812209673460527323
j:  3 , s2:  -0.00002427501259713651514549021418419413021711
j:  4 , s2:  0.000002526386902276474842073713864467689855687
j:  5 , s2:  -0.0000002843964533935719011627505046640421102775
j:  6 , s2:  0.00000002087338535959227300224647267964356971442
j:  7 , s2:  -0.000000001797321695705282603398930576647079128807
j:  8 , s2:  0.0000000001105716283451282215450868950557968496009

LeftTail1:  0.7677952195007321931882055465876641482851
RightTail1:  0.2322047804992678068117944534123358517149

Cornish-Fisher expansion: general approximation of the qtf and isf#

ctx.cornish_fisher(x, cumulants)#

where ctx is ipm, dec, mpm, or gmp.

Calculates the Cornish-Fisher approximation of the quantile by using the cumulants of the distribution.

The quantiles of many distribution functions can be approximated by the classical Cornish and Fisher expansion Cornish and Fisher [193], Fisher and Cornish [351]), provided that the cumulants \(\kappa_r\) of the distribution are known (see also Wikipedia [1179]):

\[x = z + \sum_{k=1}^{\infty} n^{-k/2} \xi_k(z)\]

where \(\xi_k(z)\) is the collection of all the terms that belong to the \(k\)-th power of \(n^{-1/2}\). It seems to be easier to describe algorithm AS269 (see Lee and Lin [426]), which is used to calculate the \(k\)-th adjustment term \(\xi_k(z)\), by giving a program which implements it, rather than to describe it using mathematical formulas:

def CornishFisherLee(self, r, x0, ac, adj):
    a = iv.matrix(r+3, 1); d = iv.matrix(r+3, 1)
    h = iv.matrix(3 * (r+3), 1)
    p = iv.matrix((3 * (r+3)) * ((r+3) + 1) // 2, 1)
    x = iv.mpf(x0); cc = -1
    for j in range(1, r+1):
        a[j] = cc * ac[j] / ((j + 1) * (j + 2))
        cc = -cc
    h[1] = -x; h[2] = x * x - 1
    for j in range(3, (3 * r)+1):
        h[j] = -(x * h[j - 1] + (j - 1) * h[j - 2])
    d[1] = -a[1] * h[2]; adj[1] = d[1]; p[1] = d[1]
    p[3] = a[1]; fac = iv.mpf(1); ja = 0
    for j in range(2, r+1):
        fac = fac * j
        bc = iv.mpf(1)
        ja = ja + 3 * (j - 1)
        jb = ja
        for k in range(1, j):
            bcd = bc * d[k]
            bca = bc * a[k]
            jb -= 3 * (j - k)
            for m in range(1, (3 * (j - k))+1):
                jbl = jb + m; jal = ja + m
                p[jal + 1] += bcd * p[jbl]
                p[jal + k + 2] += bca * p[jbl]
            bc = (bc * (j - k)) / k
        p[ja + j + 2] += a[j]; d[j] = 0
        for m in range(2, (3 * j)+1):
            d[j] -= p[ja + m] * h[m - 1]
        p[ja + 1] = d[j]
        adj[j] = d[j] / fac

However, if one prefers mathematical formulas, Jaschke [397] gives the following description of the algorithm: “Lee developed a recurrence formula for the \(k\)-th adjustment term \(\xi_k(z)\) in the Cornish-Fisher expansion, which is implemented in the algorithm AS269 :

\[\xi_k(H) = a_k H^{*(k+1)} - \sum_{j=1}^{k-1} \frac{j}{k} \left( \xi_{k-j}(H) - \xi_{k-j} \right) * \left( \xi_j - a_j H^{*(j+1)} \right) * H, \quad \text{where } a_k = \frac{\kappa_{k+2}}{(k+2)!}.\]

\(\xi_k(H)\) is a formal polynomial expression in \(H\) with the usual algebraic relations between the summation \("+"\) and the “multiplication” \("*"\). Once \(\xi_k(H)\) is multiplied out in \(*\text{-powers}\) of \(H\), each \(H^{*k}\) is to be interpreted as the Hermite polynomial \(H_k\) and then the whole term becomes a polynomial in \(z\) with the “normal” multiplication \("\cdot"\). \(\xi_k\) denotes the scalar that results when the “normal” polynomial \(\xi_k(H)\) is evaluated at the fixed quantile \(z\), while \(\xi_k(H)\) denotes the expression in the \((+,*)\text{-algebra}\). In this formula the adjustment of a given order is calculated recursively, the values of lower order adjustments being used in each stage of calculation.”

Sheppard correction, cumulants#

ctx.sheppard(cumulants)#

where ctx is ipm, dec, mpm, or gmp.

Returns the Sheppard correction term for a vector of cumulants.

The Sheppard correction is typically applied to the cumulants computed from normally distributed data which have been binned. If \(\kappa_r^*\) is the \(r\)-th cumulant of an ungrouped distribution and \(\kappa_r\) the \(r\)-th cumulant of an grouped distribution with class interval \(c\), the corrected cumulants (under rather restrictive conditions) are

\[\begin{split}\kappa_r^*=\begin{cases} \kappa_r & \text{for } r \text{ odd} \\ \kappa_r - \frac{B_r}{r} c^r & \text{for } r \text{ even}. \end{cases}\end{split}\]

where \(B_r\) is the \(r\)-th Bernoulli number. The same correction can be applied to the cumulants of the continuity-corrected test statistic of a random variable following a lattice distribution with step \(c\). The (continuous) Edgeworth and Cornish-Fisher expansions using these corrected cumulants will then achieve the usual \(o(n^{1-r/2})\) error terms of these expansions (see Kolassa and McCullagh [418]), despite the fact that the approximated CDF is discontinuous with jumps of order \(O(n^{-1/2})\).

A similar correction can be applied to the cumulant generating function (and its derivatives), as it is needed by the Luggannani-Rice saddle point approximation.

The Sheppard correction, cumulant generating function#

ctx.sheppard_per_cgf(cgf)#

where ctx is ipm, dec, mpm, or gmp.

Returns the Sheppard correction term for cumulant generating function of a discrete (“lattice”) random variable.

The Sheppard correction is typically applied to the cumulants computed from normally distributed data which have been binned. If \(\kappa_r^*\) is the \(r\)-th cumulant of an ungrouped distribution and \(\kappa_r\) the \(r\)-th cumulant of an grouped distribution with class interval \(c\), the corrected cumulants (under rather restrictive conditions) are

\[\begin{split}\kappa_r^*=\begin{cases} \kappa_r & \text{for } r \text{ odd} \\ \kappa_r - \frac{B_r}{r} c^r & \text{for } r \text{ even}. \end{cases}\end{split}\]

where \(B_r\) is the \(r\)-th Bernoulli number. The same correction can be applied to the cumulants of the continuity-corrected test statistic of a random variable following a lattice distribution with step \(c\). The (continuous) Edgeworth and Cornish-Fisher expansions using these corrected cumulants will then achieve the usual \(o(n^{1-r/2})\) error terms of these expansions (see Kolassa and McCullagh [418]), despite the fact that the approximated CDF is discontinuous with jumps of order \(O(n^{-1/2})\).

A similar correction can be applied to the cumulant generating function (and its derivatives), as it is needed by the Luggannani-Rice saddlepoint approximation.

Luggannini-Rice expansion: general approximation of the pdf, cdf, and sf#

ctx.luggannini_rice(x, s, poly_t cgf_series)#

where ctx is ipm, dec, mpm, or gmp.

Calculates the Luggannini-Rice saddlepoint approximation of the pdf by using the cumulant generating function and its derivatives.

See also Butler [172], Daniels [202], Daniels [203], Lugannani and Rice [432], Wikipedia [1214], Wikipedia [1230].

Suppose the moment generating function of the random variable \(X\) is analytic and given by \(M(t)\) for \(t\) in some open neighbourhood of zero. Let \(K(t)\) = \(\log (M(t))\) be the cumulant generating function of \(X\), and denote by \(K^{(j)}(t)\) the \(j^{\text{th}}\) derivative of \(K(t)\).

Let \(F(x)\) and \(f(x)\) denote the CDF and pdf of \(X\), respectively, and \(\Phi(\cdot)\) and \(\phi(\cdot)\) denote the CDF and pdf of the normal distribution, respectively. Let \(s\) be the solution to the saddlepoint equation \(K^{(1)}(s)=x\), which in general needs to be solved numerically. Then

\[f(x) \thickapprox \phi(w) \sum_{k=0}^{\infty} C_k,\]
\[F(x) \thickapprox \Phi(w) + \phi(w) \sum_{k=0}^{\infty} \left( A_k - B_k \right), \quad \text{and}\]
\[1 - F(x) \thickapprox \Phi(-w) - \phi(w) \sum_{k=0}^{\infty} \left( A_k - B_k \right), \quad \text{where }\]
\[w = \text{sgn}(s) \sqrt{2 (s K^{(1)}(s) - K(s))}; \quad t = \sqrt{K^{(2)}(s)}; \quad u = s t.\]

The coefficients \(A_k\), \(B_k\) and \(C_k\) are calculated as follows:

\[A_0 = \frac{1}{u}, \quad A_j = A_0 \sum_{n=0}^{2j} \left(\frac{1}{u} \right)^{2j-n} \sum_{m=0}^{n} d_{m,n} \left(-2 \right)^{m+j} \left(\tfrac{1}{2} \right)_{m+j},\]

where \(d_{m,n}\) is computed (with \(d_{0,0}=1\) and \(d_{0,n}=0\) for \(n>0\)) from

\[d_{m+1,n+1} = \frac{1}{n+1} \sum_{j=1}^{n-m+1} \frac{j}{(j+2)!} \frac{K^{(j+2)}(s)}{ t^{j+2}} \: d_{m,n-j+1}, \quad 0 \le m \le n.\]
\[B_0 = \frac{1}{w}, \quad B_j =\frac{-2 (j-\tfrac{1}{2}) B_{j-1}}{w^2}.\]
\[C_0 = \frac{1}{t}, \quad C_j = \sum_{m=1}^{2j} d_{m,2j} \left(-2 \right)^{m+j} \left(\tfrac{1}{2} \right)_{m+j}.\]

Jensen expansion: general approximation of the qtf and isf#

ctx.jensen_inv(q, s, x0, poly_t cgf_series)#

where ctx is ipm, dec, mpm, or gmp.

Calculates the Jensen approximation of the qtf and isf

Using the same notation as above, Jensen [398] derived the following approximation (see also Wang [865]):

\[F(x) \thickapprox \Phi(r(s)), \quad \text{where}\]
\[w(s) = \text{sgn}(s) \sqrt{2 (s K'(s) - K(s))},\]
\[u(s) = s \sqrt{K''(s)} ,\]
\[r(s) = w(s) + \frac{v(s)}{w(s)} , \quad \text{where } v(s)=\log \frac{u(s)}{w(s)}.\]

To obtain an approximation to the \(\alpha\)-quantile of \(F(x)\), we start with a reasonable estimate, say, \(x_0\), and convert this to the corresponding saddlepoint, say \(s_0\), such that \(K'(s_0)=x_0\).

Let \(z_\alpha\) be the \(\alpha\)-quantile of the standard normal distribution. We wish to solve for \(s\) in the equation

\[h(s) = w(s) + \frac{v(s)}{w(s)} - z_\alpha = 0.\]

Starting with \(n=0\), we compute successive approximations using Newton iterations as

\[s_{n+1} = s_n - \frac{h(s_n)}{h'(s_n)}, \quad \text{where}\]
\[h'(s) = w'(s) + \frac{w(s) u'(s) - u(s) w'(s)\left(v(s)+1\right)}{u(s) w(s)^2},\]
\[w'(s) = \text{sgn}(s) \frac{s K''(s)}{\sqrt{2 s K'(s) - 2 K(s)}} = \text{sgn}(s) \frac{s K''(s)}{w(s)},\]
\[u'(s) = \frac{s K^{(3)}(s) + 2 K''(s)}{2 \sqrt{K''(s)}}.\]

The final estimate \(x_{n+1}\) of the \(\alpha\)-quantile of \(F(x)\) is then obtained as \(x_{n+1} = K'(s_{n+1})\).

Box-Davis expansion: general approximation of the pdf, cdf and sf#

ctx.box_davis_expansion(x, f, rho, omega)#

where ctx is ipm, dec, mpm, or gmp.

Calculates the Box-Davis approximation of the pdf, cdf and sf

The Box-Davis expansion supports a general class of random variables, say \(X\), with absolutely continuous distribution functions, for which the logarithm of the characteristic function, say \(\log \left[ C_X(t) \right]\), may be validly represented by an asymptotic series

\[\log \left[ C_X(t) \right] \sim -\frac{1}{2} f \log(1-2it) + \sum_{r=1}^{\infty} \omega_r(\rho) \left[(1-2it)^{-r} -1\right],\]

corresponding to a limiting chi-squared distribution with \(f\) degrees of freedom, where \(\rho\) is chosen in such a way that \(\omega_1 = 0\), whenever this is possible.

The first comprehensive review of this type of expansion for the CDF is due to Box (see Box [40]). The corresponding expansion for the inverse CDF is due to Davis (see Davis [211]), who also gave algorithms for variables with characteristic functions different from those considered by Box. We refer therefore to this system of asymptotic expansions as the Box-Davis expansion. Distributions covered by this expansion include the beta distribution, the distribution of the product of independent beta variables, Wilks’ Lambda, Mauchley’s test for sphericity, Bartletts test for equality of independent covariance matrices, Hotelling’s \(T^2\) and Pillai’s \(V\). For more information, see Anderson [9] and Ginzberg [364].

Let \(\theta\) be a suitably chosen constant (which will usually be proportional to sample size). We define

\[a_0 = 1; \quad a_j = \frac{1}{j} \sum_{l=1}^j l \omega_l a_{j-l}, \quad j>0.\]
\[\log(K_B) = -\sum_{j=1}^{r} \omega_j + O(\theta^{-r-1}),\]

Then \(X\) follows asymptotically (with \(\theta \rightarrow \infty\)) a \(\chi^2\)-distribution with \(f\) degrees of freedom, and the pdf, CDF of \(M\) can be developed into asymptotic expansions including terms of order \(O(\theta^{-r})\), depending only on \(f\), \(\rho, K_B\) and \(a_j\):

\[\text{pdf}_X(x)) = \rho K_B \sum_{j=0}^{r} a_j f_{\chi_{f+2j}^2}(\rho x) + O(\theta^{-r-1})\]
\[\text{cdf}_X(x) = K_B \sum_{j=0}^{r} a_j F_{\chi_{f+2j}^2}(\rho x) + O(\theta^{-r-1})\]
\[1-\text{cdf}_X(x) = K_B \sum_{j=0}^{r} a_j \left[1 - F_{\chi_{f+2j}^2}(\rho x) \right] + O(\theta^{-r-1})\]

Gupta and Tang (see Gupta and Tang [372]) show that the above series are convergent as \(r \rightarrow \infty\) for fixed \(x < 4\pi \theta\).

When calling the following functions, only \(f\), \(\rho\) and \(\omega_j\) need to be supplied by the user. \(K_B\) and \(a_j\) are calculated internally.

Box-Davis expansion: general approximation of the qtf and isf (for \(\omega_1 = 0\))#

ctx.box_davis_expansion_inv(x, f, rho, omega)#

where ctx is ipm, dec, mpm, or gmp.

The Box-Davis expansion can be inverted efficiently in a two-step procedure. In the first step, an initial estimate is obtained, using the procedure given below. In the second step, the initial estimate is used as starting point of a Newton iteration, using the functions of the previous section.

Assuming \(\omega_1 = 0\), an approximation of the quantile function of \(X\), \(\text{qtf}_X(q))\), in terms of the quantiles of a \(\chi^2\)-distribution with \(f\) degrees of freedom, say \(u = \chi_{f,\alpha}^2\), can be given as \(\text{qtf}_X(q)) = (u + 2S)/\rho + O(\theta^{-8})\), where \(S\) is a polynomial in \(u\), whose coefficients depend only on \(f\) and \(\omega_j\) (see Davis [211]).

\begin{eqnarray} S & = & \omega_2 P_2(u) \nonumber \\ \nonumber & + & \omega_3 P_3(u) \\ \nonumber & + & \omega_4 P_4(u) + \tfrac{1}{2} \omega_2^2P_{2,2}(u) \\ \nonumber & + & \omega_5 P_5(u) + \omega_3 \omega_2P_{3,2}(u) \\ \nonumber & + & \omega_6 P_6(u) + \omega_4 \omega_2P_{4,2}(u) + \tfrac{1}{2} \omega_3^2P_{3,3}(u) + \tfrac{1}{6} \omega_2^3P_{2,2,2}(u) \\ \nonumber & + & \omega_7 P_7(u) + \omega_5 \omega_2P_{5,2}(u) + \omega_4 \omega_3P_{4,3}(u) + \tfrac{1}{2} \omega_3 \omega_2^2P_{3,2,2}(u), \nonumber \end{eqnarray}
\[P_r(u) = \frac{u}{f_1} + \ldots + \frac{u^r}{f_r}, \quad f_r = f(f+2) \ldots (f+2r-2), \nonumber\]
\begin{eqnarray} P_{2,2} & = & - \frac{8 u^4 (f+3)}{(f_2 f_4)} + \frac{8 u^3}{(f_2 f_3)} + \frac{6 u^2}{(f f_2)} + \frac{2 u}{f^2} \nonumber \end{eqnarray}
\begin{eqnarray} P_{3,2} & = & - \frac{12 u^5 (f+4)}{(f_2 f_5)} - \frac{2u^4 (f-6)}{(f_2 f_4)} + \frac{2 u^3 (3 f+10)}{(f_2 f_3)} + \frac{6 u^2}{(f f_2)} + \frac{2 u}{f^2} \nonumber \end{eqnarray}
\begin{eqnarray} P_{4,2}(u) & = & -\frac{16 u^6 (f+5)}{f_2 f_6} -\frac{4 u^5 (f-4}{f_2 f_5} ) + \frac{2 u^4 (3 f+14)}{f_2 f_4} +\frac{2 u^3 (3 f+10)}{f_2 f_3} + \frac{6 u^2}{f f_2} + \frac{2 u}{f^2}, \nonumber \end{eqnarray}
\begin{eqnarray} P_{3,3}(u) & = & -\frac{6 u^6 (3 f^2+30 f+80)}{f_3 f_6} - \frac{6 u^5 (f_2+2 f-16)}{f_3 f_5} + \frac{4 u^4 (f+12)}{f_2 f_4} + \frac{4 u^3 (3 f+8)}{f_2 f_3} + \frac{6 u^2}{f f_2} + \frac{2 u}{f^2}, \nonumber \end{eqnarray}
\begin{eqnarray} P_{2,2,2}(u) & = & \frac{32 u^6 (7 f^2+62 f+120)}{f_2^2 f_6} - \frac{32 u^5 (2 f^2+37 f+96)}{f_2^2 f_5} - \frac{8 u^4 (23 f^2+124 f+132)}{f_2^2 f_4} \nonumber \\ \nonumber & \: & - \frac{8 u^3 (f-10)}{f f_2 f_3} + \frac{28 u^2}{f^2 f_2} + \frac{4 u}{f^3}, \nonumber \end{eqnarray}
\begin{eqnarray} P_{5,2} & = & - \frac{20 u^7 (f+6)}{(f_2 f_7)} - \frac{2 u^6 (3 f-10)}{(f_2 f_6)} +\frac{2 u}{f^2} + 2 \sum_{r=2}^5 \frac{u^r (3f+4r-2)}{f_2 f_r} \nonumber \end{eqnarray}
\begin{eqnarray} P_{4,3}(u) & = & - \frac{24 u^7 (f_2+12 f+40)}{f_3 f_7} - \frac{2 u^6 (5 f_2+18 f-80)}{f_3 f_6} + \frac{2 u^5 (f_2+42 f+176)}{f_3 f_5} \nonumber \\ \nonumber & \: & +\frac{4 u^4 (3 f+16)}{f_2 f_4} + \frac{4 u^3 (3 f+8)}{f_2 f_3} + \frac{6 u^2}{f f_2} + \frac{2 u}{f^2}, \end{eqnarray}
\begin{eqnarray} P_{3,2,2}(u) & = & \frac{192 u^7 (2 f^3+31 f^2+154 f+240)}{f_2 f_3 f_7} - \frac{16 u^6 (4 f^3+153 f^2+1106 f+2160)}{f_2 f_3 f_6} \nonumber \\ \nonumber & \: & - \frac{8 u^5 (35 f^3420 f^2+1540 f+1632)}{f_2 f_3 f_5} - \frac{4 u^4 (25 f^2+80 f+12)}{f_2^2 f_4} + \frac{4 u^3 (7 f+38)}{f f_2 f_3} + \frac{28 u^2}{f^2 f_3} + \frac{4 u}{f^3}. \\ \nonumber \end{eqnarray}

When calling the following functions, only \(f\), \(\rho\) and \(\omega_j\) need to be supplied by the user. \(K_B\) and \(a_j\) are calculated internally.

Contents