asymptotic analysis of the bootstrap methods 20030721

Hidetoshi Shimodaira 

Department of Mathematical and Computing Sciences 
Tokyo Institute of Technology 
Ookayama, Meguro, Tokyo 152-8552, JAPAN. 

EMAIL: shimo@is.titech.ac.jp 
URL: http://www.is.titech.ac.jp/~shimo/

Asymptotic Analysis of the Bootstrap Methods

This documentation consists of three parts, namely "Exponential Family of Distributions", "Tube-Coordinates and (Subscript[z, c])-formula", and "Bootstrap Methods".  Each part is an independent Mathematica session, and should be run separately.

The calculation is third-order accurate, correct  up to O((n^(-1))) terms ignoring the error of O((n^(-3/2))). We use "o" to indicate O((n^(-1/2))) quantities. So, (o^2) indicates O((n^(-1))) and (o^3) indicates O((n^(-3/2))), etc. We repeatedly use a simplification function to keep only up to (o^2) terms.

In the first two parts, the tensor notation is heavily used. The add-on package "MathTensor" is required to run the Mathematica session by yourself.

Exponential Family of Distributions

In this part, we will provide a canonical form the density function of the exponential family of distributions. First, some basic results for normal density is obtained. Then, the distribution function of the exponential family of distributions specified in the natural parameter will be transformed to an expression using the expectation parameter and the derivatives of the potential function.

Startup

This section initializes the Mathematica session.

packages

<<  Statistics`ContinuousDistributions`

<< MathTens.m (Windows)

error messages

Off[General::spell1]

Off[General::spell]

distribution functions

gammadist[x_, m_, α_] := 
      PDF[GammaDistribution[m, α], x]
Gammadist[x_, m_, α_] := 
      CDF[GammaDistribution[m, α], x]
f[x_] := PDF[NormalDistribution[0, 1], x]
F[x_] := CDF[NormalDistribution[0, 1], x]
Q[x_] := Quantile[NormalDistribution[0, 1], x]
Chidist[x_, {di_, nc_}] := 
      CDF[NoncentralChiSquareDistribution[di, nc], x]

Normal distribution

This section first calculates the moments of the normal variables, which will be used to calculate the expected value of the exponential of normal variables.

the moments of the multivariate normal distribution

We consider the multivariate normal random vector x = ((Subscript[x, 1]),  ... , (Subscript[x, dim])) of dim-dimensions with mean b = ((Subscript[b, 1]),  ... , (Subscript[b, dim])), and the identity covariance matrix. The density function is f(x) = f[(Subscript[x, 1]) - (Subscript[b, 1])] f[(Subscript[x, 2]) - (Subscript[b, 2])] ...  
          f[(Subscript[x, dim]) - (Subscript[b, dim])]. In this subsection, we define "ruleintx" to calculate the central moments of x, such as E((Subscript[x, a]) (Subscript[x, b])), E((Subscript[x, a]) (Subscript[x, b]) (Subscript[x, c]) (Subscript[x, d])), and E((Subscript[x, a]) (Subscript[x, b]) (Subscript[x, c]) (Subscript[x, d]) (Subscript[x, e]) 
          (Subscript[x, f])) for b = 0.

the central moments of the standard normal variable in one dimension

\(The\ following\ intx2f[
          n] = \ \(E \((x\^\(2  n\))\) = ∫\_\(-∞\)\%∞\( x\^\(2  n\)\) 
            f[x] dx\ gives\ \(\(\(\((2  
                        n)\)!\)\/\(\(2\^n\) \(n!\)\)\)\(.\)\)\)\)

intx2f[n_] = 
      FullSimplify[(∫_(-∞)^∞) (x^(2  n)) 
            f[x] \[DifferentialD]x, 
        n ≥ 0 ∧ n ∈ Integers]

Table[intx2f[n], {n, 10}]

define tensors

This (Subscript[x, a])or (x^a)denotes a component of the random vector x.

DefineTensor[tx, "x", {{1}, 1}]

This (Subscript[b, a])or (b^a)denotes a component of the mean vector b.

DefineTensor[sb, "b", {{1}, 1}]

The following a0, a1, a2, a3, a4, a5, and a6 are used for coefficients in a series expansion with respect to x.

DefineTensor[ta0, "a0", {{}, 1}]

DefineTensor[ta1, "a1", {{1}, 1}]

DefineTensor[ta2, "a2", {{2, 1}, 1}]

DefineTensor[ta3, "a3", {{1, 2, 3}, 1}]

SetSymmetric[ta3[la, lb, lc]]

DefineTensor[ta4, "a4", {{1, 2, 3, 4}, 1}]

SetSymmetric[ta4[la, lb, lc, ld]]

DefineTensor[ta6, "a6", {{1, 2, 3, 4, 5, 6}, 1}]

SetSymmetric[ta6[la, lb, lc, ld, le, lf]]

the central moments for the multivariate case

Here we assume b = 0, and define the second central moment (Subscript[α, ab]) = E((Subscript[x, a]) (Subscript[x, b])), the fourth central moment (Subscript[α, abcd]) = 
        E((Subscript[x, a]) (Subscript[x, b]) (Subscript[x, c]) (Subscript[x, d])), etc. The central moments of odd degrees are zero. Considering the symmetry of the normal distribution, these central moments of even degrees are expressed by the following partition indicators.

kd2[la_, lb_] = Kdelta[la, lb];

kd4[la_, lb_, lc_, ld_] = 
        3  Symmetrize[
            Kdelta[la, lb] Kdelta[lc, ld], {la, lb, lc, ld}];

kd6[la_, lb_, lc_, ld_, le_, lf_] =  
        15  Symmetrize[
            Kdelta[la, lb] Kdelta[lc, ld] Kdelta[le, lf], {la, lb, lc, ld, 
              le, lf}];

Then,  (Subscript[α, ab]) = kd2[a, b], (Subscript[α, abcd]) = kd4[a, b, c, d], (Subscript[α, abcdef]) = kd6[a, b, c, d, e, f], etc.

RuleUnique[ruleintx1, tx[ui_], 0]

RuleUnique[ruleintx2, tx[ui_] tx[uj_], kd2[ui, uj]]

RuleUnique[ruleintx3, tx[ui_] tx[uj_] tx[uk_], 0]

RuleUnique[ruleintx4, tx[ui_] tx[uj_] tx[uk_] tx[ul_], 
      kd4[ui, uj, uk, ul]]

RuleUnique[ruleintx5, tx[ui_] tx[uj_] tx[uk_] tx[ul_] tx[um_], 
      0]

RuleUnique[ruleintx6, tx[ui_] tx[uj_] tx[uk_] tx[ul_] tx[um_] tx[un_], 
      kd6[ui, uj, uk, ul, um, un]]

ruleintx = {ruleintx6, ruleintx5, ruleintx4, ruleintx3, ruleintx2, 
          ruleintx1};

Check if this rule works for one-dimension case; \(TraditionalForm\`α\_11, \ α\_1111, \ α\_111111\)are given by

{kd2[1, 1], kd4[1, 1, 1, 1], kd6[1, 1, 1, 1, 1, 1]}

They should be equal to E((x^2)),  E((x^4)),  E((x^6)).

{intx2f[1], intx2f[2], intx2f[3]}

(Subscript[α, 1122]) = E(((Subscript[x, 1])^2)) E(
            ((Subscript[x, 2])^2)) = ((E((x^2)))^2)

{kd4[1, 1, 2, 2], intx2f[1]^2}

(Subscript[α, 112233]) = E(((Subscript[x, 1])^2)) E(((Subscript[x, 2])^2)) E(
            ((Subscript[x, 3])^2)) = ((E((x^2)))^3)

{kd6[1, 1, 2, 2, 3, 3], intx2f[1]^3}

(Subscript[α, 112222]) = E(((Subscript[x, 1])^2)) E(
            ((Subscript[x, 2])^4)) = E((x^2)) E((x^4))

{kd6[1, 1, 2, 2, 2, 2], intx2f[1] intx2f[2]}

In the below, we rather use superscript to indicate the components of x. For example, we use x = ((x^a)), instead of x=((Subscript[x, a])). (Subscript[a2, ij]) (x^i) (x^j) denotes a quadratic form with symmetric matrix (Subscript[a2, ij]).

ta2[li, lj] tx[ui] tx[uj]

The expectation E((Subscript[a2, ij]) (x^i) (x^j)) is

ApplyRules[%, ruleintx]

AbsorbKdelta[%]

Tsimplify[%]

Similarly, for the symmetric 4-form (Subscript[a4, ijkl]) (x^i) (x^j) (x^k) (x^l)

ta4[li, lj, lk, ll] tx[ui] tx[uj] tx[uk] tx[ul]

ApplyRules[%, ruleintx]

AbsorbKdelta[%]

Similarly, for the symmetric 6-form \(TraditionalForm\`\(a6\_ijklmn\) \(x\^i\) \(x\^j\) \(x\^k\) \(x\^l\) \(x\^m\) x\^n\)

ta6[li, lj, lk, ll, lm, ln] tx[ui] tx[uj] tx[uk] tx[ul] tx[um] 
      tx[un]

ApplyRules[%, ruleintx]

AbsorbKdelta[%]

the expectation of the exponential of polynomial functions

Here we would like to calculate the log of the expectation of  the exponential of poly(x) = 
        a0 + (Subscript[a1, i]) (x^i) + (Subscript[a2, ij]) (x^i) 
            (x^j) + (Subscript[a3, ijk]) (x^i) (x^j) 
            (x^k) + (Subscript[a4, ijkl]) (x^i) (x^j) (x^k) (x^l), where a1, a2, and a3 are of order O((n^(-1/2))), and a4 is O((n^(-1))). Note that (x^i) here indicates the i-th element of x=((x^1),  ... , (x^dim))  instead of the i-th power of x.  logeexppoly = log E {exp(poly(x))} is first obtained for b = 0, and next for general b ≠ 0up to O((n^(-1)))terms.

central case (b=0)

First poly(x)-a0 is set in foo1. The constant term a0 can be removed from poly in the following argument, but only be added in logexppoly later.

foo1 = 
      o ta1[li] tx[ui] + o ta2[li, lj] tx[ui] tx[uj] + 
        o ta3[li, lj, lk] tx[ui] tx[uj] tx[uk] + (o^(2 )) 
          ta4[li, lj, lk, ll] tx[ui] tx[uj] tx[uk] tx[ul]

Make rule "rulepoly" to substitute "poly" for foo1.

RuleUnique[rulepoly, poly, foo1]

Considering

Series[Exp[x], {x, 0, 2}]

we can calculate exp(poly(x)) as follows, and stored in foo4 up to O((n^(-1))) terms.

foo3 = 
        ApplyRules[1 + poly + (1/2) (poly^2), rulepoly];

foo4 = 
      Sum[Tsimplify[Coefficient[foo3, o, i]] (o^i), {i, 0, 2}]

Then, we take the expectation of foo4, calculating E {exp(poly(x))}, and stored in foo5 below.

ApplyRules[foo4, ruleintx]

AbsorbKdelta[%]

foo5 = Tsimplify[%]

Considering

Series[Log[1 + x], {x, 0, 2}]

we take the log of foo5, and stored in foo8 up to O((n^(-1))) terms below.

RuleUnique[foo6rule, foo6, foo5 - 1]

foo7 = ApplyRules[foo6 - (foo6^2)/2, foo6rule];

foo8 = 
      Sum[CanAll[Tsimplify[Coefficient[foo7, o, i]]] (o^i), {i, 0, 
          2}]

This is logeexppoly for b=0.

logeexppolyb0 = ta0 + foo8

The following expression may be easier to read for us, but violating the summation convention rule of subscripts.

logeexppolyb0 /. {u1 -> l1, u2 -> l2, 
        u3 -> l3}

noncentral case (b≠0)

We assume the mean of x is zero, but replacing (Subscript[x, i])with (x^i) + (b^i) to calculate the case for b ≠ 0. foo11 is the same as poly(x)-a0, but with this substitution.

foo11 = foo1 /. {tx[ui_] -> tx[ui] + sb[ui]}

Obtain the canonical expression of the tensor usage, and redefine "rulepoly".

foo12 = CanAll[Expand[foo11]]

RuleUnique[rulepoly, poly, foo12]

we can calculate exp(poly(x)) as follows, and stored in foo14 up to O((n^(-1))) terms.

foo13 = 
        ApplyRules[1 + poly + (1/2) (poly^2), rulepoly];

foo14 = 
      Sum[Tsimplify[Coefficient[foo13, o, i]] (o^i), {i, 0, 2}]

Then, we take the expectation of foo14, calculating E {exp(poly(x))}, and stored in foo15 below.

ApplyRules[foo14, ruleintx]

AbsorbKdelta[%]

foo15 = Tsimplify[%]

we take the log of foo15, and stored in foo18 up to O((n^(-1))) terms below.

RuleUnique[foo16rule, foo16, foo15 - 1]

foo17 = ApplyRules[foo16 - (foo16^2)/2, foo16rule];

foo18 = 
      Sum[CanAll[Tsimplify[Coefficient[foo17, o, i]]] (o^i), {i, 0, 
          2}]

This is logeexppoly for b≠0.

logeexppoly = ta0 + foo18

InputForm[logeexppoly]

The following expression may be easier to read for us, but violating the summation convention rule of subscripts.

Collect[
      logeexppoly, {o, sb[l1] sb[l2] sb[l3] sb[l4], sb[l1] sb[l2] sb[l3], 
        sb[l1] sb[l2], sb[l1]}]

% /. {u1 -> l1, u2 -> l2, u3 -> l3, u4 -> l4, 
        u5 -> l5}

Exponential family

The density function of the exponential family is first specified by using the natural parameter vector. This standard form is transformed into our canonical expression of the density with respect to the expectation parameter vector. Since the exponential family is not uniquely expressed up to the affine transformation, we assume without loss of generality that origin of the expectation parameter coincides with that of the natural parameter, and the covariance matrix of the random variable is identity at the origin. We consider only the continuous random variable throughout. The canonical form will be stored in "logdensityy".

the standard form

Let f(y; θ) = 
        exp((θ^a) (Subscript[y, a]) - h(y) - ψ(θ)) denote the density function of random variable y = ((Subscript[y, 1]),  ... , (Subscript[y, dim])) with the natural parameter vector \(TraditionalForm\`θ = \((θ\^1,  ... , θ\^dim)\)\). The log of the density function is specified by logdensity=log f(y;θ) using the cumulant function ψ(θ) and the measure function h(y).

simplification functions

Some simplification functions are defined here for later use

tsimp[exp_] := CanAll[AbsorbKdelta[CanAll[exp]]]

geto2[exp_] := 
      Sum[Simplify[Coefficient[exp, o, i]] (o^i), {i, -1, 2}]

tgeto2[exp_] := 
      Sum[tsimp[Coefficient[exp, o, i]] (o^i), {i, -1, 2}]

Define differential operator (for type-a index)

difa[exp_, ru_, 
        ala_] := (exp  -  (exp /. {ru[al_] -> 0})) /. {ru[
                al1_] ru[al2_] ru[al3_] ru[al4_] -> 
            ru[al1] ru[al2] ru[al3] Kdelta[al4, ala] + 
              ru[al1] ru[al2] ru[al4] Kdelta[al3, ala] + 
              ru[al1] ru[al4] ru[al3] Kdelta[al2, ala] + 
              ru[al4] ru[al2] ru[al3] Kdelta[al1, ala], 
ru[
                al1_] ru[al2_] ru[al3_] -> 
            ru[al1] ru[al2] Kdelta[al3, ala] + 
              ru[al1] ru[al3] Kdelta[al2, ala] + 
              ru[al2] ru[al3] Kdelta[al1, ala], 
          ru[al1_] ru[al2_] ->  
            ru[al1] Kdelta[al2, ala] + ru[al2] Kdelta[al1, ala], 
          ru[al1_] -> Kdelta[al1, ala]}

define tensors

\(TraditionalForm\`θ = \((θ\^1,  ... , θ\^dim)\)\) is the natural parameter vector

DefineTensor[st, "θ", {{1}, 1}]

η = ((Subscript[η, 1]),  ... , (Subscript[η, dim])) is the expectation parameter vector

DefineTensor[se, "η", {{1}, 1}]

\(TraditionalForm\`ϕ3\^abc = \(\(∂\^3\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c\)\)\( | \_0\)\)\) is the third derivative of the potential function ϕ(η) at the origin η=0.

DefineTensor[tp3, "ϕ3", {{1, 2, 3}, 1}]

SetSymmetric[tp3[la, lb, lc]]

\(TraditionalForm\`ϕ4\^abcd = \(\(∂\^4\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c ∂η\_d\)\)\( | \_0\)\)\) is the fourth derivative of the potential function ϕ(η) at the origin η=0.

DefineTensor[tp4, "ϕ4", {{1, 2, 3, 4}, 1}]

SetSymmetric[tp4[la, lb, lc, ld]]

\(TraditionalForm\`ψ3\_abc = \(\(∂\^3 ψ\/\(∂θ\^a ∂θ\^b ∂θ\^c\)\)\( | \_0\)\)\) is the third derivative of the cumulant function ψ(θ) at the origin θ=0.

DefineTensor[tq3, "ψ3", {{1, 2, 3}, 1}]

SetSymmetric[tq3[la, lb, lc]]

\(TraditionalForm\`ψ4\_abcd = \(\(∂\^4 ψ\/\(∂θ\^a ∂θ\^b ∂θ\^c ∂θ\^d\)\)\( | \_0\)\)\) is the fourth derivative of the cumulant function ψ(θ) at the origin θ=0.

DefineTensor[tq4, "ψ4", {{1, 2, 3, 4}, 1}]

SetSymmetric[tq4[la, lb, lc, ld]]

y = ((Subscript[y, 1]),  ... , (Subscript[y, dim])) is the random variable

DefineTensor[ry, "y", {{1}, 1}]

\(TraditionalForm\`h1\^a = \(\(∂h\/∂y\_a\)\( | \_0\)\)\) is the first derivative of h(y) at the origin y = 0.

DefineTensor[th1, "h1", {{1}, 1}]

\(TraditionalForm\`h2\^ab = \(\(∂\^2 
                h\/\(∂y\_a ∂y\_b\)\)\( | \_0\)\)\) is the second derivative of h(y) at the origin y = 0.

DefineTensor[th2, "h2", {{2, 1}, 1}]

(Subscript[ih, 2  ab]) = ((((h2^ab)))^(-1)) is the inverse matrix of (h2^ab).

DefineTensor[tih2, "ih2", {{2, 1}, 1}]

\(TraditionalForm\`h3\^abc = \(\(∂\^3 
                h\/\(∂y\_a ∂y\_b ∂y\_c\)\)\( | \_0\)\)\) is the third derivative of h(y) at the origin y = 0.

DefineTensor[th3, "h3", {{1, 2, 3}, 1}]

SetSymmetric[th3[la, lb, lc]]

\(TraditionalForm\`h4\^abcd = \(\(∂\^4 
                h\/\(∂y\_a ∂y\_b ∂y\_c ∂y\_d\)\)\( | \_0\)\)\) is the fourth derivative of h(y) at the origin y = 0.

DefineTensor[th4, "h4", {{1, 2, 3, 4}, 1}]

SetSymmetric[th4[la, lb, lc, ld]]

the log of the density function

logdensity=log f(y;θ) is specified here.

logdensity = st[ua] ry[la] - psi[theta] - h[y]

Since ∫f(y; θ) dy = 
        1, the cumulant function is defined formally by ψ(θ) = 
        log ∫exp((θ^a) (Subscript[y, a]) - 
                  h(y)) dy. The expectation parameter vector is defined by η = E(
            y;  θ) = ∫y f(
                y; θ) dy. The potential function ϕ(η) is defined by ϕ(η) = 
        (max_θ){(θ^a) (Subscript[η, a]) - ψ(θ)}. The two parametrizations are related to each other by \(TraditionalForm\`η\_a = ∂ψ\/∂θ\^a\)and \(TraditionalForm\`θ\^a = ∂ϕ\/∂η\_a\).

Without losing generality, we assume \(TraditionalForm\`\(\(∂ϕ\/∂η\_a\)\( | \_0\)\) = 0\) and \(TraditionalForm\`\(\(∂\^2\ ϕ\/\(∂η\_a ∂η\_b\)\)\( | \_0\)\) = δ\^ab\).

the expression of ψ(θ) in terms of η

In this subsection,  we derive the expression of ψ(θ) using η and  the ϕ derivatives. The result will be stored in "psieta".

derivation

Define "ruleeta1" for the Taylor series of (Subscript[η, a]) is \(TraditionalForm\`η\_a = \(\(∂\ ψ\/∂\ θ\^a\)\(=\)\)\)\(TraditionalForm\`∂\ ψ\/∂\ θ\^a\( | \_0\)\(+\(∂\^2\ ψ\/\(∂\ θ\^a ∂\ θ\^b\)\)\)\( | \_0\)θ\^b + \(1\/2\) ∂\^3\ ψ\/\(∂\ θ\^a ∂\ θ\^b ∂\ θ\^c\)\( | \_0\)\(θ\^b\) θ\^c + \(1\/6\) ∂\^4\ ψ\/\(∂\ θ\^a ∂\ θ\^b ∂\ θ\^c ∂\ θ\^d\)\( | \_0\)\(θ\^b\) \(θ\^c\) θ\^d\), where ψ(Subscript[3, abc]) = O((n^(-1/2))) and ψ(Subscript[4, abcd]) = O((n^(-1))).

RuleUnique[ruleeta1, se[la_], 
      Kdelta[la, lb] st[ub] + (1/2) o tq3[la, lb, lc] st[ub] 
          st[uc] + (1/6) (o^2) tq4[la, lb, lc, ld] st[ub] st[uc] 
          st[ud]]

ApplyRules[se[la], ruleeta1]

Define "ruletheta1" for the Taylor series of \(TraditionalForm\`θ\^a = \(\(∂\ ϕ\/∂\ η\_a\)\(=\)\)\)\(TraditionalForm\`∂\ ϕ\/∂η\_a\( | \_0\)\(+\(∂\^2\ ϕ\/\(∂η\_a ∂η\_b\)\)\)\( | \_0\)η\_b + \(1\/2\) ∂\^3\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c\)\( | \_0\)\(η\_b\) η\_c + \(1\/6\) ∂\^3\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c ∂η\_d\)\( | \_0\)\(η\_b\) \(η\_c\) η\_d\), where ϕ(3^abc) = O((n^(-1/2))) and ϕ(4^abcd) = O((n^(-1))).

RuleUnique[ruletheta1, 
      st[ua_], 
Kdelta[ua, ub] se[lb] + (1/2) 
          o tp3[ua, ub, uc] se[lb] se[lc] + (1/6) (o^2) 
          tp4[ua, ub, uc, ud] se[lb] se[lc] se[ld]]

foo21 = ApplyRules[st[ua], ruletheta1]

Calculate \(TraditionalForm\`∂\^2\ ϕ\/\(∂η\_a ∂η\_b\) = \(∂\ θ\^a\/∂η\_b = \(\(δ\^ab\)\(+\)\)\)\)\(TraditionalForm\`∂\^3\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c\)\( | \_0\)η\_c + \(1\/2\) ∂\^4\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c ∂η\_d\)\( | \_0\)\(η\_c\) η\_d\) by taking the partial differentiation of (θ^a) with respect to (Subscript[η, b]). Define "rulephi2" for \(TraditionalForm\`∂\^2\ ϕ\/\(∂η\_a ∂η\_b\)\).

foo22 = tsimp[difa[foo21, se, ub]]

DefineTensor[dp2, {{2, 1}, 1}]

RuleUnique[rulephi2, dp2[ua_, ub_], foo22]

Apply (((I + A))^(-1)) = 
        I - A + (A^2) + O((n^(-3/2))) for A = O((n^(-1/2))) to \(TraditionalForm\`∂\^2\ ϕ\/\(∂η\_a ∂η\_b\)\)

foo23 = 
      tgeto2[ApplyRules[
          Kdelta[ua, 
              ub] - (dp2[ua, ub] - 
                Kdelta[ua, ub]) + (dp2[ua, uc] - 
                  Kdelta[ua, uc]) (dp2[lc, ub] - Kdelta[lc, ub]), 
          rulephi2]]

Substitute ruleeta1 for (Subscript[η, a]) to get \(∂\^2 ψ\/\(∂θ\^a ∂θ\^b\)\)=\(\((∂\^2\ ϕ\/\(∂η\_a ∂η\_b\))\)\^\(-1\)\) using ϕ derivatives as well as ψ derivatives.

foo24 = tgeto2[ApplyRules[foo23, ruleeta1]]

The above expression is compared with \(TraditionalForm\`\(\(∂\^2 ψ\/\(∂θ\^a ∂θ\^b\)\)\(=\)\)\)(Subscript[δ, ab])+\(TraditionalForm\`∂\^3 ψ\/\(∂θ\^a ∂θ\^b ∂θ\^c\)\( | \_0\)θ\^c + ∂\^4 ψ\/\(∂θ\^a ∂θ\^b ∂θ\^c ∂θ\^d\)\( | \_0\)\(θ\^c\) θ\^d\), and the coefficients are obtained below.

foo25 = 
      Collect[CoefficientList[foo24 /. {st[l1_] -> x}, x], o, 
        tsimp[CanAll[# /. {u1 -> uc, u2 -> ud, u3 -> ue, 
                  u4 -> uf, l1 -> lc, l2 -> ld, l3 -> le, 
                  l4 -> lf}]] &]

foo26 = 
      Simplify[{foo25[[2]]/o, foo25[[3]] (2/(o^2))}]

RuleUnique[ruletq3, tq3[ua_, ub_, uc_], foo26[[1]]]

RuleUnique[ruletq4, tq4[ua_, ub_, uc_, ud_], 
      ApplyRules[foo26[[2]], ruletq3]]

We have obtained \(TraditionalForm\`\(\(∂\^3 ψ\/\(∂θ\^a ∂θ\^b ∂θ\^c\)\)\( | \_0\)\)\), and \(TraditionalForm\`\(\(∂\^4 ψ\/\(∂θ\^a ∂θ\^b ∂θ\^c ∂θ\^d\)\)\( | \_0\)\)\)  as follows.

ApplyRules[tq3[ua, ub, uc], ruletq3]

ApplyRules[tq4[ua, ub, uc, ud], ruletq4]

Let us write down the Taylor series ψ(θ) = ψ(0)+\(TraditionalForm\`∂\ ψ\/∂\ θ\^a\( | \_0\)θ\^a + \(1\/2\) ∂\^2\ ψ\/\(∂\ θ\^a ∂\ θ\^b\)\( | \_0\)\(θ\^a\) θ\^b + \(1\/6\) ∂\^3\ ψ\/\(∂\ θ\^a ∂\ θ\^b ∂\ θ\^c\)\( | \_0\)\(θ\^a\) \(θ\^b\) θ\^c + \(1\/24\) ∂\^4\ ψ\/\(∂\ θ\^a ∂\ θ\^b ∂\ θ\^c ∂\ θ\^d\)\( | \_0\)\(θ\^a\) \(θ\^b\) \(θ\^c\) θ\^d\)

psitheta = 
      psi[0] + (1/2) Kdelta[la, lb] st[ua] st[ub] + 
        (1/6) o tq3[la, lb, lc] st[ua] st[ub] st[uc] + (1/24) 
          (o^2) tq4[la, lb, lc, ld] st[ua] st[ub] st[uc] st[ud]

result

Now express ψ(θ) in terms of η and ϕ derivatives. This gives eq.(3.8) of SH02.

psieta = 
      CanAll[ApplyRules[
          CanAll[tgeto2[ApplyRules[psitheta, ruletheta1]]], {ruletq3, 
            ruletq4}]]

psieta // InputForm

the expression of h(y) in terms of ϕ derivatives

In this subsection, we derive the expression of h(y) using ϕ derivatives. The result will be stored in "hinyphi".

derivation

Consider Taylor series h(y) = 
        h(0) + (h1^a) (Subscript[y, a]) + (1/2) (h2^ab) (Subscript[y, a]) 
            (Subscript[y, b]) + (1/6) (h3^abc) (Subscript[y, a]) (Subscript[y, b]) 
            (Subscript[y, c]) + (1/24) (h4^abcd) (Subscript[y, a]) (Subscript[y, b]) (Subscript[y, c]) 
            (Subscript[y, d]), where (h1^a) = O((n^(1/2))),  (h2^ab) = O(1),  
      (h3^abc) = O((n^(-1/2))),  
      (h4^abcd) = O((n^(-1))).

hiny = 
      h[0] + (o^(-1)) th1[ua] ry[la] + (1/2) th2[ua, ub] ry[la] 
          ry[lb] + (1/6) o th3[ua, ub, uc] ry[la] ry[lb] 
          ry[lc] + (1/24) (o^2) th4[ua, ub, uc, ud] ry[la] ry[lb] 
          ry[lc] ry[ld]

Define foo31 and foo32 below, which will satisfy foo31 + foo32 = (θ^a) (Subscript[y, a]) - h(y)  +  
          O((n^(-1/2)))as shown later.

foo31 = ((-1)/2) 
        th2[ua, ub] (ry[la] - 
            tih2[la, lc] (st[uc] - (o^(-1)) th1[uc])) (ry[lb] - 
            tih2[lb, ld] (st[ud] - (o^(-1)) th1[ud]))

foo32 = -h[0] +  (1/2) 
          tih2[la, 
            lb] (st[ua] - (o^(-1)) th1[ua]) (st[
                ub] - (o^(-1)) th1[ub])

Define a rule to make ((Subscript[h2, a])^b) (ih2^ac) = (δ^bc), thus ih2 is the inverse matrix of h2.

RuleUnique[absorbth2, th2[la_, ub_] tih2[uc_, ua_], 
      Kdelta[ub, uc]]

Apply this rule to foo31+foo32 to get foo33.

foo33 = 
      CanAll[AbsorbKdelta[
          ApplyRules[CanAll[foo31 + foo32], absorbth2]]]

Confirming foo33=foo31 + foo32 = (θ^a) (Subscript[y, a]) - 
          h(y)  +  O((n^(-1/2))).

CanAll[foo33 - (ry[la] st[ua] - hiny)]

Since (((2  Pi))^(((-dim)/2))) ((det(
              (h2^ab)))^((1/2))) exp(foo31) is the multivariate normal density function with mean specified above and (Subscript[ih2, ab]) covariance, ∫exp((θ^a) (Subscript[y, a]) - h(y)  - foo32 +  
                  O((n^(-1/2))))dy=\(TraditionalForm\`∫\(exp(
              foo31)\) dy = \(\((2  π)\)\^\(dim\/2\)\) \(det(h2\^ab)\)\^\(-\(1\/2\)\)\). Then, ψ(θ)=log∫exp((θ^a) (Subscript[y, a]) - 
                h(y)) dy = (dim/2) log(
              2  π) - (1/2) det((h2^ab)) + 
          foo32 +  O((n^(-1/2))).

foo34 = (dim/2) Log[2  Pi] - (1/2) deth2 + 
        CanAll[Expand[foo32]]

Considering \(TraditionalForm\`\(\(∂\ ψ\/∂\ θ\^a\)\( | \_0\)\) = 
        0, \(\(∂\^2\ ψ\/\(∂\ θ\^a ∂\ θ\^b\)\)\( | \_0\)\) = δ\_ab\), we find that (ih2^ab) = (δ^ab) + O((n^(-1/2))) and (Subscript[h1, a]) = O((n^(-1/2)))from the above expression of ψ(θ) = foo34 + O((n^(-1/2))).  Let us write (h2^ab) = (δ^ab) + (ah2^ab)with (ah2^ab) = O((n^(-1/2))).

DefineTensor[tah2, "ah2", {{2, 1}, 1}]

ruleth2ah2 = 
        th2[ua_, ub_] -> Kdelta[ua, ub] + o tah2[ua, ub];

Rewrite the Taylor series of h(y) by noting the first derivative h1 is in fact O((n^(-1/2))).

hiny = 
      h[0] + o th1[ua] ry[la] + (1/2) th2[ua, ub] ry[la] 
            ry[lb] + (1/6) o th3[ua, ub, uc] ry[la] ry[lb] 
            ry[lc] + (1/24) (o^2) th4[ua, ub, uc, ud] ry[la] ry[lb] 
            ry[lc] ry[ld] /. ruleth2ah2

In the following, exp((θ^a) (Subscript[y, a]) - h(y)) = exp(
            foo35) = exp(foo36) exp(
              foo37 + O((n^(-3/2)))), where ∫exp(foo31) dy = 1.

foo35 = ry[la] st[ua] - hiny

foo36 = ((-dim)/2) 
          Log[2  Pi] + ((-1)/2) (ry[la] -  st[la]) (ry[ua] -  
              st[ua])

foo37 = CanAll[AbsorbKdelta[CanAll[foo35 - foo36]]]

Noting ψ(θ) = log ∫exp(
                  foo35) dy = 
          log E(exp(foo37 + O((n^(-3/2))))), where the expectation is taken for the multivariate normal with mean ((θ^1),  ... , (θ^dim)) and covariance identity. We apply "logeexppoly" to foo37. First, we get foo38={a0,a1,a2,a3,a4} for the coefficients of a0 + (a1^i) (Subscript[x, i]) + (a2^ij) (Subscript[x, i]) 
          (Subscript[x, j]) + (a3^ijk) (Subscript[x, i]) (Subscript[x, j]) 
          (Subscript[x, k]) + (a4^ijkl) (Subscript[x, i]) (Subscript[x, j]) (Subscript[x, k]) (Subscript[x, l]), and ψ(θ) = foo39below.

foo38 = 
      CoefficientList[foo37 /. ry[la_] -> x, x]/{1, o, o, o, 
          (o^2)}

foo39 = 
      CanAll[logeexppoly /. {sb -> st, ta0 -> foo38[[1]], 
            ta1[u1_] -> foo38[[2]], 
            ta2[u1_, u2_] -> foo38[[3]], 
            ta3[u1_, u2_, u3_] -> foo38[[4]], 
            ta4[u1_, u2_, u3_, u4_] -> foo38[[5]]}]

Now substitute \(TraditionalForm\`θ\^a = \(\(∂\ ϕ\/∂\ η\_a\)\(=\)\)\)\(TraditionalForm\`\(δ\^ab\) η\_b + \(1\/2\) \(ϕ\^abc\) \(η\_b\) η\_c + \(1\/6\) \(ϕ\^abcd\) \(η\_b\) \(η\_c\) η\_d\) for (θ^a) in foo39 to get ψ(θ)=foo40 in terms of η below. The coefficients for the polynomial of η is stored in foo41.

foo40 = CanAll[tgeto2[ApplyRules[foo39, ruletheta1]]];

foo41 = 
      Collect[CoefficientList[foo40 /. {se[l1_] -> x}, x], o, 
        tsimp[# /. {u1 -> ua, u2 -> ub, u3 -> uc, 
                u4 -> ud, l1 -> la, l2 -> lb, l3 -> lc, 
                l4 -> ld}] &]

Since foo40=ψ(θ)=psieta, foo41 must be equal to foo42 below.

foo42 = 
      CoefficientList[psieta /. se[la_] -> x, x] /. {u1 -> ua, 
          u2 -> ub, u3 -> uc, u4 -> ud}

At first, we compare the coefficient for (η^a) (η^b).

foo41[[3]]

foo42[[3]]

Thus, ah2 is in fact O((n^(-1))) instead of O((n^(-1/2))).

hiny = hiny /. { tah2[ua_, ub_] -> o tah2[ua, ub]}

We rewrite foo43=foo42-foo41.

foo43 = 
      Map[tgeto2, (foo42 - foo41) /.  {tah2[u1_, u2_] -> 
              o tah2[u1, u2]}]

We solve these five equations == 0. At first, foo43[[4]]==0 gives (h3^abc) = ϕ(3^abc).

foo44 = Solve[foo43[[4]] == 0, th3[ua, ub, uc]]

RuleUnique[rule44, th3[ua_, ub_, uc_], foo44[[1, 1, 2]]]

foo45 = 
      Solve[ApplyRules[foo43[[2]], rule44] == 0, 
        th1[ua]]

RuleUnique[rule45, th1[ua_], foo45[[1, 1, 2]]]

foo46 = 
      Solve[ApplyRules[foo43[[5]], {rule44, rule45}] == 0, 
        th4[ua, ub, uc, ud]]

RuleUnique[rule46, th4[ua_, ub_, uc_, ud_], 
      foo46[[1, 1, 2]]]

foo47 = 
      Simplify[Solve[
          ApplyRules[foo43[[3]], {rule44, rule45, rule46}] == 0, 
          tah2[ua, ub]]]

RuleUnique[rule47, tah2[ua_, ub_], foo47[[1, 1, 2]]]

foo48 = 
      Simplify[Solve[
          ApplyRules[
              foo43[[1]], {rule44, rule45, rule46, rule47}] == 0, 
          h[0]]]

RuleUnique[rule48, h[0], foo48[[1, 1, 2]]]

result

Now, it is time to express h(y) in terms of the ϕ derivatives.

hinyphi = 
      CanAll[AbsorbKdelta[
          ApplyRules[
            hiny, {rule44, rule45, rule46, rule47, rule48}]]]

hinyphi // InputForm

The coefficients of h(y) as a polynomial of y are as follows.

foo49 = 
      Collect[CoefficientList[hinyphi /. ry[l1_] -> x, 
            x] /. {u1 -> ua, u2 -> ub, u3 -> uc, 
            u4 -> ud, l1 -> la, l2 -> lb, l3 -> lc, 
            l4 -> ld}, o, tsimp[#] &]

the canonical form

Now we get the canonical form of log f(y; η) in "logdensityy".

the summary of the previous sections

The standard form of the density function is specified by logdensity=log f(y;θ).

logdensity

By applying "ruletheta1" to logdensity, θ can be expressed in terms of η and ϕ derivatives.

foo51 = ApplyRules[logdensity, ruletheta1]

The cumulant function ψ(θ) can be expressed in terms of η and ϕ derivatives as shown in psieta.

psieta

The measure function h(y) can be expressed in terms of y and ϕ derivatives as shown in hinyphi.

hinyphi

result

By substituting psieta and hinyphi for psi[theta] and h[y] respectively in foo51=log f(y;θ), we obtain the canonical form of log f(y;η) as follows.

foo52 = 
      Collect[foo51 /. {psi[theta] -> psieta, h[y] -> hinyphi}, o, 
        tsimp]

This is the canonical form  of log f(y;η).

logdensityy =  Expand[foo52];

logdensityy // InputForm

We also show the metric \(TraditionalForm\`∂\^2\ ϕ\/\(∂η\_a ∂η\_b\)\) by using rulephi2, and name it "phi2eta" for later use.

phi2eta = ApplyRules[dp2[ua, ub], rulephi2]

InputForm[phi2eta]

Tube-Coordinates and (Subscript[z, c])-formula

In this part, we first give an expression of the smooth surface in (R^dim) which specifies the boundary of the region of interest in the η-space. The tube-coordinate system is then defined as a pair of the coordinate system on the surface and the coordinate along the normal direction. The signed distance from the boundary is slightly modified for generalization, and named as  a modified signed distance characterized  by a coefficient vector c. The (Subscript[z, c])-formula is derived as an expression of  the distribution function of the modified signed distance, which is obtained up to O((n^(-1))) terms ignoring the error of O((n^(-3/2))).

Startup

This section initializes the Mathematica session.

packages

<<  Statistics`ContinuousDistributions`

<< MathTens.m (Windows)

error messages

Off[General::spell1]

Off[General::spell]

distribution functions

gammadist[x_, m_, α_] := 
      PDF[GammaDistribution[m, α], x]
Gammadist[x_, m_, α_] := 
      CDF[GammaDistribution[m, α], x]
f[x_] := PDF[NormalDistribution[0, 1], x]
F[x_] := CDF[NormalDistribution[0, 1], x]
Q[x_] := Quantile[NormalDistribution[0, 1], x]
Chidist[x_, {di_, nc_}] := 
      CDF[NoncentralChiSquareDistribution[di, nc], x]

Exponential family

This section summarizes the results for the exponential family derived in the previous part.

the expectation of the exponential of a polynomial function of the normal vector

Here we give the log of the expectation of  the exponential of poly(x) = 
        a0 + (Subscript[a1, i]) (x^i) + (Subscript[a2, ij]) (x^i) 
            (x^j) + (Subscript[a3, ijk]) (x^i) (x^j) 
            (x^k) + (Subscript[a4, ijkl]) (x^i) (x^j) (x^k) (x^l), where a1, a2, and a3 are of order O((n^(-1/2))), and a4 is O((n^(-1))). x is a multivariate normal random vector x=((x^1),  ... , (x^dim)) of dim-dimensions with mean b = ((b^1),  ... , (b^dim)), and the identity covariance matrix.  logeexppoly = log E {exp(poly(x))} is obtained up to O((n^(-1)))terms.

define tensors

This (Subscript[b, a])or (b^a)denotes a component of the mean vector b.

DefineTensor[sb, "b", {{1}, 1}]

The following a0, a1, a2, a3, a4, a5, and a6 are used for coefficients in a series expansion with respect to x.

DefineTensor[ta0, "a0", {{}, 1}]

DefineTensor[ta1, "a1", {{1}, 1}]

DefineTensor[ta2, "a2", {{2, 1}, 1}]

DefineTensor[ta3, "a3", {{1, 2, 3}, 1}]

SetSymmetric[ta3[la, lb, lc]]

DefineTensor[ta4, "a4", {{1, 2, 3, 4}, 1}]

SetSymmetric[ta4[la, lb, lc, ld]]

DefineTensor[ta6, "a6", {{1, 2, 3, 4, 5, 6}, 1}]

SetSymmetric[ta6[la, lb, lc, ld, le, lf]]

"logeexppoly"

logeexppoly = 
      ta0 + o*(sb[l1]*ta1[u1] + ta2[l1, u1] + sb[l1]*sb[l2]*ta2[u1, u2] + 
              3*sb[l1]*ta3[l2, u1, u2] + 
              sb[l1]*sb[l2]*sb[l3]*ta3[u1, u2, u3]) + 
        o^2*((ta1[l1]*ta1[u1])/2 + 2*sb[l1]*ta1[l2]*ta2[u1, u2] + 
              ta2[l1, l2]*ta2[u1, u2] + 
              2*sb[l1]*sb[l2]*ta2[l3, u1]*ta2[u2, u3] + 
              3*ta1[l1]*ta3[l2, u1, u2] + 
              6*sb[l1]*ta2[l2, u1]*
                ta3[l3, u2, u3] + (9*ta3[l1, l2, u1]*ta3[l3, u2, u3])/2 + 
              3*sb[l1]*sb[l2]*ta1[l3]*ta3[u1, u2, u3] + 
              6*sb[l1]*ta2[l2, l3]*ta3[u1, u2, u3] + 
              3*ta3[l1, l2, l3]*ta3[u1, u2, u3] + 
              9*sb[l1]*sb[l2]*ta3[l3, l4, u3]*ta3[u1, u2, u4] + 
              6*sb[l1]*sb[l2]*sb[l3]*ta2[l4, u1]*ta3[u2, u3, u4] + 
              9*sb[l1]*sb[l2]*ta3[l3, l4, u1]*
                ta3[u2, u3, 
                  u4] + (9*sb[l1]*sb[l2]*sb[l3]*sb[l4]*ta3[l5, u1, u3]*
                    ta3[u2, u4, u5])/2 + 3*ta4[l1, l2, u1, u2] + 
              6*sb[l1]*sb[l2]*ta4[l3, u1, u2, u3] + 
              sb[l1]*sb[l2]*sb[l3]*sb[l4]*ta4[u1, u2, u3, u4])

the canonical form of the density function

Here we give the canonical form of log f(y;η) for the exponential family of distributions, and store it in logdensityy. The metric phi2eta=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_a ∂\ η\_b\)\)is also given here.

define tensors

\(TraditionalForm\`θ = \((θ\^1,  ... , θ\^dim)\)\) is the natural parameter vector

DefineTensor[st, "θ", {{1}, 1}]

η = ((Subscript[η, 1]),  ... , (Subscript[η, dim])) is the expectation parameter vector

DefineTensor[se, "η", {{1}, 1}]

y = ((Subscript[y, 1]),  ... , (Subscript[y, dim])) is the random variable

DefineTensor[ry, "y", {{1}, 1}]

\(TraditionalForm\`ϕ3\^abc = \(\(∂\^3\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c\)\)\( | \_0\)\)\) is the third derivative of the potential function ϕ(η) at the origin η=0.

DefineTensor[tp3, "ϕ3", {{1, 2, 3}, 1}]

SetSymmetric[tp3[la, lb, lc]]

\(TraditionalForm\`ϕ4\^abcd = \(\(∂\^4\ ϕ\/\(∂η\_a ∂η\_b ∂η\_c ∂η\_d\)\)\( | \_0\)\)\) is the fourth derivative of the potential function ϕ(η) at the origin η=0.

DefineTensor[tp4, "ϕ4", {{1, 2, 3, 4}, 1}]

SetSymmetric[tp4[la, lb, lc, ld]]

"logdensityy"

logdensityy = -(dim*Log[2*Pi])/2 - (ry[l1]*ry[u1])/2 + 
        ry[l1]*se[u1] - (se[l1]*se[u1])/
          2 + (o*ry[l1]*tp3[l2, u1, u2])/
          2 - (o*ry[l1]*ry[l2]*ry[l3]*tp3[u1, u2, u3])/
          6 + (o*ry[l1]*se[l2]*se[l3]*tp3[u1, u2, u3])/
          2 - (o*se[l1]*se[l2]*se[l3]*tp3[u1, u2, u3])/
          3 + (o^2*tp3[l1, l2, l3]*tp3[u1, u2, u3])/
          6 - (o^2*ry[l1]*ry[l2]*tp3[l3, l4, u1]*tp3[u2, u3, u4])/
          4 - (o^2*tp4[l1, l2, u1, u2])/
          8 + (o^2*ry[l1]*ry[l2]*tp4[l3, u1, u2, u3])/
          4 - (o^2*ry[l1]*ry[l2]*ry[l3]*ry[l4]*tp4[u1, u2, u3, u4])/
          24 + (o^2*ry[l1]*se[l2]*se[l3]*se[l4]*tp4[u1, u2, u3, u4])/
          6 - (o^2*se[l1]*se[l2]*se[l3]*se[l4]*tp4[u1, u2, u3, u4])/
          8

"phi2eta"

phi2eta  =  
      Kdelta[ua, ub] + 
        o*se[l1]*
          tp3[u1, ua, ub] + (o^2*se[l1]*se[l2]*tp4[u1, u2, ua, ub])/
          2

Tube-coordinates

First, the expression of the surface is specified in the Taylor series. Then, the tangent vectors and the normal vector are obtained. The tube-coordinates (u,v) are defined and used instead of the η-parametrization. Here u is dim-1 dimensional vector specifying a point on the surface, and v is the signed distance. The density function f(u, v | v0) is obtained from f(y | η), where the parameter value is specified as η = (0,  ... , 0, v0) without losing the generality.

preliminary

Some functions and tensors are defined here.

indices

Here the dimension of the space is denoted by "9" for the index, whereas it is denoted by "dim" for the regular number.  type-a index may run from 1 to 8, although it is not used explicitly, but only assumes type-a index cannot be 9.

AddIndexTypes

simplification functions

type-a index cannot be "9".

RuleUnique[rulekdelta1, Kdelta[a_, 9], 0, IndexaQ[a]]

RuleUnique[rulekdelta2, Kdelta[a_, -9], 0, IndexaQ[a]]

using the above fact.

tsimp[exp_] := 
        CanAll[AbsorbKdelta[
            CanAll[ApplyRules[
                exp, {rulekdelta1, rulekdelta2}]]]];

raising -9 to 9 for a unique expression.

ruletps9 = {tp3[-9, l1_, l2_] -> tp3[9, l1, l2], 
          tp4[-9, l1_, l2_, l3_] -> 
            tp4[9, l1, l2, l3]};

further simplification

tsimpp[exp_] := 
        Tsimplify[
          CanAll[AbsorbKdelta[
              CanAll[ApplyRules[exp /. ruletps9 /. ruletps9 /. 
                      ruletps9 /. ruletps9, {rulekdelta1, 
                    rulekdelta2}]]]]];

Ignore O((n^(-3/2))) terms for scalar

geto2[exp_] := 
      Sum[Simplify[Coefficient[exp, o, i]] (o^i), {i, -1, 2}]

Ignore O((n^(-3/2))) terms for tensor.

tgeto2[exp_] := 
      Sum[tsimp[Coefficient[exp, o, i]] (o^i), {i, -1, 2}]

Series expansion for scalar ignoring O((n^(-3/2))) terms.

gets2[exp_] := geto2[Series[exp, {o, 0, 2}]]

Series expansion for tensor  ignoring O((n^(-3/2))) terms.

tgetrule[tx_] := 
      Module[{coefs, coef0, coef1, coef2, rule0, rule1, rule2}, 
        coefs = CoefficientList[tx, o]; 
        RuleUnique[rule0, coef0, coefs[[1]]]; 
        RuleUnique[rule1, coef1, coefs[[2]]]; 
        RuleUnique[rule2, coef2, 
          coefs[[3]]]; {{coef0, coef1, coef2}, {rule0, rule1, 
            rule2}}]

tgets2[exp_, x_, tx_] := 
      Module[{xcoef, xrule}, {xcoef, xrule} = tgetrule[tx]; 
        ApplyRules[
          gets2[exp /. {x -> 
                  Sum[xcoef[[i]] (o^(i - 1)), {i, 3}]}], 
          xrule]]

tgets2[exp_, x_, tx_, y_, ty_] := 
      Module[{xcoef, xrule, ycoef, yrule}, {xcoef, xrule} = 
          tgetrule[tx]; {ycoef, yrule} = tgetrule[ty]; 
        ApplyRules[
          gets2[exp /. {x -> Sum[xcoef[[i]] (o^(i - 1)), {i, 3}], 
                y -> Sum[ycoef[[i]] (o^(i - 1)), {i, 3}]}], 
          Join[xrule, yrule]]]

Define an operator to separate the regular index into the type-a index and dim.

sepa[foo_, l_] := 
        Module[{aup, alo}, 
UpLoa[{aup}, {alo}]; 
          MakeSumRange[foo, {l, alo, -9}]];

Define differential operator (for type-a index)

difa[exp_, ru_, 
        ala_] := (exp  -  (exp /. {ru[al_] -> 0})) /. {ru[
                al1_] ru[al2_] ru[al3_] ru[al4_] -> 
            ru[al1] ru[al2] ru[al3] Kdelta[al4, ala] + 
              ru[al1] ru[al2] ru[al4] Kdelta[al3, ala] + 
              ru[al1] ru[al4] ru[al3] Kdelta[al2, ala] + 
              ru[al4] ru[al2] ru[al3] Kdelta[al1, ala], 
ru[
                al1_] ru[al2_] ru[al3_] -> 
            ru[al1] ru[al2] Kdelta[al3, ala] + 
              ru[al1] ru[al3] Kdelta[al2, ala] + 
              ru[al2] ru[al3] Kdelta[al1, ala], 
          ru[al1_] ru[al2_] ->  
            ru[al1] Kdelta[al2, ala] + ru[al2] Kdelta[al1, ala], 
          ru[al1_] -> Kdelta[al1, ala]}

define tensors

u = ((Subscript[u, 1]),  ... , (Subscript[u, dim - 1])) is the parametrization of the surface. Type-a index indicate the suffix such as u = ((Subscript[u, b'])), where b' = 1,  ... , dim - 1.

DefineTensor[ru, "u", {{1}, 1}]

(d^(a' b')) = O((n^(-1/2))) is the curvature matrix (second derivative) of the surface at the origin.

DefineTensor[td, "d", {{2, 1}, 1}]

(e^(a' b' c')) = O((n^(-1)))is the third derivative at the origin.

DefineTensor[te, "e", {{1, 2, 3}, 1}]

SetSymmetric[te[la, lb, lc]]

\(TraditionalForm\`B\_b\^a' = ∂\ η\_b\/∂\ u\_\(a'\)\), b=1,...,dim is the tangent vector for a'=1,...,dim-1. Here the regular type index runs 1,...,dim, whereas the type-a index runs 1,...,dim-1.

DefineTensor[tB, "B", {{1, 2}, 1}]

tB[lb, aua]

(ϕ^(a' b'))(u)=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^a'\)(u)\) \(\(B\_q\^b'\)(u)\)\) is the metric in the tangent space.

DefineTensor[tpa, "ϕ", {{2, 1}, 1}]

the coordinates around the smooth surface

The surface is specified by (Subscript[η, b'])(u) = (Subscript[u, b']), b' = 1,  ... , dim - 1 and (Subscript[η, dim])(
          u) = -(d^(a' b')) (Subscript[u, a']) 
            (Subscript[u, b']) - (e^(a' b' c')) (Subscript[u, a']) (Subscript[u, b']) 
            (Subscript[u, c']).  They are stored in "foo1" and "foo2" or corresponding "rule1" and "rule2".  The region of interest is specified by (Subscript[η, dim]) ≤ (Subscript[η, dim])(u). The tangent vectors are given by foo3=\(TraditionalForm\`B\_\(b'\)\^a' = ∂\ η\_\(b'\)\/∂\ u\_\(a'\)\) and  foo4=\(TraditionalForm\`B\_dim\^a' = ∂\ η\_dim\/∂\ u\_\(a'\)\), or corresponding "rule3" and "rule4".  We also obtain phi2bu= (ϕ^(a' b'))(u)=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^a'\)(u)\) \(\(B\_q\^b'\)(u)\)\) , which is the metric in the tangent space. The elements of the normal vector are denoted as ((Subscript[B, a])^dim), a=1,...,dim, which are given in foo15=((Subscript[B, a'])^dim)(u) and foo16= ((Subscript[B, dim])^dim)(u), or in the corresponding "rule15" and "rule16". The reparametrization between η-coordinates and (u,v)-coordinates are specified by (Subscript[η, a])(u, 
          v) = (Subscript[η, a])(u) + ((Subscript[B, a])^dim)(u) v, and given in foo21 = (Subscript[η, a'])(u, v) and foo22 = (Subscript[η, dim])(u, v), or in "rule21", "rule22", and "rule22b".

smooth surface

Define rules for the surface.

RuleUnique[rule1, se[ala_], Kdelta[ala, aub] ru[alb], 
      IndexaQ[ala]]

RuleUnique[rule2, 
      se[-9], -o td[aua, aub] ru[ala] ru[alb] - (o^2) 
          te[aua, aub, auc] ru[ala] ru[alb] ru[alc]]

(Subscript[η, b'])(u) = (Subscript[u, b']), b' = 1,  ... , dim - 1

foo1 = ApplyRules[se[alb], {rule1, rule2}]

(Subscript[η, dim])(
        u) = -(d^(a' b')) (Subscript[u, a']) 
          (Subscript[u, b']) - (e^(a' b' c')) (Subscript[u, a']) (Subscript[u, b']) 
          (Subscript[u, c'])

foo2 = ApplyRules[se[-9], {rule1, rule2}]

tangent vectors

\(TraditionalForm\`B\_\(b'\)\^a' = ∂\ η\_\(b'\)\/∂\ u\_\(a'\)\)

foo3 = tsimp[difa[foo1, ru, aua]]

\(TraditionalForm\`B\_dim\^a' = ∂\ η\_dim\/∂\ u\_\(a'\)\)

foo4 = tsimp[difa[foo2, ru, aua]]

Define rules for the tangent vectors.

RuleUnique[rule3, tB[alb_, aua_], foo3, 
      IndexaQ[alb] ∧ IndexaQ[aua]]

RuleUnique[rule4, tB[-9, aua_], foo4, IndexaQ[aua]]

check if they work

ApplyRules[tB[alb, aua], {rule3, rule4}]

ApplyRules[tB[-9, aua], {rule3, rule4}]

Next, we will calculate (ϕ^(a' b'))(u)=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^a'\)(u)\) \(\(B\_q\^b'\)(u)\)\) below.

Apply the separate operator to phi2eta=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_a ∂\ η\_b\)\), and evaluate it on the surface to get phi2u=\(TraditionalForm\`\(\(\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_a ∂\ η\_b\)\)\( | \_\(η(
              u)\)\)\)\).

phi2eta

tsimp[sepa[sepa[phi2eta, l1], l2]]

phi2u = tgeto2[ApplyRules[%, {rule1, rule2}]]

Using phi2u above, we write phi2bu= (ϕ^(a' b'))(u)=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^a'\)(u)\) \(\(B\_q\^b'\)(u)\)\) as follows. First, the summation range of phi2bu is separated into type-a and dim. Then,  ((Subscript[B, p])^a)'(u) are substituted by their expressions.

phi2u tB[la, aua] tB[lb, aub]

CanAll[sepa[sepa[%, la], lb]]

phi2bu = tgeto2[ApplyRules[%, {rule3, rule4}]]

We may symmetrize the coefficients of phi2bu.

phi2bucoef = 
      Simplify[CoefficientList[phi2bu /. {ru[al1_] -> x}, x]/{1, o, 
              (o^2)}] /. {au1 -> auc, au2 -> aud}

phi2bucoef[[3]] = 
      tsimp[Symmetrize[
          Symmetrize[phi2bucoef[[3]], {aua, aub}], {auc, 
            aud}]]

phi2bu = 
      phi2bucoef[[1]]  +  
        o phi2bucoef[[2]] ru[alc]  +  (o^2) 
          phi2bucoef[[3]] ru[alc] ru[ald]

Further simplification is possible. Here the phi2bu is essentially the same as above although the warning messages appear.

tsimpp[phi2bu]

phi2bu = %;

the normal vector

The elements of the normal vector are denoted as ((Subscript[B, a])^dim), a=1,...,dim. But for the moment, we use "norma" for ((Subscript[B, a'])^dim), and "normb" for ((Subscript[B, dim])^dim). First of all, we assume the following expressions of these values using unknown na2, na3, nb2, nb3.

DefineTensor[tna2, "na2", {{1, 2}, 1}]

DefineTensor[tna3, "na3", {{1, 3, 2}, 1}]

norma  = 
      o  tna2[ala, aub] ru[alb] + (o^2) tna3[ala, aub, auc] ru[alb] 
          ru[alc]

RuleUnique[rule5, tB[ala_, 9], norma, IndexaQ[ala]]

DefineTensor[tnb1, "nb1", {{1}, 1}]

DefineTensor[tnb2, "nb2", {{2, 1}, 1}]

normb  = 
      1 + o  tnb1[aub] ru[alb] + (o^2) tnb2[aub, auc] ru[alb] 
          ru[alc]

RuleUnique[rule6, tB[-9, 9], normb]

The inner product of the normal vector (((Subscript[B, a]))^dim) and a tangent vector ((Subscript[B, a])^a)' is foo7=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^a'\)(u)\) \(\(B\_q\^dim\)(u)\)\).

phi2u tB[la, aua] tB[lb, 9]

CanAll[sepa[sepa[%, la], lb]]

foo7 = tgeto2[ApplyRules[%, {rule3, rule4, rule5, rule6}]]

The squared norm of the normal vector (((Subscript[B, a]))^dim)  is foo8= \(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^dim\)(u)\) \(\(B\_q\^dim\)(u)\)\).

phi2u tB[la, 9] tB[lb, 9]

CanAll[sepa[sepa[%, la], lb]]

foo8 = tgeto2[ApplyRules[%, {rule3, rule4, rule5, rule6}]]

Now, we will solve na2, na3, nb2, nb3 from the equations foo7==0, foo8==1. First we get the coefficients in foo7 and foo8, and relabel nonfree indexes.

foo9 = 
      Simplify[CoefficientList[foo7 /. {ru[al_] -> x}, x]/{1, o, 
              (o^2)}] /. {au1 -> aub, au2 -> auc}

foo10 = 
      Simplify[CoefficientList[foo8 - 1 /. {ru[al_] -> x}, x]/{1, o, 
              (o^2)}] /. {au1 -> aub, au2 -> auc}

aoo = Solve[foo9[[2]] == 0, tna2[aua, aub]]

RuleUnique[rule11, tna2[aua_, aub_], aoo[[1, 1, 2]], 
      IndexaQ[aua] ∧ IndexaQ[aub]]

aoo = Solve[foo10[[2]] == 0, tnb1[aub]]

RuleUnique[rule12, tnb1[aub_], aoo[[1, 1, 2]], 
      IndexaQ[aub]]

aoo = 
      Solve[ApplyRules[foo9[[3]], {rule11, rule12}] == 0, 
        tna3[aua, aub, auc]]

aoo2 = tsimp[Symmetrize[aoo[[1, 1, 2]], {aub, auc}]]

RuleUnique[rule13, tna3[aua_, aub_, auc_], aoo2, 
      IndexaQ[aua] ∧ IndexaQ[aub] ∧ IndexaQ[auc]]

aoo = 
      Solve[ApplyRules[foo10[[3]], {rule11, rule12, rule13}] == 0, 
        tnb2[aub, auc]]

aoo2 = tsimp[Symmetrize[aoo[[1, 1, 2]], {aub, auc}]]

aoo2 = tsimpp[%]

RuleUnique[rule14, tnb2[aub_, auc_], aoo2, 
      IndexaQ[aub] ∧ IndexaQ[auc]]

ApplyRules[tnb2[aua, aub], rule14]

We get foo15=((Subscript[B, a'])^dim)(u) and foo16= ((Subscript[B, dim])^dim)(u) below.

norma

foo15 = 
      tgeto2[ApplyRules[norma, {rule11, rule12, rule13, rule14}]]

foo15 = tsimpp[foo15]

RuleUnique[rule15, tB[ala_, 9], foo15, IndexaQ[ala]]

ApplyRules[tB[ala, 9], rule15]

normb

foo16 = 
      tgeto2[ApplyRules[normb, {rule11, rule12, rule13, rule14}]]

RuleUnique[rule16, tB[-9, 9], foo16]

ApplyRules[tB[-9, 9], rule16]

In the below, we confirm if the normal vector is orthogonal to the tangent vectors and if the length of the normal vector is 1.

The inner product of the normal vector (((Subscript[B, a]))^dim) and a tangent vector ((Subscript[B, a])^a)' is foo17=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^a'\)(u)\) \(\(B\_q\^dim\)(u)\)\).

phi2u tB[la, aua] tB[lb, 9]

CanAll[sepa[sepa[%, la], lb]];

foo17 = 
      tgeto2[ApplyRules[%, {rule3, rule4, rule15, rule16}]]

The squared norm of the normal vector (((Subscript[B, a]))^dim)  is foo18= \(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η(
              u)\)\)\(\(B\_p\^dim\)(u)\) \(\(B\_q\^dim\)(u)\)\).

phi2u tB[la, 9] tB[lb, 9]

CanAll[sepa[sepa[%, la], lb]];

foo18 = 
      tgeto2[ApplyRules[%, {rule3, rule4, rule15, rule16}]]

tsimpp[foo18]

(u,v)-coordinate system

foo20 = (Subscript[η, a])(u, 
          v) = (Subscript[η, a])(u) + ((Subscript[B, a])^dim)(u) v.

foo20 = se[la] + tB[la, 9] v

foo21 = (Subscript[η, a'])(u, v).

foo21 = 
      ApplyRules[
        foo20 /. {la -> ala}, {rule1, rule2, rule15, rule16}]

RuleUnique[rule21, se[ala_], foo21, IndexaQ[ala]]

foo22 = (Subscript[η, dim])(u, v).

foo22 = 
      ApplyRules[
        foo20 /. {la -> -9}, {rule1, rule2, rule15, 
          rule16}]

RuleUnique[rule22, se[-9], foo22]

RuleUnique[rule22b, se[9], foo22]

change of variables

The Jacobian of the change of variables η↔(u,v) is J=det((( ∂ η(u, v))/( ∂(u, v)))). The asymptotic expression of  log det J is obtained up to O((n^(-1))) term in "logdetJ". The density function f(u, v | v0) is obtained from f(y | η) as shown in "logdensityuv", where the parameter value is specified as η = (0,  ... , 0, v0).

Jacobian

Here we will calculate the log of the Jacobian J=det((( ∂ (Subscript[η, a])(u, v))/( ∂(u, v)))).

First, we obtain the expression of (({{(( ∂ (Subscript[η, a'])(u, v))/( ∂(Subscript[u, b']))), (( ∂ (Subscript[η, a'])(u, v))/( ∂v))}, {(( ∂ (Subscript[η, dim])(u, v))/( ∂(Subscript[u, b']))), (( ∂ (Subscript[η, dim])(u, v))/( ∂v))}}))=(({{foo23, foo25}, {foo24, foo26}})).

foo23 =(( ∂ (Subscript[η, a'])(u, v))/( ∂(Subscript[u, b'])))

foo23 = tsimp[difa[foo21, ru, aub]]

foo24 =(( ∂ (Subscript[η, dim])(u, v))/( ∂(Subscript[u, b'])))

foo24 = tsimp[difa[foo22, ru, aub]]

foo25 =(( ∂ (Subscript[η, a'])(u, v))/( ∂v))

foo25 = tsimp[D[foo21, v]]

foo26 =(( ∂ (Subscript[η, dim])(u, v))/( ∂v))

foo26 = tsimp[D[foo22, v]]

foo27 = (foo25/foo26)

foo27 = tsimpp[tgets2[y/x, x, foo26, y, foo25]]

foo28=foo27 foo24 =(( foo24 foo25)/foo26)

foo28 = tsimpp[tgets2[x y, x, foo27, y, foo24]]

foo29=foo23-foo28-((Subscript[δ, a'])^b)'=foo23-(( foo24 foo25)/foo26)-((Subscript[δ, a'])^b)'.

foo29 = 
      tsimp[tgets2[x -  y, x, foo23, y, foo28]] - 
        Kdelta[ala, aub]

Now, (({{foo23, foo25}, {foo24, foo26}})) is transformed to (({{I + foo29, 0}, {foo24, foo26}})) by the simple conventions preserving the determinant.

Then, log J = 
        log det(I + foo29)  +  log(1+(foo26-1)) = trace(foo29) - (1/2) trace((foo29^2))  +  
        foo26 - 1  -  (1/2) (((foo26 - 1))^2).

foo30=trace(foo29)

foo30 = tsimpp[Kdelta[aua, alb] foo29]

foo31 = trace((foo29^2))

foo31 = 
      tsimp[tgets2[x y, x, foo29, y, 
          foo29 /. {ala -> aua, aub -> alb}]]

foo32 = (((foo26 - 1))^2)

foo32 = tsimp[tgets2[(x^2), x, foo26 - 1]]

Finally,  logdetJ=log det J = foo30 - (1/2)foo31 +foo26 - 1 - (1/2) foo32.

logdetJ = 
      tgeto2[tsimpp[
          foo30 - (1/2) foo31 + foo26 - 1 - (1/2) foo32]]

Collect[logdetJ /. {ru[al_] -> u}, {o, v, u}, Simplify]

density function f(u,v|v0)

log f(y | η) is given as follows.

tgeto2[logdensityy]

First, we set η=(0,...,0,v0). In other words, (Subscript[η, a']) = 0, and (Subscript[η, dim]) = v0. Here we use symbol v0 for λ.

RuleUnique[rule35, se[la_], v0 Kdelta[la, 9]]

foo35 = tsimpp[ApplyRules[logdensityy, rule35]]

Then, change of variables y=η(u,v).

foo36 = 
      tsimpp[CanAll[
          sepa[CanAll[sepa[CanAll[sepa[CanAll[sepa[foo35, l1]], l1]], l1]], 
            l1]]]

foo37 = 
      tsimpp[tgeto2[
            ApplyRules[foo36 /. ry -> se, {rule21, rule22, rule22b}]] + 
          logdetJ]

Now we get logdensityuv=foo38=f(u,v|v0).

foo38 = 
        Collect[foo37, {o, v, ru[al1] ru[al2] ru[al3] ru[al4], 
            ru[al1] ru[al2] ru[al3], ru[al1] ru[al2], ru[al1]}];

logdensityuv = foo38

InputForm[foo38]

We extract the coefficients for u, v from the density function f(u,v|v0).

foo39 = 
        Collect[CoefficientList[foo38 /. {ru[al_] -> u}, {v, u}], {o, 
            v0}, tsimpp[
              CanAll[# /. {au1 -> aua, au2 -> aub, au3 -> auc, 
                    au4 -> aud, au5 -> aue, al1 -> ala, 
                    al2 -> alb, al3 -> alc, al4 -> ald, 
                    al5 -> ale}]] &];

Dimensions[foo39]

foo40 = foo39;

(u^0)-coefficients

(v^0) (u^0)

foo40[[1, 1]]

(v^1) (u^0)

foo40[[2, 1]]

(v^2) (u^0)

foo40[[3, 1]]

(v^3) (u^0)

foo40[[4, 1]]

(v^4) (u^0)

foo40[[5, 1]]

(u^1)-coefficients

(v^0) (u^1)

foo40[[1, 2]]

(v^1) (u^1)

foo40[[2, 2]]

(v^2) (u^1)

foo40[[3, 2]]

(v^3) (u^1)

foo40[[4, 2]]

(v^4) (u^1)

foo40[[5, 2]]

(u^2)-coefficients

(v^0) (u^2)

foo40[[1, 3]]

(v^1) (u^2)

foo40[[2, 3]]

(v^2) (u^2)

foo40[[3, 3]]

(v^3) (u^2)

foo40[[4, 3]]

(v^4) (u^2)

foo40[[5, 3]]

(u^3)-coefficients

(v^0) (u^3)

foo40[[1, 4]]

(v^1) (u^3)

foo40[[2, 4]]

(v^2) (u^3)

foo40[[3, 4]]

(v^3) (u^3)

foo40[[4, 4]]

(v^4) (u^3)

foo40[[5, 4]]

(u^4)-coefficients

(v^0) (u^4)

foo40[[1, 5]]

(v^1) (u^4)

foo40[[2, 5]]

(v^2) (u^4)

foo40[[3, 5]]

(v^3) (u^4)

foo40[[4, 5]]

(v^4) (u^4)

foo40[[5, 5]]

(Subscript[z, c])-formula

We modify the signed distance v to obtain a modified signed distance specified by w = 
        v + (∑_(r = 0)^∞) cbr[
                r] (v^r) + (Subscript[u, a']) (∑_(r = 0)^∞) br[
                  a', r] (v^r). The density function of w and its cumulants are obtained up to O((n^(-1))) terms. The distribution function of w is calculated by applying the Cornish-Fisher expansion to the cumulants of w. We would also take account of the scaling by the factor tau as well as the local coordinates at the projection in the below.

modified signed distance w

The inverse series of the modified signed distance specifies  v = 
        w - (∑_(r = 0)^∞) cr[
                r] (w^r) - (Subscript[u, a']) (∑_(r = 0)^∞) br[
                  a', r] (w^r). We have assumed that cbr[0] and cbr[2] are O((n^(-1/2))) and cbr[1], cbr[3], and all br[a',r] are O((n^(-1))). The other coefficients are cbr[r] = O((n^(-3/2))) for r≥4. Then the same order applies to cr[r]'s. The modified signed distance w is characterized by the coefficient vector c=(cr[0],cr[1],cr[2],cr[3]) up to O((n^(-1))) terms, since we can ignore the linear term in u as explained later. The change of variable is given in "rulevinuw" for v expressed in terms of u and w. The Jacobian is given in logjvw=log ∂ (v/∂) w. The joint density of (u,w) is obtained by f(u,w|v0)=f(u,v(w,u)|v0)J , and log f(u,w|v0) is stored in "logdensityuw".  We then calculate  log f(
            w | v0) =    log ∫f(u, 
                w | v0) du as an application of "logeexppoly" to logdensityuw, and stored in "logdensityw". In fact, the linear term in u, namely (Subscript[u, a']) (∑_(r = 0)^∞) br[
              a', r] (w^r) does not contribute to the argument for deriving the distribution function of w as seen in logdensityw. By using "logeexppoly" again, we obtain the cumulants of w as shown in "cumulantw".

define the modified signed distance as a series of v.

Define w = foo45 as a function of v below.

foo44 = {o, (o^2), o, (o^2)};

foo45 = v + Sum[
          foo44[[i + 1]] cbr[i] (v^i) ,  {i, 0, 
            3} ]

Then consider the inversion v=foo46 as a function of w  below.

foo46 = 
      w -  Sum[foo44[[i + 1]] cr[i] (w^i) ,  {i, 0, 3} ]

The relations between the two sets of the coefficients are given below.

func47[cbs_] := 
      Table[gets2[
          cbs[[i + 1]] - 
            Sum[(i - j + 1) cbs[[i - j + 1 + 1]] 
                cbs[[j + 1]], {j, 0, i}]], {i, 0, 
          Length[cbs] - 2}]

foo47 = 
      func47[{ o cbr[0] ,  (o^2) cbr[1] ,  o cbr[2],  (o^2) cbr[3], 
            0}]/foo44

rule47 = Table[cr[i] -> foo47[[i + 1]], {i, 0, 3}]

The relation is actually obtained by solving the following coefficients==0.

Simplify[
      CoefficientList[gets2[foo45 /. {v -> foo46}], w]]

Checking if the relation is correct by seeing the identity.

gets2[foo45 /. {v -> foo46} /. rule47]

Consider the following O((n^(-1))) term.

func48[v_] = 
      ru[ala] (o^2) (∑_(r = 0)^∞) br[aua, r] 
            (v^r)

Since v = w + O((n^(-1/2))), func48[v] = func48[w] + O((n^(-3/2))), and we can ignore the difference between func48[v] and func48[w]. So, if we redefine w = foo45+func48[v], and v=foo46-func48[w], the inversion relation still holds.  We call "v" as the signed distance, and "w" as a modified signed distance characterized by the coefficients cr[r] and br[r].

Jacobian J = ∂ (v/∂) w of the transformation from v to w is given below. Here D[func48[w],w] is denoted as (o^2) (Subscript[u, a']) (dbr^a)'.

DefineTensor[tdbr, "dbr", {{1}, 1}]

foo49 = D[foo46, w] - (o^2) ru[ala]  tdbr[aua]

We need the log of the Jacobian for later use. logjvw=log ∂ (v/∂) w.

logjvw = 
      geto2[(foo49 - 1) - (1/2) (((foo49 - 1))^2)]

Similarly, we write func48[w] as (o^2) (Subscript[u, a']) (br^a)',and vinuw is v expressed by u and w.

DefineTensor[tbr, "br", {{1}, 1}]

vinuw  =  CanAll[foo46 - (o^2) ru[ala] tbr[aua]]

RuleUnique[rulevinuw, v, vinuw]

density function of w

First, we obtain the joint density of (u,w) i.e., f(u,w|v0)=f(u,v(w,u)|v0)J from the joint density of (u,v), the transformation v=v(w,u), and the Jacobian.

foo50 = 
        tsimpp[tgeto2[
            ApplyRules[logdensityuv, rulevinuw] + logjvw]];

This is the log of the density. logdensityuw = log f(u,w|v0).

logdensityuw = 
      Collect[foo50, {ru[al1] ru[al2] ru[al3] ru[al4], 
          ru[al1] ru[al2] ru[al3], ru[al1] ru[al2], ru[al1], o, w, 
          v0}]

This is the log of the standard multivariate normal density in dim-1 dimension.

foo51 = -((dim - 1)/2) Log[2  Pi]  -  (1/2) ru[ala] 
          ru[aua]

Get the coefficients of the polynomials of (Subscript[u, a']) from log f(u, w | v0) - foo51.

foo52 =  
        Collect[CoefficientList[(logdensityuw - 
                  foo51) /. {ru[al_] -> u}, u], o, 
          tsimpp[CanAll[# /. {au1 -> aua, au2 -> aub, 
                    au3 -> auc, au4 -> aud, au5 -> aue, 
                    al1 -> ala, al2 -> alb, al3 -> alc, 
                    al4 -> ald, 
                    al5 -> ale}]] &];

Dimensions[foo52]

constant term in terms of u

foo52[[1]]

coefficient of (Subscript[u, a'])

foo52[[2]]

coefficient of (Subscript[u, a']) (Subscript[u, b'])

foo52[[3]]

coefficient of (Subscript[u, a']) (Subscript[u, b']) (Subscript[u, c'])

foo52[[4]]

coefficient of (Subscript[u, a']) (Subscript[u, b']) (Subscript[u, c']) 
        (Subscript[u, d'])

foo52[[5]]

We now calculate log f(
            w | v0) =    log ∫f(u, 
                w | v0) du as an application of "logeexppoly" to foo52. It is denoted by logdensityw.

foo53 = Collect[foo52/{1, o, o, o, (o^2)}, o];

RuleUnique[rule53u0, ta0, foo53[[1]]]

RuleUnique[rule53u1, ta1[aua_], foo53[[2]]]

RuleUnique[rule53u2, ta2[aua_, aub_], 
      tsimp[Symmetrize[foo53[[3]], {aua, aub}]]]

RuleUnique[rule53u3, ta3[aua_, aub_, auc_], 
      tsimp[Symmetrize[foo53[[4]], {aua, aub, auc}]]]

RuleUnique[rule53u4, ta4[aua_, aub_, auc_, aud_], 
      tsimp[Symmetrize[foo53[[5]], {aua, aub, auc, aud}]]]

foo54 = 
        tgeto2[ApplyRules[logeexppoly /. {l1 -> al1, u1 -> au1, 
                  l2 -> al2, u2 -> au2, l3 -> al3, 
                  u3 -> au3, l4 -> al4, 
                  u4 -> au4} /. {sb[la_] -> 0}, {rule53u0, 
              rule53u1, rule53u2, rule53u3, rule53u4}]];

This is log f(w | v0).

logdensityw = Collect[foo54, {w, o, v0}, tsimpp]

logdensityw // InputForm

The coefficients of (w^i)in log f(w | v0)are shown below.

foo55 = 
        Collect[CoefficientList[logdensityw, w], {o, v0}, 
          tsimpp];

Coefficient of (w^0).

foo55[[1]]

Coefficient of (w^1).

foo55[[2]]

Coefficient of (w^2).

foo55[[3]]

Coefficient of (w^3).

foo55[[4]]

Coefficient of (w^4).

foo55[[5]]

cumulants of w

Consider the normal density with mean v0+t and variance 1 for w. The log of the density is foo60.

foo60 = ((-1)/2) 
            Log[2  Pi] - (1/2) (((w - (v0 + t)))^2);

Then,we define foo61 by  (e^(w t)) f(
            w | v0) = (e^foo60) exp( foo61). This foo61 is a polynomial function of w.

foo61 = 
      Collect[w t  +  logdensityw  -  foo60, {w, o}, 
        Simplify]

Get the coefficients of (w^i)'s for foo61, and store them in foo62 below.

foo62 = 
        Collect[CoefficientList[foo61, w], {o, v0}, tsimpp];

Length[foo62]

foo62[[1]]

foo62[[2]]

foo62[[3]]

foo62[[4]]

foo62[[5]]

Apply "logeexppoly" to foo62.  We get foo64= \(TraditionalForm\`log \(∫\_\(-∞\)\%∞ e\^\(w\ t\)\ \(f(w | v0)\) dw\) = \ 
        log \(∫\_\(-∞\)\%∞\ \(e\^foo60\) \(exp(\ foo61)\)\ \(\(dw\)\(.\)\)\)\)

foo63 = Simplify[foo62/{1, o, o, o, (o^2)}];

RuleUnique[rule63u0, ta0, foo63[[1]]]

RuleUnique[rule63u1, ta1[aua_], foo63[[2]]]

RuleUnique[rule63u2, ta2[aua_, aub_], foo63[[3]]]

RuleUnique[rule63u3, ta3[aua_, aub_, auc_], foo63[[4]]]

RuleUnique[rule63u4, ta4[aua_, aub_, auc_, aud_], 
      foo63[[5]]]

foo64 = 
        tgeto2[ApplyRules[
            logeexppoly /. {sb[la_] -> v0 + t}, {rule63u0, rule63u1, 
              rule63u2, rule63u3, rule63u4}]];

Get the coefficients of (t^i), i = 0, 1, 2, 3, 4,  and multiply i! so that we get (Subscript[κ, i]).

foo65 = 
        Collect[CoefficientList[foo64, t]* {1, 1, 2, 6, 24}, {o, v0}, 
          tsimpp];

Length[foo65]

(Subscript[κ, 0]) should be zero.

foo65[[1]]

foo65[[1]] = 
      CanAll[% /. {al1 -> alc, au1 -> auc, al2 -> alb, 
            au2 -> aub}]

(Subscript[κ, 1])

foo65[[2]]

(Subscript[κ, 2])

foo65[[3]]

(Subscript[κ, 3])

foo65[[4]]

(Subscript[κ, 4])

foo65[[5]]

cumulantw = Drop[foo65, 1];

cumulantw // InputForm

distribution function of w

The Cornish-Fisher expansion for the standardized random variable is taken from Johnson and Kotz (1994) as shown in "cfexpx" below. We first obtain the same expansion for nonstandardized variable as shown in "cfexpw", and apply it to the cumulants of w. This gives the distribution function of w, and we obtain zformula = (Φ^(-1))(
          Pr {W ≤ w | v0}), where (Φ^(-1)) is the quantile function of the standard normal distribution. The same expression, but without MathTensor notation, is also given in "zform". Finally, the scaling by the factor "tau" is applied to these results, and zc-formula is stored in "zformulatau" as well as in "zformtau".

Cornish-Fisher expansion  (p.66 of Johnson and Kotz 1994).

cfexpx = U((Subscript[X, a])) in p.66 of JK94 is for the standardized distribution. We assume the cumulants are kx1 = 0, kx2 = 1, kx3 = O((n^(-1/2))),  
      kx4 = O((n^(-1))), kx5, 
      kx6, ...=O((n^(-3/2))).  The following expression is defined by Pr {X  ≤  x}  =  Φ(
          cfexpx), or cfexpx = (Φ^(-1))(
          Pr {X ≤ x}).

cfexpx = 
        x  -  (1/6) ((x^2) - 1) kx[3] - (1/24) ((x^3) - 3  x) 
            kx[4]  +  (1/36) (4  (x^3) - 7  x) kx[3]^2;

For w,  the cumulants are kw1 = O(1), 
              kw2 = 1 + O((n^(-1/2))), kw3=O((n^(-1/2))), kw4 = O((n^(-1))). We first standardize w and apply cfexpx to the standardized w to get cfexpw = (Φ^(-1))(
          Pr {W ≤ w}).

rule70 = {x -> ((w - kw[1])/Sqrt)[kw[2]]};

rule71 = {kx[1] -> 0, kx[2] -> 1, 
          kx[3] -> kw[3]/kw[2]^(3/2), 
          kx[4] -> kw[4]/kw[2]^2};

rule72 = {kw[1] -> aw1, kw[2] -> 1 + o aw2, 
          kw[3] -> o aw3, kw[4] -> (o^2) aw4};

rule73 = 
      Simplify[Solve[(x /. rule70) == x, 
            w][[1]]]

cfexpw = gets2[cfexpx /. Join[rule70, rule71] /. rule72]

Cornish-Fisher expansion of w

We apply cfexpw to cumulantw.

aw = 
        Simplify[(cumulantw - {0, 1, 0, 0})/{1, o, o, 
              (o^2)}];

RuleUnique[rulecumaw1, aw1, aw[[1]]]

RuleUnique[rulecumaw2, aw2, aw[[2]]]

RuleUnique[rulecumaw3, aw3, aw[[3]]]

RuleUnique[rulecumaw4, aw4, aw[[4]]]

func75[exp_, rule_] := tgeto2[ApplyRules[exp, rule]]

foo75 = 
        func75[func75[func75[func75[cfexpw, rulecumaw1], rulecumaw2], 
            rulecumaw3], rulecumaw4];

The following zformula is defined as zformula = (Φ^(-1))(
          Pr {W ≤ w}).

zformula = Collect[foo75, {o, w, v0}, tsimpp]

zformula // InputForm

Get the coefficients of (w^i) (v0^j) for zformula.

foo76 = 
        Collect[CoefficientList[zformula, {w, v0}], o, tsimpp];

Dimensions[foo76]

(w^0) (v0^0)

foo76[[1, 1]]

(w^0) (v0^1)

foo76[[1, 2]]

(w^0) (v0^2)

foo76[[1, 3]]

(w^0) (v0^3)

foo76[[1, 4]]

(w^1) (v0^0)

foo76[[2, 1]]

(w^1) (v0^1)

foo76[[2, 2]]

(w^1) (v0^2)

foo76[[2, 3]]

(w^1) (v0^3)

foo76[[2, 4]]

(w^2) (v0^0)

foo76[[3, 1]]

(w^2) (v0^1)

foo76[[3, 2]]

(w^2) (v0^2)

foo76[[3, 3]]

(w^2) (v0^3)

foo76[[3, 4]]

(w^3) (v0^0)

foo76[[4, 1]]

(w^3) (v0^1)

foo76[[4, 2]]

(w^3) (v0^2)

foo76[[4, 3]]

(w^3) (v0^3)

foo76[[4, 4]]

zc-formula using a simplified notation

In the below, the tensor symbols are replaced by regular symbols.

RuleUnique[rule80a, td[al1_, au1_], Daa, PairaQ[al1, au1]]

RuleUnique[rule80b, td[al1_, au2_] td[al2_, au1_], Dab2, 
      PairaQ[{al1, al2}, {au1, au2}]]

RuleUnique[rule80c, tp3[9, 9, 9], P999]

RuleUnique[rule80d, tp4[9, 9, 9, 9], P9999]

RuleUnique[rule80e, tp3[9, 9, al1_] tp3[9, 9, au1_], P99a2, 
      PairaQ[al1, au1]]

RuleUnique[rule80f, td[au1_, au2_] tp3[9, al1_, al2_], DabP9ab, 
      PairaQ[{al1, al2}, {au1, au2}]]

RuleUnique[rule80g, tp3[9, al1_, au2_] tp3[9, al2_, au1_], P9ab2, 
      PairaQ[{al1, al2}, {au1, au2}]]

RuleUnique[rule80h, tp4[9, 9, al1_, au1_], P99aa, 
      PairaQ[al1, au1]]

rule81 = {cr[0] -> c0, cr[1] -> c1, cr[2] -> c2, 
          cr[3] -> c3};

The following expression of "zform" is equivalent to "zformula", but without the tensor symbols so that we can use it without MathTensor.

zform = 
      Collect[ApplyRules[
            zformula, {rule80a, rule80b, rule80c, rule80d, rule80e, rule80f, 
              rule80g, rule80h}] /. rule81, {o, w, v0}]

zform // InputForm

scaling by the factor "tau".

When scaling the problem by the factor τ, the expression of the z-formula changes. First, the O(1) term such as w and v0 are multiplied by (τ^(-1)). On the other hand, O((n^(-1/2))) terms such as (ϕ^ijk)and (d^ab)are multiplied by τ, O((n^(-1)))terms such as (ϕ^ijkl)and (e^abc)are by (τ^2)  The cr[r] coefficient for modified signed distance is multiplied by (τ^(r - 1)).

rule85a = {w -> w/tau, v0 -> v0/tau};

rule85b = {tp3[ala_, alb_, alc_] -> tau tp3[ala, alb, alc], 
          td[ala_, alb_] -> tau td[ala, alb]};

rule85c = {tp4[ala_, alb_, alc_, ald_] -> (tau^2) 
              tp4[ala, alb, alc, ald], 
          te[ala_, alb_, alc_] -> (tau^2) 
              te[ala, alb, alc]};

rule85d = {cr[r_] -> cr[r] (tau^(r - 1))};

rule85 = Join[rule85a, rule85b, rule85c, rule85d];

zformulatau = Collect[zformula /. rule85, {tau, o, w, v0}]

zformulatau // InputForm

zformtau = 
      Collect[ApplyRules[
            zformulatau, {rule80a, rule80b, rule80c, rule80d, rule80e, 
              rule80f, rule80g, rule80h}] /. rule81, {tau, o, w, 
          v0}]

zformtau // InputForm

local coordinates at the projection

We consider a local coordinate Δη around a point η(u0,0) on the surface, where u0 indicates any specified value of u.  This will be used for u0 specifying the projection of y onto the surface. The change of variable η↔Δη  is specified by  (Subscript[η, a]) = (Subscript[η, a])(u0, 
            0)  +  ((Subscript[B, a])^b)(u0) Δ(Subscript[η, b]) for each u0.  The expression for (Subscript[η, a']) is given in "rule93", and that for (Subscript[η, 9]) is in "rule94". The surface is expressed as \(TraditionalForm\`Δη\_dim = \(-\ 
              d\&^\^\(a' b'\)\) \(Δη\_\(a'\)\) Δη\_\(b'\) - \ \(e\&^\^\(a' b' 
                  c'\)\) \(Δη\_\(a'\)\) \(Δη\_\(b'\)\) Δη\_\(c'\)\), where the coefficients are shown in foo101 and foo102. Next,  the expression for \(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η = η(u0) + 
                BΔη\)\)\ \(\(B\_p\^a\)(u0)\) \(\(B\_q\^b\)(
            u0)\)\) is obtained and stored in foo114[ua,ub]. This is equated with \(TraditionalForm\`ϕ\&^\^ab + \(ϕ\&^\^abc\) Δη\_c + \(1\/2\) \(ϕ\&^\^abcd\) \(Δη\_c\) Δη\_d\), and the coefficients  \(TraditionalForm\`ϕ\&^\^ab\), \(TraditionalForm\`ϕ\&^\^abc\), and \(TraditionalForm\`ϕ\&^\^abcd\) are obtained in foo121. The inverse of the metric=\(TraditionalForm\`ϕ\&^\_\(a' b'\) = \((ϕ\&^\^\(a' b'\))\)\^\(-1\)\)is in rule131, which is used for \(TraditionalForm\`foo132 = \(ϕ\&^\_\(a' b'\)\) 
          d\&^\^\(a' b'\)\) . These conversion rules are summarized in "rulesproj". zc-formula evaluated at η(u0,0) is shown in zformulau0, and that for scaling tau is in zformulatauu0.

the expression for the surface in the local coordinates

We use (Subscript[u0, a'])for the (Subscript[u, a'])-coordinate of the projection.

DefineTensor[ru0, "u0", {{1}, 1}]

Define local coordinates at the projection, and denote Δ(Subscript[η, a]). The local coordinate Δ(Subscript[η, a]) is in the (B^a)direction.

DefineTensor[re, "Δη", {{1}, 1}]

First, separate the type-a and 9 indices in the local parametrization at the projection. (Subscript[η, a]) in foo90 indicates the projection point, whereas foo90 itself indicates the (Subscript[η, a])-coordinates for a general point.

se[la] +  tB[la, ub] re[lb]

foo90 = CanAll[sepa[se[la] +  tB[la, ub] re[lb], lb]]

foo91 is foo90 for a=a'.

foo91 = ApplyRules[foo90 /. {la -> ala}, rule1]

foo92 is foo90 for a=9.

foo92 = ApplyRules[foo90 /. {la -> -9}, rule2]

Then, (B^a) is expanded by its expression. (Subscript[u, a'])is now substituted by (Subscript[u0, a']) to change the origin to the projection. Here we obtain foo93 = (Subscript[η, a']) and foo94 = (Subscript[η, 9]).

foo93 = 
      ApplyRules[foo91, {rule3, rule4, rule15, rule16}] /. 
        ru -> ru0

RuleUnique[rule93, se[ala_], foo93, IndexaQ[ala]]

foo94 = 
      ApplyRules[foo92, {rule3, rule4, rule15, rule16}] /. 
        ru -> ru0

RuleUnique[rule94, se[-9], foo94]

We will equate foo93 with foo95, and foo94 with foo96 below to derive the expression of the surface in the local coordinates. foo95 and foo96 define the surface with the origin at 0. First we consider (Subscript[η, a'])direction.

foo95 = tsimpp[ApplyRules[se[ala], rule1]]

Thus, (Subscript[u, a'])= foo93 as a function of Δ(Subscript[η, a]). We make it rule95.

RuleUnique[rule95, ru[ala_], foo93]

Consider (Subscript[η, 9])direction.

foo96 = ApplyRules[se[-9], rule2]

Using the rule95, we get an expression of foo96 in terms of Δ(Subscript[η, a])'s.

foo97 = tgeto2[ApplyRules[foo96, rule95]]

Equate foo94==foo97 to solve the expression of Δ(Subscript[η, 9]) in terms of Δ(Subscript[η, a']).

foo98 = 
        CoefficientList[tsimpp[foo94 - foo97], re[-9]];

Length[foo98]

foo98[[1]]

foo98[[2]]

The following foo99 gives an expression of Δ(Subscript[η, 9])in terms of Δ(Subscript[η, a']).

foo99 =  
      Collect[tgets2[-x/y, x, foo98[[1]], y, 
          foo98[[2]]], {re[al1] re[al2] re[al3], re[al1] re[al2], 
          ru0[al3]}]

get the coefficients of Δ(Subscript[η, a']) Δ(Subscript[η, b']) and  Δ(Subscript[η, a']) \(\(Δη\_\(b'\)\) Δη\_\(c'\)\)  for Δ(Subscript[η, 9]).

foo100 = 
      Collect[CoefficientList[foo99 /. re[ala_] -> x, 
            x] /. {au1 -> aua, au2 -> aub, au3 -> auc, 
            au4 -> aud, al1 -> ala, al2 -> alb, 
            al3 -> alc, al4 -> ald}, o, tsimpp]

This is (d^(a' b')) at the projection. We denote it  \(TraditionalForm\`d\&^\^\(a' b'\)\)=foo101.

foo101 = 
      Collect[Simplify[-foo100[[3]]/o], {ru0[al1], o}, 
        tsimp]

RuleUnique[rule101, td[aua, aub], foo101]

This is (e^(a' b' c')) at the projection. We denote it \(TraditionalForm\`e\&^\^\(a' b' c'\)\) = foo102.

foo102 = Simplify[-foo100[[4]]/(o^2)]

Now the surface is expressed in the local coordinates as \(TraditionalForm\`Δη\_9 = \(-\ 
              d\&^\^\(a' b'\)\) \(Δη\_\(a'\)\) Δη\_\(b'\) - \ \(e\&^\^\(a' b' 
                  c'\)\) \(Δη\_\(a'\)\) \(Δη\_\(b'\)\) Δη\_\(c'\)\).

the expressions for the potential derivatives

the metric

DefineTensor[tp2, "ϕ2", {{2, 1}, 1}]

the inverse of the metric

DefineTensor[tr2, "iϕ2", {{2, 1}, 1}]

phi2eta=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_a ∂\ η\_b\)\) is derived earlier.

phi2eta

foo110=phi2eta but (Subscript[η, p])is separated into (Subscript[η, p']) and (Subscript[η, 9]) in the summation.

foo110 = CanAll[sepa[sepa[phi2eta, l1], l2]]

foo111=foo110 but (Subscript[η, p']) and (Subscript[η, 9]) are substituted by their expressions in the local coordinates.

foo111 = tgeto2[ApplyRules[foo110, {rule93, rule94}]]

foo112=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_a ∂\ η\_b\)\( | \_\(η = η(u0) + 
                BΔη\)\)\ \(\(B\_a\^c\)(u0)\) \(\(B\_b\^d\)(
            u0)\)\)

foo112 = foo111 tB[la, uc] tB[lb, ud]

Separate the regular indices a and b for the summation into type-a indexes and 9's.

foo113 = CanAll[sepa[sepa[foo112, la], lb]];

foo114[ua,ub]=\(TraditionalForm\`\(\(∂\^2\( ϕ(η)\)\)\(\ \)\)\/\(∂\ η\_p ∂\ η\_q\)\( | \_\(η = η(u0) + 
                BΔη\)\)\ \(\(B\_p\^a\)(u0)\) \(\(B\_q\^b\)(
            u0)\)\) is the same as foo112, but the subscripts are changed. In the below, we substitute (B^a)'s by their expressions at the projection.

foo114[ua_, ub_] = foo113 /. {uc -> ua, ud -> ub}

foo115=foo114[aua,aub]

foo115 = 
      tgeto2[ApplyRules[foo114[aua, aub], {rule3, rule4, rule15, rule16}] /. 
          ru -> ru0]

check if foo115 is consistent with ph2bu obtained earlier.

tsimpp[
      tsimp[foo115 /. {re[la_] -> 0, ru0 -> ru}] - 
        phi2bu]

foo116=foo114[9,aua]

foo116 = 
      tgeto2[ApplyRules[foo114[9, aua], {rule3, rule4, rule15, rule16}] /. 
          ru -> ru0]

foo117=foo114[9,9]

foo117 = 
      tgeto2[ApplyRules[foo114[9, 9], {rule3, rule4, rule15, rule16}] /. 
          ru -> ru0]

Now obtain the coefficients in terms of Δ(Subscript[η, a]) for foo115,foo116,foo117.

goo118 = {{1, 1/o, 2/(o^2)}, {1/o, 1/(o^2), 2/(o^3)}, {2/(o^2), 2/(o^3), 
            4/(o^4)}};

goo118 // MatrixForm

foo118 = 
        Collect[CoefficientList[
              foo115 /. {re[-9] -> D9, re[ala_] -> DA}, {DA, 
                D9}] * goo118, o, 
          tsimpp[CanAll[# /. {au1 -> auc, au2 -> aud, 
                    au3 -> aue, au4 -> auf, au5 -> aug, 
                    al1 -> alc, al2 -> ald, al3 -> ale, 
                    al4 -> alf, al5 -> alg}]] &];

foo119 = 
        Collect[CoefficientList[
              foo116 /. {re[-9] -> D9, re[ala_] -> DA}, {DA, 
                D9}]*goo118, o, 
          tsimpp[CanAll[# /. {au1 -> auc, au2 -> aud, 
                    au3 -> aue, au4 -> auf, au5 -> aug, 
                    al1 -> alc, al2 -> ald, al3 -> ale, 
                    al4 -> alf, al5 -> alg}]] &];

foo120 = 
        Collect[CoefficientList[
              foo117 /. {re[-9] -> D9, re[ala_] -> DA}, {DA, 
                D9}]*goo118, o, 
          tsimpp[CanAll[# /. {au1 -> auc, au2 -> aud, 
                    au3 -> aue, au4 -> auf, au5 -> aug, 
                    al1 -> alc, al2 -> ald, al3 -> ale, 
                    al4 -> alf, al5 -> alg}]] &];

Dimensions[foo118]

Dimensions[foo119]

Dimensions[foo120]

foo121 = {foo118, foo119, foo120};

Geometric quantities at the projection

\(TraditionalForm\`ϕ\&^\^\(a' b'\)\)

foo121[[1, 1, 1]]

RuleUnique[rule121Pab, tp2[aua_, aub_], 
      foo121[[1, 1, 1]]]

\(TraditionalForm\`ϕ\&^\^\(9  a'\)\)

foo121[[2, 1, 1]]

RuleUnique[rule121P9a, tp2[9, aua_], foo121[[2, 1, 1]]]

\(TraditionalForm\`ϕ\&^\^99\)

foo121[[3, 1, 1]]

RuleUnique[rule121P99, tp2[9, 9], foo121[[3, 1, 1]]]

\(TraditionalForm\`ϕ\&^\^\(a' b' c'\)\)

foo121[[1, 2, 1]]

RuleUnique[rule121Pabc, tp3[aua_, aub_, auc_], 
      foo121[[1, 2, 1]]]

\(TraditionalForm\`ϕ\&^\^\(9  a' b'\)\)

foo121[[1, 1, 2]]

(foo121[[2, 2, 1]] /. auc -> aub) - 
      foo121[[1, 1, 2]]

RuleUnique[rule121P9ab, tp3[9, aua_, aub_], 
      foo121[[1, 1, 2]]]

\(TraditionalForm\`ϕ\&^\^\(99  a'\)\)

foo121[[2, 1, 2]]

(foo121[[3, 2, 1]] /. auc -> aua) - 
      foo121[[2, 1, 2]]

RuleUnique[rule121P99a, tp3[9, 9, aua_], 
      foo121[[2, 1, 2]]]

\(TraditionalForm\`ϕ\&^\^999\)

foo121[[3, 1, 2]]

RuleUnique[rule121P999, tp3[9, 9, 9], foo121[[3, 1, 2]]]

\(TraditionalForm\`ϕ\&^\^\(a' b' c' d'\)\)

foo121[[1, 3, 1]]

RuleUnique[rule121Pabcd, tp4[aua_, aub_, auc_, aud_], 
      foo121[[1, 3, 1]]]

\(TraditionalForm\`ϕ\&^\^\(9  a' b' c'\)\)

foo121[[1, 2, 2]]

foo121[[2, 3, 1]] /. {auc -> aub, aud -> auc}

RuleUnique[rule121P9abc, tp4[9, aua_, aub_, auc_], 
      foo121[[1, 2, 2]]]

\(TraditionalForm\`ϕ\&^\^\(99  a' b'\)\)

foo121[[1, 1, 3]]

foo121[[2, 2, 2]] /. {auc -> aub}

foo121[[3, 3, 1]] /. {auc -> aua, aud -> aub}

RuleUnique[rule121Pabcd, tp4[9, 9, aua_, aub_], 
      foo121[[1, 1, 3]]]

\(TraditionalForm\`ϕ\&^\^\(999  a'\)\)

foo121[[2, 1, 3]]

foo121[[3, 2, 2]] /. auc -> aua

RuleUnique[rule121P999a, tp4[9, 9, 9, aua_], 
      foo121[[2, 1, 3]]]

\(TraditionalForm\`ϕ\&^\^9999\)

foo121[[3, 1, 3]]

RuleUnique[rule121P9999, tp4[9, 9, 9, 9], 
      foo121[[3, 1, 3]]]

zero terms

{foo121[[1, 3, 2]], foo121[[1, 2, 3]], 
      foo121[[1, 3, 3]]}

{foo121[[2, 3, 2]], foo121[[2, 2, 3]], 
      foo121[[2, 3, 3]]}

{foo121[[3, 3, 2]], foo121[[3, 2, 3]], 
      foo121[[3, 3, 3]]}

A=\(TraditionalForm\`ϕ\&^\^\(a' b'\) - δ\^\(a' b'\)\)

DefineTensor[t125, "A", {{2, 1}, 1}]

foo125 = 
      ApplyRules[tp2[aua, aub], rule121Pab] - Kdelta[aua, aub]

RuleUnique[rule125, t125[aua_, aub_], foo125]

foo130=(A^2)

foo130 = 
      tgeto2[ApplyRules[t125[aua, auc] Kdelta[alc, ald] t125[aud, aub], 
            rule125]] /. {aua -> ala, aub -> alb}

foo131=(((I + A))^(-1)) = I - A + (A^2)=the inverse of the metric=\(TraditionalForm\`ϕ\&^\_\(a' b'\) = \((ϕ\&^\^\(a' b'\))\)\^\(-1\)\).

foo131 = 
      Collect[Kdelta[ala, alb] - ApplyRules[t125[ala, alb], rule125] + 
          foo130, o, tsimpp]

RuleUnique[rule131, tr2[ala_, alb_], foo131]

\(TraditionalForm\`foo132 = \(iϕ2\_\(a' b'\)\) 
          d\&^\^\(a' b'\)\) substitutes (Subscript[δ, a' b']) (d^(a' b'))at the projection.

tr2[ala, alb] td[aua, aub]

foo132 = 
      tsimpp[tgeto2[
            ApplyRules[o tr2[ala, alb] td[aua, aub], {rule101, rule131}]]/
          o]

RuleUnique[rule132, td[ala_, aua_], foo132, 
      PairaQ[ala, aua]]

Here we summarize the substitution rules for the projection for the derivatives. Here only discuss O((n^(-1/2)))terms. This includes (ϕ^999),  (d^(a' a')).

rulesproj = {rule121P999, rule132};

zc-formula

zformulau0 = 
      Collect[tgeto2[ApplyRules[zformula, rulesproj]], {o, ru0[al1], w, 
          v0}]

zformulatauu0 = 
      Collect[tgeto2[ApplyRules[zformulatau, rulesproj]], {tau, o, ru0[al1], 
          w, v0}]

InputForm[zformulatauu0]

We confirm that zformula depends on (Subscript[u0, a'])only linearly, and the term is only O((n^(-1))).

Collect[Coefficient[zformulatauu0, ru0[al1]], {o, tau}]

Bootstrap Methods

In this part, the asymptotic accuracies of several bootstrap methods are discussed. We use the zc-formula obtained in the previous part.

Startup

This section initializes the Mathematica session.

packages

<<  Statistics`ContinuousDistributions`

error messages

Off[General::spell1]

distribution functions

gammadist[x_, m_, α_] := 
      PDF[GammaDistribution[m, α], x]
Gammadist[x_, m_, α_] := 
      CDF[GammaDistribution[m, α], x]
f[x_] := PDF[NormalDistribution[0, 1], x]
F[x_] := CDF[NormalDistribution[0, 1], x]
Q[x_] := Quantile[NormalDistribution[0, 1], x]
Chidist[x_, {di_, nc_}] := 
      CDF[NoncentralChiSquareDistribution[di, nc], x]

Asymptotic Analysis of Bootstrap Methods

This section calculates the distribution functions of several bootstrap methods for showing their asymptotic accuracies in terms of the unbiasedness of hypothesis testing of the region R. The calculations are based on the zc-formula obtained in the previous part. We first define a pivot statistic and shows its third-order accuracy. The bootstrap probability is first-order accurate, the double bootstrap is second-order accurate, and the two-level bootstrap is second-order accurate.

preliminary

The zc-formula is given in zc[w,cc,v0,tau]=zformtau=(Subscript[z, c])(w; v0, tau) = (Φ^(-1))[
          Pr {W ≤ w; v0, tau}], where cc={c0,c1,c2,c3} specifies the modified signed distance w. The signed distance v is expressed in terms of w by  v = w - (∑_(r = 0)^3) (Subscript[c, r]) (w^r), or in the inverse series w = v - (∑_(r = 0)^3) (Subscript[cb, r]) (v^r). We write the coefficients (Subscript[cb, 0]) = cb0,  (Subscript[cb, 1]) = cb1,  or (Subscript[c, 0]) = c0,  (Subscript[c, 1]) = c1, etc. The rule to calculate {c0,c1,c2,c3} in terms of {cb0,cb1,cb2,cb3} is given in "rulecc2cb" or in cb2cc function. The function zq2cc[exp] calculates the {c0,c1,c2,c3} from the polynomial of w in terms of v.

simplification functions

Ignore O((n^(-3/2))) terms for scalar

geto2[exp_] := 
      Sum[Simplify[Coefficient[exp, o, i]] (o^i), {i, -1, 2}]

Series expansion for scalar ignoring O((n^(-3/2))) terms.

gets2[exp_] := geto2[Series[exp, {o, 0, 2}]]

geto[exp_, n_] := 
      Sum[Simplify[Coefficient[exp, o, i]] (o^i), {i, -1, n}]

gets[exp_, n_] := geto[Series[exp, {o, 0, n}], n]

zc-formula

zc[w,cc,v0,tau]=zformtau is (Subscript[z, c])(w; v0, tau) = (Φ^(-1))[
          Pr {W ≤ w; v0, tau}], where (Subscript[u0, a'])is assumed zero. The coefficients cc={c0,c1,c2,c3} specify the modified signed distance w. We do not need to use the more general zformulatauu0 in which (Subscript[u0, a']) ≠ 0, because the u0 terms contribute only O((n^(-3/2))) in our calculation below.

zc[w_, {c0_, c1_, c2_, c3_}, v0_, tau_] = 
      tau*(o*(-Daa - P999/6) + 
              o^2*((Dab2 - (Daa*P999)/6 + (5*P999^2)/72 - 
                          P9999/24 + P99a2/8)*
                      v0 + (Dab2 - 
                          DabP9ab + (Daa*P999)/6 + (13*P999^2)/72 - 
                          P9999/8 + P99a2/2 - P99aa/4 + P9ab2/2)*
                      w)) + (-v0 + w + 
              o*(-c0 - (P999*v0^2)/3 + (P999*v0*w)/
                      6 + (-c2 + P999/6)*w^2) + 
              o^2*(-(c0*P999*v0)/
                      6 + (P999^2/18 - P9999/8 + P99a2/8)*
                      v0^3 + (-c1 - (c0*P999)/
                            3 + (P9999/24 - P99a2/8)*v0^2)*
                      w + (-(c2*P999)/6 - P999^2/24 + P9999/24)*
                      v0*w^2 + (-c3 - (c2*P999)/3 - P999^2/72 + 
                          P9999/24)*w^3))/tau

modified signed distance

{c0,c1,c2,c3} are represented by {cb0,cb1,cb2,cb3}

rulecc2cb = {c0 -> cb0, c1 -> cb1 - 2 cb0  cb2, 
          c2 -> cb2, c3 -> -2 (cb2^2) + cb3};

coef2cb[coef_] := 
      Expand[Simplify[(coef - {0, 1, 0, 0})/{o, (o^2), o, 
              (o^2)}]]

cb2cc[{cb0_, cb1_, cb2_, cb3_}] := 
        Expand[Simplify[{cb0, cb1 - 2 cb0 cb2, 
              cb2, -2 (cb2^2) + cb3}]];

zq2cc[exp_] := 
      cb2cc[coef2cb[PadRight[CoefficientList[exp, v], 4]]]

the pivot and some existing bootstrap methods

The pivot is defined as z8[v]=\(TraditionalForm\`\(z\&^\_∞\)(0, 
          v) = \(-\(\(Φ\^\(-1\)\)(
            Pr {V ≥ v; v0 = 0))\)\)\)=zc[v,{0,0,0,0},0,1]. (Subscript[cb, r])'s and (Subscript[c, r])'s are in cbz8 and ccz8 respectively for z8[v].  We define zq[v]=\(TraditionalForm\`\(z\&^\_q\)(0, v)\)=z8[v] + q0 + q1 v + q2 (v^2)+ q3 (v^3), where the coefficients qq={q0,q1,q2,q3} specify the zq[v]. zq2qq calculates qq from any z-value. (Subscript[cb, r])'s and (Subscript[c, r])'s are in cbzq and cczq respectively for zq[v]. The distribution function of zq is obtained as Pr {zq[V] ≤ w; v0, tau}=Φ{zfzq[w,{q0,q1,q2,q3},v0,tau]}. We observe that zfzq[w,{0,0,0,0},0,1]=w, and thus the distribution function of z8 under v0=0 is Pr{Z8≤w;0,1)=Φ(w).

The bootstrap probability is \(TraditionalForm\`\(α\&~\_1\)(0, v, tau1) = 
        Pr {V ≤ 0; \ v0 = v, tau1}\) for y=η(0,v), and the corresponding z-value is z1[v,tau1]=\(TraditionalForm\`\(z\&~\_1\)(0, v, 
          tau1) = \(-\(\(Φ\^\(-1\)\)(\(α\&~\_1\)(0, v, 
              tau1))\)\)\)=-zc[0,{0,0,0,0},v,tau1]. For tau1=1, we define \(TraditionalForm\`\(\(z\&^\_0\)\(=\)\)\)z0[v]=z1[v,1], which can be regarded as another w. For general tau1, w1=tau1 z1[v,tau1] is regarded as another w with (Subscript[cb, r])'s being cbw1 and (Subscript[c, r])'s being ccw1, and the distribution function is expressed as Pr {W1 ≤ w;  v0, tau} = Φ(
          zfw1). For tau=1 and tau1=1, we have \(TraditionalForm\`Pr {Z\&^\_0 ≤ w; v0, 
            tau = 1} = Φ {zfz0[w, v0]}\), which becomes zfz0[w, 0] = w + O((n^(-1/2))) under v0=0, showing the first-order accuracy of z0[v].

The z-value of the double bootstrap probability is \(TraditionalForm\`zd[
          v] = \(-\(Φ\^\(-1\)\)[
            Pr {Z\&^\_0 ≤ \(z\&^\_0\)(v); v0 = 0}]\)\), and we observe that z8[v]=zd[v], showing the double bootstrap asymptotically equivalent to the third-order accurate pivot statistic up to O((n^(-1))) terms.

The ABC formula is given in abcformula[v,ac]. The z-value of the two-level bootstrap method is calculated in za[v]. Its (Subscript[q, i])'s are in qqza. The distribution function of za[v]  under tau=1 is Pr{za[V]≤w;v0,1}=Φ[zfzq[w,qqza,v0,1]]. This becomes w + O((n^(-1))) for v0=0, showing the second-order accuracy of the two-level bootstrap.

pivot statistic

We define \(TraditionalForm\`\(α\&^\_∞\)(0, v) = 
        Pr {V ≥ v; v0 = 0)\), and the corresponding z-value \(TraditionalForm\`\(z\&^\_∞\)(0, 
          v) = \(-\(\(\(Φ\^\(-1\)\)(\(α\&^\_∞\)(
                0, v))\)\(.\)\)\)\) We denote this z-value as z8[v]=zc[v,{0,0,0,0},0,1].

z8[v_] = Collect[zc[v, {0, 0, 0, 0}, 0, 1], {o, v}]

cbz8 = coef2cb[CoefficientList[z8[v], v]]

ccz8 = cb2cc[cbz8]

We slightly alter z8[v] and denoted zq[v]=\(TraditionalForm\`\(z\&^\_q\)(0, v)\) below.

zq[v_, {q0_, q1_, q2_, q3_}] = 
        z8[v]  +  o q0 +  (o^2) q1 v  + 
          o q2 (v^2)  +  (o^2) q3 (v^3);

Here we define a function to collect (Subscript[q, i])'s for later use

zq2qq[zz_] := 
      Expand[Simplify[
          PadRight[CoefficientList[zz - z8[v], v], 4]/{o, (o^2), o, 
              (o^2)}]]

check if this is correct.

zq2qq[zq[v, {q0, q1, q2, q3}]]

Now we continue to calculate the distribution function of zq[v]

cbzq = coef2cb[CoefficientList[zq[v, {q0, q1, q2, q3}], v]]

cczq = cb2cc[cbzq]

The distribution function of zq is obtained here. Pr {zq[V] ≤ w; v0, tau}=Φ{zfzq[w,{q0,q1,q2,q3},v0,tau]}.

zfzq[w_, {q0_, q1_, q2_, q3_}, v0_, tau_] = 
      Collect[zc[w, cczq, v0, tau], {o, w, v0}, Expand]

For tau=1, zfzq becomes

zfzq[w, {q0, q1, q2, q3}, v0, 1]

When v0=0, tau=1,  zfzq becomes

zfzq[w, {q0, q1, q2, q3}, 0, 1]

In particular, the distribution function of z8 under v0=0 is Pr{Z8≤w;0,1)=Φ(w).

zfzq[w, {0, 0, 0, 0}, 0, 1]

bootstrap probability

We define \(TraditionalForm\`\(α\&~\_1\)(0, v, tau1) = 
        Pr \(\({V ≤ 0; \ v0 = v, tau1}\)\(.\)\)\) This is the bootstrap probability for y=η(0,v).  The z-value is \(TraditionalForm\`\(z\&~\_1\)(0, v, 
          tau1) = \(-\(\(Φ\^\(-1\)\)(\(α\&~\_1\)(0, v, 
              tau1))\)\)\). This becomes z1[v,tau1]=-zc[0,{0,0,0,0},v,tau1]. For a general y=η(u,v) with u≠0, the expression changes only by the linear term in u and the difference is only O((n^(-1))).

z1[v_, tau1_] = 
      Collect[-zc[0, {0, 0, 0, 0}, v, tau1], {tau1, o, v}, 
        Expand]

The following w1=tau1 z1[v,tau1] is regarded as another w.

w1[v_, tau1_] = Collect[tau1 z1[v, tau1], {o, v, tau1}]

Here we obtain the coefficients cc for w1. The scale tau1 is specified instead of tau.

cbw1 = coef2cb[CoefficientList[w1[v, tau1], v]]

ccw1 = cb2cc[cbw1]

The distribution function of w1 under v0 and scale tau is  Pr {W1 ≤ w;  v0, tau} = Φ(
          zfw1).

zfw1[w_, tau1_, v0_, tau_] = 
      Collect[zc[w, ccw1, v0, tau], {o, w, tau1, v0, tau}, Expand]

We can use the function zq2cc to obtain ccw1 directly from w1[v].

zq2cc[w1[v, tau1]] - ccw1

We can do the same as above in another way though zq.

qqw1 = zq2qq[w1[v, tau1]]

Collect[zfzq[w, qqw1, v0, tau], {o, w, tau1, v0, tau}, Expand] - 
      zfw1[w, tau1, v0, tau]

The usual bootstrap probability is defined from  \(TraditionalForm\`\(α\&~\_1\)(0, v, tau)\) with tau=1. We denote it as \(TraditionalForm\`\(α\&^\_0\)(0, v) = \(α\&~\_1\)(0, v, 
          1)\). The corresponding z-value is denoted by \(TraditionalForm\`\(\(\(z\&^\_0\)(0, v)\)\(.\)\)\)

z0[v_] = z1[v, 1]

The distribution function of \(TraditionalForm\`\(z\&^\_0\)(u, v)\) is denoted by \(TraditionalForm\`Pr {Z\&^\_0 ≤ w; v0, 
            tau = 1} = Φ {zfz0[w, v0]}\).

zfz0[w_, v0_] = Collect[zfw1[w, 1, v0, 1], {o, w, v0}]

Under v0=0, zfz0 becomes

zfz0[w, 0]

double bootstrap

The z-value of the double bootstrap probability is \(TraditionalForm\`zd[
          v] = \(-\(\(\(Φ\^\(-1\)\)[
              Pr {Z\&^\_0 ≤ \(z\&^\_0\)(v); 
                  v0 = 0}]\)\(.\)\)\)\)

zd[v_] = Collect[geto2[zfz0[z0[v], 0]], {o, v}]

This is equivalent to the pivot z8[v]

z8[v] - zd[v]

two-level bootstrap

The ABC conversion formula of Efron (1987) and DiCiccio and Efron (1992) is

abcformula[v_, 
          ac_] = ((z0[v] - z0[0])/(1 - ac (z0[v] - z0[0]))) - 
          z0[0];

The acceleration constant "ac" is

ac = -(1/6) o P999;

The z-value of the bootstrap probability around the projection is

z0[0]

Thus the ABC formula, denoted za[v] here, becomes

za[v_] = gets2[abcformula[v, ac]]

The coefficients (Subscript[q, i])'s are obtained by

qqza = zq2qq[za[v]]

The distribution function of za[v]  under tau=1 is Pr{za[V]≤w;v0,1}=Φ[zfzq[w,qqza,v0,1]]

Collect[zfzq[w, qqza, v0, 1], {o, w, v0}, Expand]

Under v0=0, tau=1, it becomes

Collect[zfzq[w, qqza, 0, 1], {o, w}, Expand]

multistep-multiscale bootstrap method

Here we calculate the asymptotic accuracy of the three-step multiscale bootstrap method.

The pivot z8[v] is generalized to define  \(TraditionalForm\`\(z\&^\_∞\)(0, v, v0, 
          tau) = \(z8[v, v0, tau] = zc[v, {0, 0, 0, 0}, v0, tau]\)\), which reduces to z8[v,0,1]=z8[v]. We standardize z8[v,v0,tau] by w8[v,v0,tau]=z8[v,v0,tau]*isz8 + v0, where isz8 = 
        tau - (1/6) P999 tau v0  +  (1/36) (P999^2) 
            tau (v0^2), so that w8[v, v0, tau] = v + O((n^(-1/2))).  The distribution function is Pr {W8 ≤ w; v0, 
            tau} = Φ {zfz8[w, v0, tau]}. We show Pr {z8[V, v0, tau] ≤ w; v0, 
            tau} = Pr {W8  ≤  w isz8  +  v0; v0, 
              tau} = Φ {zfw8[w isz8 + v0]} = Φ(
              w). 
We obtain the inverse function of z8[v,v0,tau]=z in terms of v so that v8[z,v0,tau]=v is defined.

Let func[z]=(az+b) + rem[z], where rem[z]=(∑_(i = 0)^m) c[[i + 1]] (z^i)is of order O((n^(-1/2))) . We would like to calculate intfz=Q[(∫_(-∞)^∞) F[
              func[z]] f[z] dz], where f, F, and Q are, respectively, the density, the distribution, and quantile functions of the standard normal distribution. Let us define \(intxfg[a, b, 
          n] = \(\@\(1 + a\^2\)\/f[
                b\/\@\(1 + a\^2\)]\) \(∫\_\(-∞\)\%∞\( x\^n\) f[a\ x + b] f[x] dx\)\). Then we will show that  intfz=ddint[a,b,d,z]=\(TraditionalForm\`b\/\@\(1 + a\^2\) + ∑\+\(i = 0\)\%\(2  m + 1\)d[\([i + 1]\)] 
            intxfg[a, b, i]\/\@\(1 + a\^2\) + \(1\/2\) \(b\/\@\(1 + a\^2\)\) \((∑\+\(i = 0\)\%\(2  m + 1\)d[\([i + 1]\)] intxfg[a, b, i]\/\@\(1 + a\^2\))\)\^2\), where d=func2dd[a,b,rem,z] calculates the coefficients of the expansion rem[
            z] - (1/2) (a z + 
                b) rem[z]^2 = (∑_(i = 0)^(2  m + 1))d[[i + 
                  1]] (z^i).

Using ddint function defined above, we will calculate z2[v0, tau1, tau2]=\(\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z1[v, tau2])\) 
                  f[v, v0, tau1] dv]\)  =\(TraditionalForm\`\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z1[v8[z, v0, tau1], tau2])\) 
            f[z] dz]\). This is the z-value of the twostep-multiscale bootstrap probability.

Similarly we will calculate the z-value of the threestep-multiscale bootstrap probability defined by 
z3[v0, tau1, tau2, tau3]=\(\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z2[v, tau2, tau3])\) 
                  f[v, v0, tau1] dv]\)  =\(TraditionalForm\`\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z2[v8[z, v0, tau1], tau2, tau3])\) 
            f[z] dz]\).

The six geometric quantities (Subscript[γ, 1]),  ... , (Subscript[γ, 6]) are denoted by G1,...,G6 here. The scaling parameter (Subscript[s, 1]),  ... , (Subscript[s, 4]) are denoted by S1,...,S4 here. We define Z3G and Z8G in terms of G1,...,G6, S1,...,S4, and will show that Z3G=z3[v,tau1,tau2,tau3] and Z8G=z8[v].

a generalization of the pivot

We define a generalization of the pivot by \(TraditionalForm\`\(z\&^\_∞\)(0, v, v0, 
          tau) = \(z8[v, v0, tau] = zc[v, {0, 0, 0, 0}, v0, tau]\)\). We use z8 in the analysis of the multistep bootstrap.

z8[v_, v0_, tau_] = 
      Collect[zc[v, {0, 0, 0, 0}, v0, tau], {tau, v, o, v0}]

This reduces to the pivot when v0=0, tau=1.

Simplify[z8[v, 0, 1] - z8[v]]

We standardize z8[v,v0,tau] so that it can be regarded as w with proper coefficients. First, we find the rescaling factor.

scz8 = (1/tau) (1 + o (P999/6) v0)

isz8 = gets2[1/scz8]

The standardization of z8[v,v0,tau] up to O((n^(-1/2))) term is denoted w8[v,v0,tau]=z8[v,v0,tau]*isz8 + v0.

w8[v_, v0_, tau_] = 
      Collect[geto2[z8[v, v0, tau] isz8 + v0], {o, v, tau, v0}]

qqw8 = Collect[zq2qq[w8[v, v0, tau]], o, Simplify]

Then the distribution function of W8=w8[V,v0,tau] is Pr {W8 ≤ w; v0, 
            tau} = Φ {zfz8[w, v0, tau]}.

zfw8[w_, v0_, tau_] = 
      Simplify[geto2[zfzq[w, qqw8, v0, tau]]]

Then, Pr {z8[V, v0, tau] ≤ w; v0, 
            tau} = Pr {W8  ≤  w isz8  +  v0; v0, 
              tau} = Φ {zfw8[w isz8 + v0]}. The following equation implies that  Pr {z8[V, v0, tau] ≤ w; v0, 
            tau} = Φ(w).

geto2[zfw8[w isz8 + v0, v0, tau]]

We obtain the inverse function of z8[v,v0,tau]=z in terms of v so that v8[z,v0,tau]=v by applying the inversion formula of the modified signed distance to w8. We use v8-function in the following section.

ccw8 = Collect[zq2cc[w8[v, v0, tau]], o, Simplify];

ccw8o = ccw8 * {o, (o^2), o, (o^2)};

v8[z_, v0_, tau_] = 
      geto2[(w - geto2[Sum[ccw8o[[i]] (w^(i - 1)), {i, 4}]]) /. 
          w -> z isz8 + v0]

Check if z8[v8[z]]=z and v8[z8[v]]=v

geto2[z8[v8[z, v0, tau], v0, tau]]

geto2[v8[z8[v, v0, tau], v0, tau]]

some useful formula for normal integration

The normal distribution function Φ[x], the quantile function (Φ^(-1))[p],  and the normal density function ϕ[x].

{F[x], Q[p], f[x]}

First, we will calculate the following intxfp[a,b]

intxfp[a_, b_] = 
      FullSimplify[(∫_(-∞)^∞) F[a x + b] 
            f[x] \[DifferentialD]x, 
        aϵReals  ∧  b ∈ Reals]

Mathematica does not calculate intxfp, so we consider change of variables. First rewrite intxfp as follows.

FullSimplify[(∫_(-∞)^∞)(∫_(-∞)^(a x1 + b))f[x1] 
            f[x2] \[DifferentialD]x2 \[DifferentialD]x1, 
      aϵReals  ∧  b ∈ Reals]

Then, {x1,x2} is transformed to {y1,y2}.

foo1 = {x1 -> -((-y1 + a y2)/(Sqrt[1 + (a^2)])), 
          x2 -> -((-a y1 - 
                    y2)/(Sqrt[1 + (a^2)]))};

The range of the integration is now expressed as

foo2 = 
      FullSimplify[x2  ≤  a x1 + b /. foo1, 
        aϵReals  ∧  b ∈ Reals]

Now we get the integration intxfp[a, 
          b] = (∫_(-∞)^∞) F[a x + b] 
          f[x] dx

intxfp[a_, b_] = F[(b/(Sqrt[1 + (a^2)]))];

The following intx2f[
          n] = (∫_(-∞)^∞) (x^(2  n)) 
          f[x] dx gives (((2  
                      n)!)/(((2^n) n!))).

intx2f[n_] = 
      FullSimplify[(∫_(-∞)^∞) (x^(2  n)) 
            f[x] \[DifferentialD]x, 
        n ≥ 0 ∧ n ∈ Integers]

Table[intxf[n], {n, 10}]

\(The\ following\ intxfg[a, b, 
          n] = \(\@\(1 + a\^2\)\/f[
              b\/\@\(1 + a\^2\)]\) \(∫\_\(-∞\)\%∞\( x\^n\) f[a\ x + b] 
            f[x] dx\ will\ be\ used\ repeatedly . \ 
              This\ may\ be\ obtained\ \ by\ using\ the\ integration\ by\ parts\ and\ intx2f[n]\), \ 
    but\ Mathematica\ does\ it\ \(\(automatically\)\(.\)\)\)

foo5 = ((Sqrt[1 + (a^2)])/f)[
                (b/(Sqrt[1 + (a^2)]))] (∫_(-∞)^∞) (x^n) f[a x + b] f[x] \[DifferentialD]x;

intxfg[a_, b_, n_] = 
      Simplify[foo5, 
        n ≥ 0 ∧  n ∈ Integers ∧  
          a  ∈  Reals  ∧  b  ∈ Reals]

Simplify[Table[intxfg[a, b, n], {n, 0, 10}]]

Finally we consider some asymptotic expansions.

foo7 = gets2[((F[x + o y]  -  F[x])/f)[x]]

Thus, F[x+oy] = F[x] + f[x] (o y - (1/2) (o^2) x (y^2))=F[x] + f[x] * expandF[x+oy,x], where x=Coefficient[x+oy,o,0].

expandF[exp_, x_] := 
      gets2[(exp - x) - (1/2) (((exp - x))^2) x]

Next, note that

Simplify[
      gets2[F[x + o y  +  (1/2) (o^2) x (y^2)]] - (F[x] + 
            f[x] o y)]

Thus,  Q[F[x]+f[x] o y] = x + o y  +  (1/2) (o^2) x (y^2)=expandQ[x,y].

expandQ[x_, y_] := gets2[x +  y  +  (1/2) x (y^2)]

Combining these asymptotic expansions, we obtain the following integration. Let func[z]=(az+b) + rem[z], where rem[z]=(∑_(i = 0)^m) c[[i + 1]] (z^i)is of order O((n^(-1/2))) . We would like to calculate intfz=Q[(∫_(-∞)^∞) F[
            func[z]] f[z] dz.

The integrand F[func[z]]f[z] is  expandF[func[z],az+b]f[z] =  F[a z + b] f[z] + 
        f[a z + b] 
          f[z] (rem[z] - (1/2) (a z + b) rem[z]^2) . We denote rem[
            z] - (1/2) (a z + 
                b) rem[z]^2 = (∑_(i = 0)^(2  m + 1))d[[i + 
                  1]] (z^i). The coefficients d's are obtained from a,b, and c's.

func2dd[a_, b_, rem_, z_] := 
      Collect[CoefficientList[gets2[rem - (1/2) (a z + b) (rem^2)], 
          z], o, Simplify]

Now the integration is intfz=Q[(∫_(-∞)^∞)(F[
                a z + b] + 
              f[a z + 
                    b] (∑_(i = 0)^(2  m + 1))d[[i + 
                          1]] (z^i)) f[z] dz. Since (∫_(-∞)^∞) F[
              a z + b] f[z] dz = F[(b/(Sqrt[1 + (a^2)]))], and (∫_(-∞)^∞) (z^n) 
            f[a z + b] f[z] dz = 
        f[(b/(Sqrt[1 + (a^2)]))] intxfg[a, b, n]/(Sqrt[1 + (a^2)]), we obtain 
intfz = 
        Q[F[(b/(Sqrt[1 + (a^2)]))] + 
            f[(b/(Sqrt[1 + (a^2)]))] (∑_(i = 0)^(2  m + 1))d[[i + 
                        1]] intxfg[a, b, i]/(Sqrt[1 + (a^2)])]. Using expandQ, this becomes intfz=\(TraditionalForm\`b\/\@\(1 + a\^2\) + ∑\+\(i = 0\)\%\(2  m + 1\)d[\([i + 1]\)] 
            intxfg[a, b, i]\/\@\(1 + a\^2\) + \(1\/2\) \(b\/\@\(1 + a\^2\)\) \((∑\+\(i = 0\)\%\(2  m + 1\)d[\([i + 1]\)] intxfg[a, b, i]\/\@\(1 + a\^2\))\)\^2\).

dd2int[a_, b_, dd_, z_] := 
      Module[{intg}, 
        intg = (∑_(i = 1)^(Length[dd]))dd[[i]] 
              intxfg[a, b, i - 1]/(Sqrt[1 + (a^2)]);  
        gets2[(b/(Sqrt[1 + (a^2)])) (1 + (1/2) (intg^2)) + 
            intg]]

Now, the integration function may be written as intfz=dd2int[a,b,func2dd[a,b,rem,z],z];

two-step multiscale bootstrap

In the below, we will calculate z2[v0, tau1, tau2]=\(\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z1[v, tau2])\) 
                  f[v, v0, tau1] dv]\)  =\(TraditionalForm\`\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z1[v8[z, v0, tau1], tau2])\) 
            f[z] dz]\).

z1v8=z1[v8[z,v0,tau1],tau2]. This is regarded as a polynomial of z, and then intfz will be applied to z1v8.

z1v8 = geto2[z1[v8[z, v0, tau1], tau2]]

foo11 is the O(1) term

z1v8o0 = Coefficient[z1v8, o, 0]

z1v8ab = CoefficientList[z1v8o0, z]

We may make replacements b->z1v8ab[[1]], a->z1v8ab[[2]]  later for the normal integration with intfxg[a,b,n], i.e., z1v8o0=az+b.

Get the dd coefficients and store them in z1v8dd

z1v8dd = 
      Collect[func2dd[z1v8ab[[2]], z1v8ab[[1]], z1v8 - z1v8o0, z], o, 
        Simplify[#, tau1 > 0 ∧ tau2 > 0] &]

Length[z1v8dd]

Now the integration is z2[v0, tau1, tau2]= \(TraditionalForm\`\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z1v8)\) f[z] dz]\). .

z2[v0_, tau1_, tau2_] = 
      Collect[dd2int[z1v8ab[[2]], z1v8ab[[1]], z1v8dd, z], o, 
        FullSimplify[#, tau1 > 0 ∧ tau2 > 0] &]

Check if z2 reduces to z1 when one of the scales is zero

Simplify[Limit[z2[v, tau1, tau2], tau2 -> 0] - z1[v, tau1], 
      tau1 > 0]

Simplify[Limit[z2[v, tau1, tau2], tau1 -> 0] - z1[v, tau2], 
      tau2 > 0]

three-step multiscale bootstrap

In the below, we will calculate z3[v0, tau1, tau2, tau3]=\(\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z2[v, tau2, tau3])\) 
                  f[v, v0, tau1] dv]\)  =\(TraditionalForm\`\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z2[v8[z, v0, tau1], tau2, tau3])\) 
            f[z] dz]\).

z2v8=z2[v8[z,v0,tau1],tau2,tau3]. This is regarded as a polynomial of z, and then intfz will be applied to z2v8.

z2v8 = geto2[z2[v8[z, v0, tau1], tau2, tau3]]

foo11 is the O(1) term

z2v8o0 = Coefficient[z2v8, o, 0]

z2v8ab = CoefficientList[z2v8o0, z]

We may make replacements b->z2v8ab[[1]], a->z2v8ab[[2]]  later for the normal integration with intfxg[a,b,n], i.e., z2v8o0=az+b.

Get the dd coefficients and store them in z2v8dd

z2v8dd = 
      Collect[func2dd[z2v8ab[[2]], z2v8ab[[1]], z2v8 - z2v8o0, z], o, 
        FullSimplify[#, 
            tau1 > 0 ∧ tau2 > 0 ∧ tau3 > 0] &]

Length[z2v8dd]

Now the integration is z3[v0, tau1, tau2, tau3]= \(TraditionalForm\`\(Φ\^\(-1\)\)[∫\_\(-∞\)\%∞\( Φ(z2v8)\) f[z] dz]\). .

z3[v0_, tau1_, tau2_, tau3_] = 
      Collect[dd2int[z2v8ab[[2]], z2v8ab[[1]], z2v8dd, z], o, 
        FullSimplify[#, 
            tau1 > 0 ∧ tau2 > 0 ∧ tau3 > 0] &]

Check if z3 reduces to z2 when one of the scales is zero

Simplify[
      Limit[z3[v, tau1, tau2, tau3], tau3 -> 0] - z2[v, tau1, tau2], 
      tau1 > 0 ∧ tau2 > 0]

Simplify[
      Limit[z3[v, tau1, tau2, tau3], tau2 -> 0] - z2[v, tau1, tau3], 
      tau1 > 0 ∧ tau3 > 0]

Simplify[
      Limit[z3[v, tau1, tau2, tau3], tau1 -> 0] - z2[v, tau2, tau3], 
      tau2 > 0 ∧ tau3 > 0]

z3[v, tau1, tau2, tau3] // InputForm

simplifying z3 and z8

We define six geometric quantities;

GG = {G1 -> 
            v0 + (1/3) o P999 (v0^2) + 
              (o^2) (-((P999^2)/18) + (P9999/8) - (P99a2/8)) (v0^3), 
          G2 -> 
            o (-Daa - (P999/6)) v0 + 
              (o^2) (Dab2 - ((Daa P999)/2) + ((P999^2)/72) - (P9999/24) + 
                    (P99a2/8)) (v0^2), 
          G3 -> -(1/6) o P999 v0 + 
              (o^2) (((P999^2)/9) - (P9999/8) + (P99a2/4)) (v0^2), 
          G4 -> 
            (o^2) (-DabP9ab + ((Daa P999)/3) + ((2 (P999^2))/9) - 
                  (P9999/6) + (P99a2/2) - (P99aa/4) + (P9ab2/2)) (v0^2), 
          G5 -> 
            (o^2) (-((P999^2)/8) + (P9999/12) - (P99a2/8)) (v0^2), 
          G6 -> 
            (o^2) (-((P999^2)/8) + (P9999/24) - 
                  (P99a2/8)) (v0^2)};

We also define four scale parameters (but there are only three degrees of freedom among them).

SS = {S1 -> (1/(Sqrt[(tau1^2) + (tau2^2) + (tau3^2)])), 
          S2 -> (((tau1^2) (tau2^2) + (tau1^2) (tau3^2) + (tau2^2) (tau3^2))/((((tau1^2) + (tau2^2) + (tau3^2)))^2)), 
          S3 -> (((tau1^2) (tau2^2) (tau3^2) + (tau2^4) (tau3^2) + (tau1^4) ((tau2^2) + (tau3^2)))/((((tau1^2) + (tau2^2) + (tau3^2)))^3)), 
          S4 -> (((tau1^2) (tau2^2) (tau3^2))/((((tau1^2) + (tau2^2) + (tau3^2)))^3))};

The pivot z8[v0,1] is denoted by Z8G

z8G = ((G2 + ((G3^2)/2) + G4 + G5)/G1) + 
          G1 (1 + G3 + 4 (G3^2) + G6);

Simplify[gets2[z8G /. GG] - z8[v0]]

The three-step multiscale z-value z3[v0,tau1,tau2,tau3] is denoted by Z3G

Z3G = 
        G1 S1 (1 + G3 S2 + 4 (G3^2) (S2^2) + G5 S3 + 
                G6 S4) - ((G2 + G3 S2 + G4 S2 + 7 (G3^2) (S2^2) + 3 G5 S3 + 3 G6 S4)/(G1 S1));

Simplify[
      gets2[Z3G /. Join[GG, SS]] - z3[v0, tau1, tau2, tau3]]