Normal integral representation for the joint survival function of the cumulative sums of the components of multinomial random vectors


Abstract

This paper presents a multivariate normal integral representation for the joint survival function of the cumulative sums of the components of any multinomial random vector at interior lattice points. This result can be viewed as a multivariate analog of Equation (7) in [1], whose proof starts from the beta integral representation of binomial survival probabilities and uses Laplace’s method to improve Tusnády’s inequality. Our findings are based on a crucial relationship between the joint survival function of the cumulative sums of the components of any multinomial random vector and a Dirichlet probability over a corresponding cumulative-sum region. The main motivation is that such an explicit formula may eventually help streamline the conditional quantile-transformation arguments used in the multivariate KMT approximation of [2], a connection left for future work. We provide numerical checks of the identity for \(d = 2,3,4,5\).

Dirichlet distribution, Gaussian integral representation, Laplace’s method, multinomial distribution, multivariate normal distribution, normal integral representation

1 Introduction↩︎

One-dimensional quantile couplings between discrete distributions and Gaussian distributions are a central tool in probability and statistics. A classical example appears in the Komlós–Major–Tusnády (KMT) approximation [3][7], where the empirical distribution function is coupled with a Brownian bridge. As explained by [1], a key ingredient in that construction is the quantile coupling between a \(\mathrm{Binomial}(n,1/2)\) random variable and a normal random variable with the same mean and variance. More precisely, if \(X\sim \mathrm{Binomial}(n,1/2)\) and \(Y\sim \mathcal{N}(n/2,n/4)\), then one defines cutpoints \(-\infty = \beta_0 < \beta_1 < \cdots < \beta_n < \beta_{n+1} = \infty\) by \[\mathsf{P}(X\geq k) = \mathsf{P}(Y > \beta_k),\qquad k = 0,1,\ldots,n.\] Sharp control of the difference between the binomial quantiles and their Gaussian counterparts is the content of Tusnády-type inequalities. The original inequality of [8] and its later refinements [1], [9], [10] show that the normal approximation to these quantiles is much more accurate than what follows from a crude central limit theorem.

The starting point of [1] is not a local approximation of individual binomial probabilities, but rather an exact integral representation for the whole binomial tail: \[\label{eq:beta46integral46rep} \mathsf{P}(X\geq k) = \frac{n!}{(k-1)!(n-k)!}\int_0^{1/2} t^{k-1}(1-t)^{n-k}\hskip 1pt\mathrm{d}\hskip 0.75ptt, \qquad 1 \leq k \leq n.\tag{1}\] This beta integral representation converts the discrete tail probability into a continuous integral. It is the one-dimensional prototype for the Dirichlet integral representation developed in the present paper.

Now, let us see how the Gaussian structure enters. For \(n/2 < k < n\), let \(K = k-1\), \(N = n-1\), and \(\varepsilon= (2K-N)/N\). Define \[2H(t) = (1+\varepsilon)\ln t+(1-\varepsilon)\ln(1-t),\qquad h(s) = H\left(\frac{1-s}{2}\right)-H(1/2),\] and let \(\lambda_m\) denote the remainder term in Stirling’s formula, as in 5 below. With \[\gamma(\varepsilon) = \frac{(1+\varepsilon)\ln(1+\varepsilon)+(1-\varepsilon)\ln(1-\varepsilon)-\varepsilon^2}{2\varepsilon^4}, \qquad \varepsilon\neq 0,\] where \(\gamma(0) = 1/12\) by continuous extension, and \[a_N(\varepsilon) = \ln(1+N^{-1})+\lambda_N-\lambda_K-\lambda_{N-K}-\frac{1}{2}\ln(1-\varepsilon^2)-N\varepsilon^4\gamma(\varepsilon),\] Equation (7) of [1] can be written as \[\mathsf{P}(X\geq k) = e^{a_N(\varepsilon)}\sqrt{\frac{N}{2\pi}}\int_0^1 \exp\left\{Nh(s)-\frac{N\varepsilon^2}{2}\right\}\hskip 1pt\mathrm{d}\hskip 0.75pts.\] Near its maximum at \(s = 0\), the function \(h\) starts as \(h(s) = -\varepsilon s-s^2/2\) plus higher order terms. Thus the leading part of the integral is governed by the Gaussian expression \[\sqrt{\frac{N}{2\pi}}\int_0^\infty \exp\left\{-\frac{N(s+\varepsilon)^2}{2}\right\}\hskip 1pt\mathrm{d}\hskip 0.75pts = \overline{\Phi}(\varepsilon\sqrt{N}),\] where \(\overline{\Phi}(x) = \mathsf{P}(Z \geq x)\) for \(Z\sim \mathcal{N}(0,1)\). [1] then use this representation, together with careful Taylor bounds and inequalities for ratios of normal tail probabilities, to sharpen Tusnády’s inequality.

The goal of the present paper is to develop the corresponding integral representation in the multinomial setting, with the hope that such an explicit formula may eventually help simplify the difficult conditional quantile-transformation arguments used in the multivariate KMT approximation of [2]; see also [11]. If \(\boldsymbol{X} = (X_1,\ldots,X_d)\sim \mathrm{Multinomial}(n,\boldsymbol{p})\), the natural multivariate analog of a binomial tail is the joint survival function of the cumulative sums of the components \[\mathsf{P}(X_1+\cdots+X_i\geq k_1+\cdots+k_i, ~\forall i\in[d]),\] where \([d] \equiv \{1,\ldots,d\}\). These probabilities arise naturally when multinomial counts are constructed from a partition of the unit interval. In that construction, the event above can be expressed in terms of selected order statistics of independent uniform random variables; see 2 below. The joint density of those order statistics yields a Dirichlet-type integral over a region \(\mathcal{R}_d\) inside the simplex; see 4 below. When \(d = 1\) and \(p_1 = 1/2\), this identity reduces exactly to the beta integral representation for binomial tails in 1 .

The main result of the paper (Theorem 1) is an exact multivariate normal integral representation for this joint survival function at interior lattice points. This identity is not an asymptotic approximation. It is an exact rewriting of a Dirichlet integral, obtained by applying Stirling’s formula and completing the square in the natural covariance structure of the multinomial distribution. In this sense, it is a multivariate analog of the normal integral representation used by [1] in the binomial case to improve Tusnády’s inequality.

The rest of the paper is organized as follows. Section 2 introduces the necessary notation and derives the Dirichlet representation of the joint survival function of the cumulative sums of multinomial components. Section 3 states the multivariate normal integral representation. Section 4 gives two proofs: the first expands the logarithm of the Dirichlet density directly, while the second uses a Laplace-type decomposition of the Dirichlet integral. Section 5 provides numerical checks of the claimed identity.

2 Definitions and notation↩︎

For any integer \(d\in \mathbb{N}\), the \(d\)-dimensional simplex and its interior are defined by \[\mathcal{S}_d = \big\{\boldsymbol{s}\in [0,1]^d: \|\boldsymbol{s}\|_1 \leq 1\big\}, \qquad \mathrm{Int}(\mathcal{S}_d) = \big\{\boldsymbol{s}\in (0,1)^d: \|\boldsymbol{s}\|_1 < 1\big\},\] where \(\|\boldsymbol{s}\|_1 = \sum_{i=1}^d |s_i|\) denotes the \(\ell_1\)-norm in \(\mathbb{R}^d\). Given a set of probability weights \(\boldsymbol{p}\in \mathrm{Int}(\mathcal{S}_d)\), the \(\mathrm{Multinomial}(n,\boldsymbol{p})\) probability mass function is defined, for all \(\boldsymbol{k}\in \mathbb{N}_0^d \cap n \mathcal{S}_d\), by \[p_n(\boldsymbol{k}) = \frac{n!}{(n - \|\boldsymbol{k}\|_1)! \prod_{i=1}^d k_i!} \, p_{d+1}^{n - \|\boldsymbol{k}\|_1} \prod_{i=1}^d p_i^{k_i},\] where \(p_{d+1} = 1 - \|\boldsymbol{p}\|_1\in (0,1)\) and \(n\in \mathbb{N}\). Throughout, when \(\boldsymbol{s}\in \mathcal{S}_d\), set \(s_{d+1} = 1 - \|\boldsymbol{s}\|_1\). The covariance matrix of the multinomial distribution is well-known to be \(n \, \Sigma_{\boldsymbol{p}}\), where \(\Sigma_{\boldsymbol{p}} = \text{diag}(\boldsymbol{p}) - \boldsymbol{p} \boldsymbol{p}^{\top}\); see, e.g., [12]. From Theorem 1 of [13], it is also well known that \(\det(\Sigma_{\boldsymbol{p}}) = p_1 \dots p_d \, p_{d+1}\). The centered multivariate normal density with covariance matrix \(\Sigma_{\boldsymbol{p}}\) (the per-trial multinomial covariance matrix) is defined, for all \(\boldsymbol{x}\in \mathbb{R}^d\), by \[\phi_{\Sigma_{\boldsymbol{p}}}(\boldsymbol{x}) = \frac{\exp\big(-\boldsymbol{x}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \, \boldsymbol{x} / 2\big)}{\sqrt{(2\pi)^d \, \det(\Sigma_{\boldsymbol{p}})}}.\]

Let \(I_1 = (0,p_1]\) and \(I_j = (p_1 + \dots + p_{j-1}, p_1 + \dots + p_j]\) for all \(j\in [d]\backslash \{1\}\). If \(U_1, \dots,U_n \smash{\stackrel{\mathrm{iid}}{\sim}} \mathrm{Uniform}(0,1)\), then \[\boldsymbol{X} = (X_1,\ldots,X_d) \equiv \Big(\sum_{i=1}^n \mathbb{1}\{U_i\in I_1\},\ldots,\sum_{i=1}^n \mathbb{1}\{U_i\in I_d\}\Big)\sim \mathrm{Multinomial}(n,\boldsymbol{p}).\] For any \(\boldsymbol{k}\in \mathbb{N}^d \cap n \mathcal{S}_d\), set \(K_0 = 0\), \(K_i = k_1 + \dots + k_i\) for all \(i\in [d]\), and \(K_{d+1} = n+1\). If \(U_{(1)} \leq \dots \leq U_{(n)}\) denote the order statistics, then \[\label{eq:order46statistics46link} \begin{align} &\mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) \\[1mm] &\qquad= \mathsf{P}(U_{(K_i)} \leq p_1 + \dots + p_i, ~\forall i\in [d]) \\ &\qquad= \int_0^{p_1} \int_{u_1}^{p_1+p_2} \dots \int_{u_{d-1}}^{p_1+\dots+p_d} n! \prod_{i=1}^{d+1} \frac{(u_i - u_{i-1})^{K_i-K_{i-1}-1}}{(K_i - K_{i-1} - 1)!} \hskip 1pt\mathrm{d}\hskip 0.75ptu_1 \hskip 1pt\mathrm{d}\hskip 0.75ptu_2 \dots \hskip 1pt\mathrm{d}\hskip 0.75ptu_d, \end{align}\tag{2}\] where one defines \(u_0 = 0\) and \(u_{d+1} = 1\) in the multidimensional integral. After applying the change of variables \(s_i = u_i - u_{i-1}\) for all \(i\in [d]\), setting \(s_{d+1} = 1 - \|\boldsymbol{s}\|_1\), and writing \(k_{d+1} = K_{d+1} - K_d = (n + 1) - \|\boldsymbol{k}\|_1\), the above can be rewritten as \[\begin{align} &\mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) \\ &\qquad= \int_0^{p_1} \int_0^{(p_1 - s_1) + p_2} \dots \int_0^{\sum_{k=1}^{d-1} (p_k - s_k) + p_d} n! \prod_{i=1}^{d+1} \frac{s_i^{k_i-1}}{(k_i - 1)!} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}. \end{align}\] This identity provides a direct relationship between the joint survival function of the cumulative sums of the components of any multinomial random vector and a Dirichlet probability over the corresponding cumulative-sum region.

In turn, letting \(N = n - d\), \(J_i = k_i - 1\) for all \(i\in [d+1]\), so that \(J_{d+1} = n - \|\boldsymbol{k}\|_1\), writing \(\boldsymbol{J} = (J_1,\ldots,J_d)^{\top}\), and defining the regions \[\label{eq:def46R46d} \mathcal{R}_d = \left\{\boldsymbol{s}\in \mathcal{S}_d : (s_1,s_2,\dots,s_i)\in \left(\sum_{k=1}^i p_k\right)\mathcal{S}_i ~~\forall i\in [d]\right\}, \qquad \mathcal{R}_d^{\circ} = \mathcal{R}_d \cap \mathrm{Int}(\mathcal{S}_d),\tag{3}\] one can write \[\label{eq:multinomial46dirichlet46relation} \mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) = \int_{\mathcal{R}_d} \frac{(N + d)!}{N!} \times \frac{N!}{\prod_{i=1}^{d+1} J_i!} \prod_{i=1}^{d+1} s_i^{J_i} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}.\tag{4}\] Here \(\mathcal{R}_d^{\circ}\) denotes the part of \(\mathcal{R}_d\) lying in \(\mathrm{Int}(\mathcal{S}_d)\); it is not meant to denote the topological interior of \(\mathcal{R}_d\), since some cumulative-sum constraints may be active. Since \(\mathcal{R}_d\backslash \mathcal{R}_d^{\circ}\) has Lebesgue measure zero, either set may be used in 4 ; below we use \(\mathcal{R}_d^{\circ}\) whenever logarithms of the \(s_i\)’s appear.

3 Normal integral representation↩︎

For any integer \(m\in \mathbb{N}\), let \(\lambda_m\) denote the error term in Stirling’s approximation for \(\ln (m!)\): \[\label{eq:log46Stirling} \ln(m!) = \frac{1}{2} \ln(2\pi m) + m \ln m - m + \lambda_m,\tag{5}\] where \((12 m + 1)^{-1} \leq \lambda_m \leq (12 m)^{-1}\); see, e.g., [14].

For the normal integral representation, assume from now on that \[\boldsymbol{k}\in \mathcal{K}_{n,d} = \left\{\boldsymbol{k}\in \mathbb{N}^d \cap n \mathcal{S}_d : k_i \geq 2 ~~\forall i\in [d],~\|\boldsymbol{k}\|_1 \leq n - 1\right\}.\] This condition is equivalent to \(J_i\in \mathbb{N}\) for all \(i\in [d+1]\), and it ensures that \(N\in \mathbb{N}\) and that all logarithms below are finite. Together with the notation introduced in Section 1, define, for all \(i\in [d+1]\), \[\varepsilon_i = \frac{J_i / N - p_i}{p_i}, \qquad \widetilde{\varepsilon}_i = p_i \varepsilon_i = J_i/N - p_i,\] and write \[\boldsymbol{\varepsilon} = (\varepsilon_1,\ldots,\varepsilon_{d+1}), \qquad \widetilde{\boldsymbol{\varepsilon}} = (\widetilde{\varepsilon}_1,\ldots,\widetilde{\varepsilon}_d)^{\top}.\] Also, set \[\Lambda_N = \lambda_N - \sum_{i=1}^{d+1} \lambda_{J_i},\] and \[\begin{align} \Delta_N &= \ln \left\{\frac{(N + d)!}{N! N^d}\right\} + \Lambda_N - \frac{1}{2} \sum_{i=1}^{d+1} \ln (1 + \varepsilon_i) - N \widetilde{\gamma}(\boldsymbol{\varepsilon}), \tag{6} \\ \widetilde{\gamma}(\boldsymbol{\varepsilon}) &= \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln (1 + \varepsilon_i) - \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}}, \tag{7} \\ \gamma^{\star}(\boldsymbol{s}) &= \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \left(\frac{s_i}{p_i}\right) - \left\{\widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{p})^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p})\right\}. \tag{8} \end{align}\]

Theorem 1 below expresses the joint survival function of the cumulative sums of the components of any multinomial random vector at interior lattice points \(\boldsymbol{k}\in \mathcal{K}_{n,d}\) in terms of a multivariate normal integral. The result can be seen as a multivariate analog of Eq. (7) of [1], who improved on Tusnády’s inequality. For an overview of the literature on Tusnády’s inequality and the most recent improvement in the bulk, see [9] and [10], respectively.

Theorem 1 (Normal integral representation). For all \(\boldsymbol{k}\in \mathcal{K}_{n,d}\), one has \[\mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) = e^{\Delta_N} \int_{\mathcal{R}_d^{\circ}} \exp\left\{N \gamma^{\star}(\boldsymbol{s})\right\} N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}}\{N^{1/2} (\boldsymbol{p} - \boldsymbol{s} + \widetilde{\boldsymbol{\varepsilon}})\} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}.\] When \(\max_{i\in [d+1]} |\varepsilon_i| \leq \eta < 1\), a useful expansion for \(\widetilde{\gamma}(\boldsymbol{\varepsilon})\) can be found in 12 . An exact alternative expression for \(\gamma^{\star}(\boldsymbol{s})\) can be found in 13 .

Remark 1. The restriction \(\boldsymbol{k}\in \mathcal{K}_{n,d}\) is used only for the normal representation above. The Dirichlet identity 4 remains valid when some \(J_i=0\), with the convention \(0! = 1\), but those boundary cases are not covered by Theorem 1 because the logarithms in 6 , 7 , and 8 would not all be finite.

Remark 2. Theorem 1 also complements the local limit theorem for the multinomial distribution developed independently by [15] and [16].

4 Proofs↩︎

First proof of Theorem 1. Under the assumption \(\boldsymbol{k}\in \mathcal{K}_{n,d}\), one has \(J_i\in \mathbb{N}\) for all \(i\in [d+1]\) and \(\sum_{i=1}^{d+1} J_i = N\). For \(\boldsymbol{s}\in \mathcal{R}_d^{\circ}\), by taking the logarithm of the integrand in 4 and applying Stirling’s formula 5 , one obtains \[\label{eq:log46p46before} \begin{align} &\ln \left\{\frac{(N + d)!}{N!} \times \frac{N!}{\prod_{i=1}^{d+1} J_i!} \prod_{i=1}^{d+1} s_i^{J_i}\right\} \\ &\quad= \ln \left\{\frac{(N + d)!}{N!}\right\} + \ln (N!) - \sum_{i=1}^{d+1} \ln (J_i!) + \sum_{i=1}^{d+1} J_i \ln s_i \\ &\quad= \ln \left\{\frac{(N + d)!}{N! N^d}\right\} - \frac{1}{2} \ln \left\{\left(\frac{2\pi}{N}\right)^d \, \prod_{i=1}^{d+1} (J_i/N)\right\} + \sum_{i=1}^{d+1} J_i \ln \left(\frac{s_i}{J_i / N}\right) + \lambda_N - \sum_{i=1}^{d+1} \lambda_{J_i} \\ &\quad= \ln \left\{\frac{(N + d)!}{N! N^d}\right\} - \frac{1}{2} \sum_{i=1}^{d+1} \ln (1 + \varepsilon_i) - \frac{1}{2} \ln \left\{\left(\frac{2\pi}{N}\right)^d \, \prod_{i=1}^{d+1} p_i\right\} \\ &\qquad+ N \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \left(\frac{s_i}{p_i}\right) + N \sum_{i=1}^{d+1} (J_i / N) \ln \left(\frac{p_i}{J_i / N}\right) + \Lambda_N. \end{align}\tag{9}\] By the definition of \(\widetilde{\gamma}(\boldsymbol{\varepsilon})\), one can write exactly \[\label{eq:entropy461} \sum_{i=1}^{d+1} (J_i / N) \ln \left(\frac{p_i}{J_i / N}\right) = - \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} - \widetilde{\gamma}(\boldsymbol{\varepsilon}).\tag{10}\] If, in addition, \(\max_{i\in [d+1]} |\varepsilon_i| \leq \eta < 1\), then the Taylor expansion \[(1 + x) \ln(1 + x) = x + \frac{x^2}{2} - \frac{x^3}{6} + \frac{x^4}{12} + \mathcal{O}_{\eta}(x^5), \quad |x| \leq \eta < 1,\] and the fact that \(\widetilde{\varepsilon}_{d+1} = -\sum_{i=1}^d \widetilde{\varepsilon}_i\) show that \(\widetilde{\gamma}(\boldsymbol{\varepsilon})\) can be expanded as follows: \[\begin{align} \widetilde{\gamma}(\boldsymbol{\varepsilon}) &\equiv \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln (1 + \varepsilon_i) - \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} \\ &= \frac{1}{2} \sum_{i,j=1}^d \widetilde{\varepsilon}_i \widetilde{\varepsilon}_j \left\{\frac{1}{p_i} \mathbb{1}\{i = j\} + \frac{1}{p_{d+1}}\right\} - \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} - \frac{1}{6} \sum_{i,j,k=1}^d \widetilde{\varepsilon}_i \widetilde{\varepsilon}_j \widetilde{\varepsilon}_k \left\{\frac{1}{p_i^2} \mathbb{1}\{i = j = k\} - \frac{1}{p_{d+1}^2}\right\} \\ &\quad+ \frac{1}{12} \sum_{i,j,k,\ell=1}^d \widetilde{\varepsilon}_i \widetilde{\varepsilon}_j \widetilde{\varepsilon}_k \widetilde{\varepsilon}_{\ell} \left\{\frac{1}{p_i^3} \mathbb{1}\{i = j = k = \ell\} + \frac{1}{p_{d+1}^3}\right\} + \mathcal{O}_{d,\boldsymbol{p},\eta}\big(\|\boldsymbol{\varepsilon}\|_1^5\big). \end{align}\] Moreover, by [13], it is known that \((\Sigma_{\boldsymbol{p}}^{-1})_{ij} = p_i^{-1} \mathbb{1}\{i = j\} + p_{d+1}^{-1}\) for all \(i,j\in [d]\), so \[\label{eq:gamma46tilde46simplified} \frac{1}{2} \sum_{i,j=1}^d \widetilde{\varepsilon}_i \widetilde{\varepsilon}_j \left\{\frac{1}{p_i} \mathbb{1}\{i = j\} + \frac{1}{p_{d+1}}\right\} - \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} = 0.\tag{11}\] Under this additional condition, using the last two equations, one can rewrite \(\widetilde{\gamma}(\boldsymbol{\varepsilon})\) as follows: \[\label{eq:gamma46epsilon46alternative} \begin{align} \widetilde{\gamma}(\boldsymbol{\varepsilon}) &=- \frac{1}{6} \sum_{i,j,k=1}^d \widetilde{\varepsilon}_i \widetilde{\varepsilon}_j \widetilde{\varepsilon}_k \left\{\frac{1}{p_i^2} \mathbb{1}\{i = j = k\} - \frac{1}{p_{d+1}^2}\right\} \\ &\quad+ \frac{1}{12} \sum_{i,j,k,\ell=1}^d \widetilde{\varepsilon}_i \widetilde{\varepsilon}_j \widetilde{\varepsilon}_k \widetilde{\varepsilon}_{\ell} \left\{\frac{1}{p_i^3} \mathbb{1}\{i = j = k = \ell\} + \frac{1}{p_{d+1}^3}\right\} + \mathcal{O}_{d,\boldsymbol{p},\eta}\big(\|\boldsymbol{\varepsilon}\|_1^5\big). \end{align}\tag{12}\] Now, consider \(\boldsymbol{s} = (s_1,\ldots,s_d)\in \mathcal{R}_d^{\circ}\) as defined in 3 . For all \(i\in [d+1]\) and \(t\in (0,1)\), let \[\delta_i(t) = \frac{t - p_i}{p_i}.\] In particular, note that this definition yields \(\sum_{i=1}^{d+1} p_i \delta_i(s_i) = 0\). Upon applying Taylor’s formula with Lagrange remainder to the function \(x\mapsto \ln(1+x)\), there exist points \(s_i^{\star}\), with \(s_i^{\star}\) between \(s_i\) and \(p_i\), \(i\in [d+1]\), such that \[\begin{align} &\sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \left(\frac{s_i}{p_i}\right) = \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \{1 + \delta_i(s_i)\} \\ &\quad= \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \delta_i(s_i) - \frac{1}{2} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \{\delta_i(s_i)\}^2 + \frac{1}{3} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \frac{\{\delta_i(s_i)\}^3}{\{1 + \delta_i(s_i^{\star})\}^3} \\ &\quad= \sum_{i=1}^{d+1} \widetilde{\varepsilon}_i \delta_i(s_i) - \frac{1}{2} \sum_{i=1}^{d+1} p_i \{\delta_i(s_i)\}^2 - \frac{1}{2} \sum_{i=1}^{d+1} \widetilde{\varepsilon}_i \{\delta_i(s_i)\}^2 + \frac{1}{3} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \frac{\{\delta_i(s_i)\}^3}{\{1 + \delta_i(s_i^{\star})\}^3} \\ &\quad= \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{p})^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} \sum_{i=1}^{d+1} \widetilde{\varepsilon}_i \{\delta_i(s_i)\}^2 + \frac{1}{3} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \frac{\{\delta_i(s_i)\}^3}{\{1 + \delta_i(s_i^{\star})\}^3}, \end{align}\] where \(\widetilde{\boldsymbol{\varepsilon}} = (\widetilde{\varepsilon}_1,\ldots,\widetilde{\varepsilon}_d)\). Using the last equation, the quantity \(\gamma^{\star}(\boldsymbol{s})\), as defined in 8 , can be rewritten as \[\label{eq:gamma46star46alternative} \begin{align} \gamma^{\star}(\boldsymbol{s}) &\equiv \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \left(\frac{s_i}{p_i}\right) - \left\{\widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{p})^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p})\right\} \\ &= - \frac{1}{2} \sum_{i=1}^{d+1} \widetilde{\varepsilon}_i \{\delta_i(s_i)\}^2 + \frac{1}{3} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \frac{\{\delta_i(s_i)\}^3}{\{1 + \delta_i(s_i^{\star})\}^3}. \end{align}\tag{13}\] It readily follows that \[\label{eq:entropy462} \begin{align} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \left(\frac{s_i}{p_i}\right) &= \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{p})^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) + \gamma^{\star}(\boldsymbol{s}) \\ &= \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{J}/N)^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{J}/N) + \gamma^{\star}(\boldsymbol{s}). \end{align}\tag{14}\] Therefore, by putting 10 and 14 back into 9 , and exponentiating, one gets \[\begin{align} \frac{(N + d)!}{N!} \times \frac{N!}{\prod_{i=1}^{d+1} J_i!} \prod_{i=1}^{d+1} s_i^{J_i} &= \exp\left\{\Delta_N + N \gamma^{\star}(\boldsymbol{s})\right\} \frac{\exp\big\{-N (\boldsymbol{s} - \boldsymbol{J}/N)^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{J}/N) / 2\big\}}{\sqrt{(2\pi / N)^d \det(\Sigma_{\boldsymbol{p}})}} \\ &= \exp\left\{\Delta_N + N \gamma^{\star}(\boldsymbol{s})\right\} N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}}\{N^{1/2} (\boldsymbol{J}/N - \boldsymbol{s})\} \\[2mm] &= \exp\left\{\Delta_N + N \gamma^{\star}(\boldsymbol{s})\right\} N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}}\{N^{1/2} (\boldsymbol{p} - \boldsymbol{s} + \widetilde{\boldsymbol{\varepsilon}})\}, \end{align}\] where we recall that \(\Delta_N = \ln \{(N + d)!/(N! N^d)\} + \Lambda_N - (1/2) \sum_{i=1}^{d+1} \ln (1 + \varepsilon_i) - N \widetilde{\gamma}(\boldsymbol{\varepsilon})\) from 6 . Putting this last equation in 4 yields \[\mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) = e^{\Delta_N} \int_{\mathcal{R}_d^{\circ}} \exp\left[N \gamma^{\star}(\boldsymbol{s})\right] N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}}\{N^{1/2} (\boldsymbol{p} - \boldsymbol{s} + \widetilde{\boldsymbol{\varepsilon}})\} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}.\] This concludes the proof. ◻

Second proof of Theorem 1. For all \(\boldsymbol{s}\in \mathrm{Int}(\mathcal{S}_d)\), with \(s_{d+1} = 1 - \|\boldsymbol{s}\|_1\), define \[H(\boldsymbol{s}) = \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln s_i = \sum_{i=1}^{d+1} \frac{J_i}{N} \ln s_i.\] For all \(\boldsymbol{k}\in \mathcal{K}_{n,d}\), the integral representation of the joint survival function of the cumulative sums of the components of any multinomial random vector derived in 4 can be rewritten as \[\label{eq:normal46integral46representation46eq46begin} \mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) = \frac{(N + d)!}{N!} \times \frac{N!}{\prod_{i=1}^{d+1} J_i!} \int_{\mathcal{R}_d^{\circ}} \exp\left\{N H(\boldsymbol{s})\right\} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}.\tag{15}\] By Stirling’s formula 5 , note that \[m! = \sqrt{2\pi m} \, \exp(m \ln m - m + \lambda_m),\] where \((12 m + 1)^{-1} \leq \lambda_m \leq (12 m)^{-1}\). Since \(\sum_{i=1}^{d+1} J_i = N\) and \(\Lambda_N = \lambda_N - \sum_{i=1}^{d+1} \lambda_{J_i}\), one deduces that \[\frac{N!}{\prod_{i=1}^{d+1} J_i!} = \frac{\exp(\Lambda_N)}{\sqrt{(2\pi N)^d \prod_{i=1}^{d+1} \{p_i (1 + \varepsilon_i)\}}} \exp\left\{- N H(\boldsymbol{J}/N)\right\}.\] Therefore, one can rewrite 15 as \[\label{eq:normal46integral46representation46eq46next} \begin{align} &\mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) \\ &\qquad= \frac{(N + d)!}{N! N^d} \times \frac{N^{d/2} \exp(\Lambda_N)}{\sqrt{(2\pi)^d \prod_{i=1}^{d+1} \{p_i (1 + \varepsilon_i)\}}} \int_{\mathcal{R}_d^{\circ}} \exp\left[N \left\{H(\boldsymbol{s}) - H(\boldsymbol{J}/N)\right\}\right] \, \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}. \end{align}\tag{16}\] Decompose the integral above as \[\label{eq:decomposition} \int_{\mathcal{R}_d^{\circ}} \exp\left[N \left\{H(\boldsymbol{s}) - H(\boldsymbol{p})\right\} + N \left\{H(\boldsymbol{p}) - H(\boldsymbol{J}/N)\right\}\right] \, \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}.\tag{17}\] Using the quantity \(\widetilde{\gamma}(\boldsymbol{\varepsilon})\) as defined in 7 , one has \[\label{eq:max46diff46H} H(\boldsymbol{p}) - H(\boldsymbol{J}/N) = - \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln(1 + \varepsilon_i) = - \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} - \widetilde{\gamma}(\boldsymbol{\varepsilon}),\tag{18}\] where we recall that \(\widetilde{\boldsymbol{\varepsilon}} = (p_1 \varepsilon_1,\ldots, p_d \varepsilon_d)^{\top}\), \(\widetilde{\varepsilon}_{d+1} = p_{d+1} \varepsilon_{d+1}\), and \(\Sigma_{\boldsymbol{p}}^{-1} = \{p_i^{-1} \mathbb{1}\{i = j\} + p_{d+1}^{-1}\}_{1 \leq i,j \leq d}\). Recall from 6 that \[\Delta_N = \ln\left\{\frac{(N + d)!}{N! N^d}\right\} + \Lambda_N - \frac{1}{2} \sum_{i=1}^{d+1} \ln(1 + \varepsilon_i) - N \widetilde{\gamma}(\boldsymbol{\varepsilon}),\] so one can rewrite 16 as \[\begin{align} &\mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) \\[1mm] &\qquad= \frac{N^{d/2} \exp(\Delta_N)}{\sqrt{(2\pi)^d \det(\Sigma_{\boldsymbol{p}})}} \int_{\mathcal{R}_d^{\circ}} \exp\left[N \left\{H(\boldsymbol{s}) - H(\boldsymbol{p})\right\} - \frac{N}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}}\right] \, \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}. \end{align}\] Using 13 and 14 , one has \[\begin{align} H(\boldsymbol{s}) - H(\boldsymbol{p}) &= \sum_{i=1}^{d+1} \frac{J_i}{N} \ln (s_i) - \sum_{i=1}^{d+1} \frac{J_i}{N} \ln (p_i) = \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \ln \left(\frac{s_i}{p_i}\right) \\ &= \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{p})^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} \sum_{i=1}^{d+1} \widetilde{\varepsilon}_i \{\delta_i(s_i)\}^2 + \frac{1}{3} \sum_{i=1}^{d+1} p_i (1 + \varepsilon_i) \frac{\{\delta_i(s_i)\}^3}{\{1 + \delta_i(s_i^{\star})\}^3} \\ &= \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{p})^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{p}) + \gamma^{\star}(\boldsymbol{s}) \\ &= \frac{1}{2} \widetilde{\boldsymbol{\varepsilon}}^{\top} \Sigma_{\boldsymbol{p}}^{-1} \widetilde{\boldsymbol{\varepsilon}} - \frac{1}{2} (\boldsymbol{s} - \boldsymbol{J}/N)^{\top} \Sigma_{\boldsymbol{p}}^{-1} (\boldsymbol{s} - \boldsymbol{J}/N) + \gamma^{\star}(\boldsymbol{s}). \end{align}\] Putting this in the previous equation yields \[\begin{align} \mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]) &= e^{\Delta_N} \int_{\mathcal{R}_d^{\circ}} \exp\left\{N \gamma^{\star}(\boldsymbol{s})\right\} N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}}\{N^{1/2} (\boldsymbol{J}/N - \boldsymbol{s})\} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s} \\ &= e^{\Delta_N} \int_{\mathcal{R}_d^{\circ}} \exp\left\{N \gamma^{\star}(\boldsymbol{s})\right\} N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}}\{N^{1/2} (\boldsymbol{p} - \boldsymbol{s} + \widetilde{\boldsymbol{\varepsilon}})\} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}. \end{align}\] This concludes the proof. ◻

Remark 3. The choice of decomposition in 17 can be explained as follows. The Hessian matrix of \(H\) at \(\boldsymbol{s}\in \mathrm{Int}(\mathcal{S}_d)\) is equal to \[\left(-\frac{J_i}{N s_i^2} \mathbb{1}\{i = j\} - \frac{J_{d+1}}{N s_{d+1}^2}\right)_{1 \leq i,j \leq d},\] which is (symmetric) negative definite, so that \(H\) is strictly concave on \(\mathrm{Int}(\mathcal{S}_d)\). Given the strict concavity and the fact that \(\overline{\boldsymbol{s}} = \boldsymbol{J}/N\in \mathrm{Int}(\mathcal{S}_d)\) is the unique point satisfying \[\nabla H(\overline{\boldsymbol{s}}) = \left(\frac{J_i}{N \, \overline{s}_i} - \frac{J_{d+1}}{N \, \overline{s}_{d+1}}\right)_{1 \leq i \leq d} = \boldsymbol{0}\] in the interior of the simplex, the global maximum of \(H\) on \(\mathrm{Int}(\mathcal{S}_d)\) is achieved at \(\boldsymbol{J}/N\). If \(\boldsymbol{J}/N\in \mathcal{R}_d^{\circ}\), then 17 is the usual centering at the unconstrained maximum in Laplace’s method. If \(\boldsymbol{J}/N\notin \mathcal{R}_d^{\circ}\), then the constrained maximum over \(\mathcal{R}_d\) is attained on the boundary of \(\mathcal{R}_d\) and need not be equal to \(\boldsymbol{p}\); in that case, 17 should be viewed simply as an exact algebraic decomposition.

5 Numerical checks↩︎

We numerically check the exact identity in Theorem 1 by comparing, for randomly generated admissible triples \((n,\boldsymbol{p},\boldsymbol{k})\), the two quantities \[\begin{align} L(n,\boldsymbol{p},\boldsymbol{k}) &= \mathsf{P}(X_1 + \dots + X_i \geq k_1 + \dots + k_i, ~\forall i\in [d]), \\ R(n,\boldsymbol{p},\boldsymbol{k}) &= e^{\Delta_N} \int_{\mathcal{R}_d^{\circ}} \exp\left\{N\gamma^{\star}(\boldsymbol{s})\right\} N^{d/2} \phi_{\Sigma_{\boldsymbol{p}}} \{N^{1/2}(\boldsymbol{p}-\boldsymbol{s}+\widetilde{\boldsymbol{\varepsilon}})\} \hskip 1pt\mathrm{d}\hskip 0.75pt\boldsymbol{s}. \end{align}\] The first quantity, \(L(n,\boldsymbol{p},\boldsymbol{k})\), is evaluated by direct summation of the multinomial probabilities, \[L(n,\boldsymbol{p},\boldsymbol{k}) = \sum_{\substack{x_1,\ldots,x_{d+1}\in\mathbb{N}_0 \\ x_1+\cdots+x_{d+1} = n}} \frac{n!}{x_1!\cdots x_{d+1}!} \prod_{j=1}^{d+1}p_j^{x_j} \mathbb{1}\{x_1+\cdots+x_i \geq k_1+\cdots+k_i,~\forall i\in[d]\},\] where \(p_{d+1} = 1-\|\boldsymbol{p}\|_1\). This is an exact finite sum, up to floating-point roundoff; no Monte Carlo approximation is used.

The second quantity is evaluated by deterministic adaptive numerical integration. To integrate over \(\mathcal{R}_d\), the code uses cumulative coordinates. Let \(P_i = p_1+\cdots+p_i\) for \(i\in[d]\), and map \(\boldsymbol{u}\in[0,1]^d\) to \(\boldsymbol{s}\in\mathcal{R}_d\) recursively by \[q_0 = 0,\qquad q_i = q_{i-1}+u_i(P_i-q_{i-1}),\qquad s_i = q_i-q_{i-1},\qquad i\in[d].\] The Jacobian of this transformation is \(\prod_{i=1}^d (P_i-q_{i-1})\), and the resulting integral over \([0,1]^d\) is computed using adaptive cubature. The comparison therefore checks the equality in Theorem 1 by comparing an exact multinomial summation against a deterministic numerical integral.

For each tested case, we report the absolute and symmetric relative errors: \[\mathrm{AbsErr} = |R(n,\boldsymbol{p},\boldsymbol{k})-L(n,\boldsymbol{p},\boldsymbol{k})|, \qquad \mathrm{RelErr} = \frac{2|R(n,\boldsymbol{p},\boldsymbol{k})-L(n,\boldsymbol{p},\boldsymbol{k})|}{|R(n,\boldsymbol{p},\boldsymbol{k})|+|L(n,\boldsymbol{p},\boldsymbol{k})|}.\] The symmetric relative error is used to avoid favoring either side of the identity.

The numerical experiment was run for \(d = 2,3,4,5\); see Tables 14. For each dimension, \(40\) randomized cases were generated. In each case, \(n\) was sampled uniformly from a finite integer interval, and the full probability vector \[\boldsymbol{p}^{+} = (p_1,\ldots,p_d,p_{d+1})\] was sampled from a symmetric Dirichlet distribution with parameter \(2\), conditional on all components being at least \(0.03\). The vector \(\boldsymbol{k}\) was sampled subject to the theorem’s restrictions \(k_i\geq 2\) for all \(i\in[d]\) and \(\|\boldsymbol{k}\|_1\leq n-1\). More precisely, after assigning the baseline value \(2\) to each component of \(\boldsymbol{k}\), the remaining slack \(n-1-2d\) was randomly allocated among \(d+1\) cells using probabilities sampled from a symmetric Dirichlet distribution with parameter \(1\), and only the first \(d\) allocations were added to \(\boldsymbol{k}\). Thus every generated case satisfies \(\boldsymbol{k}\in\mathcal{K}_{n,d}\). The ranges for \(n\) and the cubature tolerances were set to \(12\leq n\leq 35\) and \(10^{-7}\), respectively, for all \(d = 2,3,4,5\).

The R code that generated Tables 14 is available online in the GitHub repository [17].

Table 1: Results from 40 randomized checks of the main theorem for \(d=2\). Here \(L\) denotes the left side and \(R\) denotes the right side of the equation in Theorem 1, and the relative error is \(2|R-L|/(|R|+|L|)\).
Case \(n\) \(\boldsymbol{p}^+\) \(\boldsymbol{k}\) Abs. error Rel. error
1 33 \((0.277, 0.209, 0.513)\) \((7, 13)\) 8.381e-12 7.330e-11
2 14 \((0.149, 0.575, 0.276)\) \((3, 6)\) 6.537e-12 2.063e-11
3 16 \((0.249, 0.221, 0.530)\) \((3, 8)\) 7.723e-12 1.173e-10
4 28 \((0.447, 0.240, 0.313)\) \((13, 4)\) 2.864e-11 5.854e-11
5 25 \((0.184, 0.664, 0.153)\) \((9, 10)\) 6.767e-12 2.426e-10
6 33 \((0.071, 0.478, 0.451)\) \((22, 7)\) 2.177e-19 1.961e-03
7 33 \((0.156, 0.614, 0.230)\) \((4, 21)\) 5.533e-10 1.023e-09
8 17 \((0.200, 0.414, 0.386)\) \((2, 5)\) 2.721e-11 3.140e-11
9 14 \((0.301, 0.293, 0.406)\) \((2, 6)\) 5.970e-10 8.992e-10
10 35 \((0.310, 0.266, 0.424)\) \((15, 12)\) 1.840e-12 2.551e-10
11 13 \((0.624, 0.159, 0.218)\) \((4, 2)\) 2.399e-11 2.414e-11
12 31 \((0.524, 0.083, 0.393)\) \((18, 10)\) 3.518e-13 1.287e-09
13 25 \((0.120, 0.451, 0.428)\) \((2, 3)\) 1.060e-10 1.294e-10
14 32 \((0.061, 0.741, 0.198)\) \((13, 11)\) 8.192e-18 4.330e-10
15 31 \((0.108, 0.343, 0.549)\) \((4, 9)\) 2.185e-10 6.059e-10
16 25 \((0.126, 0.340, 0.535)\) \((4, 14)\) 3.205e-11 4.621e-09
17 31 \((0.479, 0.127, 0.394)\) \((12, 11)\) 4.244e-10 5.094e-09
18 12 \((0.413, 0.249, 0.338)\) \((4, 5)\) 1.835e-12 5.120e-12
19 30 \((0.082, 0.467, 0.451)\) \((16, 9)\) 5.537e-16 1.097e-05
20 19 \((0.119, 0.753, 0.128)\) \((2, 15)\) 1.032e-10 2.596e-10
21 23 \((0.370, 0.303, 0.326)\) \((4, 10)\) 5.045e-10 6.227e-10
22 21 \((0.208, 0.490, 0.301)\) \((14, 6)\) 8.603e-17 5.945e-11
23 28 \((0.622, 0.095, 0.283)\) \((11, 3)\) 1.803e-12 1.815e-12
24 19 \((0.162, 0.596, 0.242)\) \((7, 5)\) 1.020e-11 4.208e-10
25 32 \((0.185, 0.386, 0.430)\) \((11, 8)\) 6.337e-11 2.984e-09
26 14 \((0.399, 0.172, 0.428)\) \((2, 6)\) 7.570e-12 1.238e-11
27 29 \((0.348, 0.340, 0.312)\) \((10, 18)\) 1.324e-13 5.022e-10
28 12 \((0.435, 0.427, 0.138)\) \((4, 2)\) 6.192e-11 7.354e-11
29 33 \((0.274, 0.558, 0.168)\) \((3, 22)\) 3.320e-10 3.658e-10
30 18 \((0.638, 0.153, 0.209)\) \((3, 13)\) 1.484e-10 6.111e-10
31 15 \((0.230, 0.357, 0.413)\) \((3, 6)\) 2.005e-11 4.305e-11
32 18 \((0.143, 0.485, 0.372)\) \((8, 8)\) 1.969e-14 7.514e-11
33 28 \((0.312, 0.524, 0.165)\) \((22, 4)\) 3.349e-17 1.200e-10
34 35 \((0.244, 0.334, 0.423)\) \((12, 15)\) 1.212e-11 1.734e-09
35 19 \((0.170, 0.455, 0.375)\) \((11, 2)\) 5.249e-15 7.843e-11
36 29 \((0.158, 0.605, 0.237)\) \((12, 6)\) 1.258e-12 1.412e-09
37 14 \((0.374, 0.268, 0.358)\) \((3, 8)\) 5.201e-11 2.590e-10
38 23 \((0.284, 0.244, 0.472)\) \((3, 3)\) 1.473e-10 1.510e-10
39 12 \((0.123, 0.500, 0.377)\) \((3, 5)\) 1.303e-11 1.042e-10
40 25 \((0.386, 0.423, 0.191)\) \((2, 13)\) 1.956e-09 1.964e-09

Errors larger than \(10^{-6}\)

are highlighted in bold.

Table 2: Results from 40 randomized checks of the main theorem for \(d=3\). Here \(L\) denotes the left side and \(R\) denotes the right side of the equation in Theorem 1, and the relative error is \(2|R-L|/(|R|+|L|)\).
Case \(n\) \(\boldsymbol{p}^+\) \(\boldsymbol{k}\) Abs. error Rel. error
1 35 \((0.145, 0.071, 0.242, 0.542)\) \((2, 15, 12)\) 6.359e-16 5.696e-10
2 16 \((0.334, 0.173, 0.383, 0.110)\) \((2, 3, 2)\) 7.017e-11 7.320e-11
3 35 \((0.291, 0.185, 0.480, 0.044)\) \((15, 2, 5)\) 1.116e-11 1.957e-10
4 21 \((0.112, 0.555, 0.267, 0.065)\) \((5, 5, 7)\) 4.691e-12 6.016e-11
5 32 \((0.341, 0.151, 0.171, 0.337)\) \((2, 12, 13)\) 1.802e-11 9.159e-10
6 30 \((0.309, 0.266, 0.067, 0.357)\) \((19, 5, 5)\) 3.024e-15 2.121e-09
7 25 \((0.041, 0.470, 0.144, 0.345)\) \((7, 7, 4)\) 1.099e-14 3.342e-10
8 13 \((0.228, 0.255, 0.277, 0.240)\) \((5, 4, 2)\) 2.516e-12 5.962e-11
9 28 \((0.080, 0.278, 0.428, 0.214)\) \((3, 11, 4)\) 1.478e-11 2.506e-10
10 21 \((0.082, 0.123, 0.693, 0.103)\) \((3, 6, 2)\) 3.335e-12 2.571e-10
11 22 \((0.179, 0.073, 0.195, 0.553)\) \((2, 9, 5)\) 1.474e-13 7.438e-11
12 23 \((0.432, 0.267, 0.070, 0.230)\) \((6, 6, 8)\) 1.693e-10 8.968e-10
13 34 \((0.351, 0.366, 0.156, 0.127)\) \((9, 4, 9)\) 1.335e-10 1.494e-10
14 18 \((0.077, 0.437, 0.324, 0.163)\) \((4, 5, 2)\) 6.677e-13 1.671e-11
15 33 \((0.214, 0.503, 0.107, 0.176)\) \((8, 2, 8)\) 5.313e-12 1.291e-11
16 18 \((0.293, 0.302, 0.184, 0.221)\) \((4, 3, 5)\) 2.904e-11 3.783e-11
17 17 \((0.250, 0.208, 0.191, 0.351)\) \((4, 5, 3)\) 4.307e-11 1.822e-10
18 18 \((0.305, 0.064, 0.107, 0.525)\) \((5, 7, 3)\) 1.099e-12 9.402e-10
19 14 \((0.198, 0.189, 0.274, 0.339)\) \((3, 6, 3)\) 1.789e-11 8.474e-10
20 29 \((0.194, 0.364, 0.309, 0.133)\) \((15, 4, 7)\) 6.473e-15 8.165e-11
21 31 \((0.115, 0.146, 0.446, 0.294)\) \((7, 6, 7)\) 3.968e-13 2.582e-11
22 20 \((0.579, 0.200, 0.053, 0.169)\) \((7, 4, 4)\) 2.662e-10 2.996e-10
23 23 \((0.545, 0.204, 0.144, 0.107)\) \((7, 5, 2)\) 3.338e-10 3.369e-10
24 31 \((0.258, 0.461, 0.151, 0.130)\) \((5, 5, 13)\) 6.744e-11 7.339e-11
25 17 \((0.282, 0.167, 0.181, 0.370)\) \((10, 2, 3)\) 2.302e-13 1.425e-10
26 34 \((0.128, 0.221, 0.188, 0.463)\) \((10, 3, 6)\) 1.840e-12 2.468e-10
27 25 \((0.238, 0.375, 0.297, 0.089)\) \((8, 6, 7)\) 2.777e-11 1.320e-10
28 15 \((0.122, 0.443, 0.211, 0.224)\) \((3, 8, 3)\) 5.945e-12 2.187e-10
29 12 \((0.101, 0.237, 0.433, 0.228)\) \((4, 3, 2)\) 1.407e-12 1.369e-10
30 19 \((0.388, 0.173, 0.193, 0.246)\) \((5, 2, 7)\) 2.510e-10 3.858e-10
31 22 \((0.089, 0.113, 0.360, 0.438)\) \((5, 9, 6)\) 4.910e-17 5.126e-11
32 29 \((0.115, 0.365, 0.371, 0.149)\) \((8, 9, 4)\) 2.164e-12 2.750e-10
33 14 \((0.664, 0.102, 0.084, 0.150)\) \((3, 6, 2)\) 6.641e-11 8.021e-11
34 35 \((0.351, 0.111, 0.314, 0.223)\) \((4, 10, 3)\) 2.316e-11 2.829e-11
35 22 \((0.330, 0.215, 0.300, 0.156)\) \((2, 4, 12)\) 1.628e-12 2.177e-12
36 20 \((0.320, 0.252, 0.333, 0.095)\) \((11, 3, 4)\) 9.725e-13 4.943e-11
37 16 \((0.257, 0.244, 0.214, 0.285)\) \((5, 6, 3)\) 1.994e-11 4.932e-10
38 35 \((0.294, 0.344, 0.323, 0.039)\) \((5, 11, 11)\) 2.180e-11 2.223e-11
39 25 \((0.480, 0.269, 0.057, 0.195)\) \((9, 8, 6)\) 8.902e-11 8.243e-10
40 13 \((0.420, 0.194, 0.320, 0.066)\) \((3, 4, 4)\) 7.452e-11 9.735e-11

Errors larger than \(10^{-6}\)

are highlighted in bold.

Table 3: Results from 40 randomized checks of the main theorem for \(d=4\). Here \(L\) denotes the left side and \(R\) denotes the right side of the equation in Theorem 1, and the relative error is \(2|R-L|/(|R|+|L|)\).
Case \(n\) \(\boldsymbol{p}^+\) \(\boldsymbol{k}\) Abs. error Rel. error
1 25 \((0.084, 0.184, 0.259, 0.406, 0.067)\) \((4, 6, 3, 10)\) 1.027e-12 2.539e-11
2 33 \((0.079, 0.323, 0.277, 0.186, 0.135)\) \((2, 15, 13, 2)\) 5.760e-14 6.661e-11
3 17 \((0.260, 0.191, 0.046, 0.292, 0.211)\) \((2, 5, 6, 3)\) 8.684e-13 7.474e-11
4 18 \((0.388, 0.248, 0.227, 0.054, 0.083)\) \((5, 3, 2, 3)\) 8.157e-10 9.312e-10
5 23 \((0.179, 0.233, 0.159, 0.297, 0.132)\) \((2, 2, 12, 5)\) 4.141e-12 3.765e-11
6 25 \((0.267, 0.074, 0.112, 0.280, 0.267)\) \((3, 3, 2, 14)\) 6.534e-11 9.657e-10
7 12 \((0.135, 0.127, 0.417, 0.057, 0.263)\) \((2, 2, 2, 4)\) 1.283e-11 8.436e-11
8 14 \((0.209, 0.137, 0.087, 0.476, 0.091)\) \((3, 2, 6, 2)\) 4.103e-13 5.670e-11
9 16 \((0.213, 0.367, 0.140, 0.211, 0.069)\) \((2, 3, 2, 2)\) 1.413e-10 1.608e-10
10 29 \((0.080, 0.260, 0.407, 0.223, 0.031)\) \((3, 9, 2, 6)\) 1.385e-11 8.685e-11
11 34 \((0.498, 0.196, 0.176, 0.079, 0.051)\) \((6, 5, 11, 3)\) 6.169e-09 6.170e-09
12 27 \((0.286, 0.235, 0.164, 0.198, 0.117)\) \((3, 7, 13, 3)\) 1.911e-12 7.615e-11
13 27 \((0.153, 0.186, 0.142, 0.139, 0.380)\) \((9, 7, 8, 2)\) 2.024e-16 2.726e-10
14 18 \((0.237, 0.193, 0.155, 0.209, 0.206)\) \((2, 2, 4, 2)\) 2.090e-10 2.363e-10
15 18 \((0.390, 0.188, 0.052, 0.207, 0.162)\) \((4, 3, 4, 6)\) 9.655e-12 5.789e-11
16 31 \((0.277, 0.045, 0.287, 0.316, 0.076)\) \((7, 3, 7, 10)\) 1.613e-10 3.290e-10
17 17 \((0.081, 0.079, 0.437, 0.128, 0.274)\) \((6, 2, 5, 3)\) 5.173e-15 6.530e-11
18 26 \((0.234, 0.285, 0.042, 0.279, 0.159)\) \((3, 4, 8, 2)\) 2.816e-10 5.472e-10
19 20 \((0.260, 0.291, 0.151, 0.107, 0.191)\) \((2, 2, 5, 6)\) 3.691e-10 4.500e-10
20 18 \((0.213, 0.041, 0.201, 0.227, 0.318)\) \((3, 2, 6, 3)\) 3.540e-12 4.186e-11
21 22 \((0.237, 0.339, 0.110, 0.246, 0.069)\) \((4, 8, 3, 3)\) 7.543e-12 1.545e-11
22 16 \((0.139, 0.227, 0.127, 0.187, 0.319)\) \((7, 3, 2, 3)\) 1.866e-14 7.125e-11
23 28 \((0.047, 0.188, 0.374, 0.321, 0.071)\) \((3, 2, 10, 11)\) 6.638e-12 6.884e-11
24 16 \((0.199, 0.236, 0.260, 0.192, 0.113)\) \((2, 4, 2, 6)\) 1.233e-11 2.237e-11
25 21 \((0.242, 0.110, 0.210, 0.385, 0.053)\) \((5, 3, 2, 9)\) 4.149e-11 1.092e-10
26 19 \((0.052, 0.110, 0.330, 0.238, 0.269)\) \((3, 9, 4, 2)\) 1.244e-16 1.745e-10
27 18 \((0.052, 0.232, 0.238, 0.339, 0.140)\) \((3, 6, 5, 2)\) 7.462e-14 3.526e-11
28 30 \((0.142, 0.149, 0.311, 0.144, 0.254)\) \((7, 3, 2, 17)\) 6.329e-14 1.229e-10
29 23 \((0.089, 0.276, 0.173, 0.290, 0.173)\) \((2, 6, 4, 6)\) 6.762e-12 1.935e-11
30 25 \((0.216, 0.128, 0.146, 0.166, 0.344)\) \((3, 4, 8, 2)\) 5.779e-11 3.415e-10
31 32 \((0.387, 0.031, 0.239, 0.064, 0.279)\) \((2, 5, 11, 13)\) 3.406e-13 8.799e-10
32 31 \((0.292, 0.158, 0.254, 0.256, 0.040)\) \((3, 4, 8, 9)\) 3.664e-09 3.690e-09
33 27 \((0.105, 0.146, 0.401, 0.158, 0.189)\) \((3, 9, 3, 4)\) 5.993e-12 2.849e-10
34 29 \((0.241, 0.237, 0.279, 0.164, 0.079)\) \((9, 7, 9, 2)\) 6.904e-12 1.598e-10
35 14 \((0.161, 0.213, 0.388, 0.159, 0.079)\) \((3, 2, 2, 4)\) 2.260e-12 6.552e-12
36 34 \((0.264, 0.100, 0.468, 0.075, 0.094)\) \((2, 4, 4, 14)\) 1.962e-09 1.972e-09
37 13 \((0.290, 0.232, 0.047, 0.267, 0.164)\) \((2, 3, 3, 4)\) 6.687e-12 2.799e-11
38 30 \((0.180, 0.468, 0.173, 0.137, 0.041)\) \((4, 7, 10, 2)\) 5.290e-10 6.651e-10
39 18 \((0.044, 0.405, 0.108, 0.333, 0.110)\) \((5, 6, 4, 2)\) 7.570e-15 1.317e-10
40 18 \((0.327, 0.057, 0.088, 0.386, 0.141)\) \((7, 5, 3, 2)\) 9.330e-14 8.192e-11

Errors larger than \(10^{-6}\)

are highlighted in bold.

Table 4: Results from 40 randomized checks of the main theorem for \(d=5\). Here \(L\) denotes the left side and \(R\) denotes the right side of the equation in Theorem 1, and the relative error is \(2|R-L|/(|R|+|L|)\).
Case \(n\) \(\boldsymbol{p}^+\) \(\boldsymbol{k}\) Abs. error Rel. error
1 34 \((0.195, 0.214, 0.239, 0.033, 0.106, 0.213)\) \((7, 3, 8, 8, 2)\) 3.309e-09 2.740e-08
2 25 \((0.181, 0.265, 0.081, 0.191, 0.203, 0.079)\) \((2, 2, 2, 6, 6)\) 1.078e-07 1.131e-07
3 22 \((0.298, 0.188, 0.160, 0.095, 0.071, 0.188)\) \((9, 3, 2, 2, 3)\) 5.178e-10 5.672e-09
4 14 \((0.078, 0.059, 0.261, 0.395, 0.139, 0.069)\) \((2, 2, 3, 3, 2)\) 7.270e-10 1.192e-08
5 28 \((0.132, 0.154, 0.278, 0.068, 0.137, 0.230)\) \((2, 6, 13, 2, 3)\) 3.399e-11 5.701e-09
6 17 \((0.044, 0.189, 0.136, 0.180, 0.248, 0.203)\) \((3, 2, 3, 2, 6)\) 9.442e-11 2.029e-08
7 28 \((0.268, 0.129, 0.068, 0.048, 0.364, 0.123)\) \((7, 10, 3, 5, 2)\) 7.993e-13 4.436e-08
8 30 \((0.215, 0.150, 0.257, 0.140, 0.156, 0.082)\) \((2, 4, 13, 4, 3)\) 3.431e-08 8.052e-08
9 13 \((0.058, 0.081, 0.185, 0.164, 0.301, 0.211)\) \((2, 3, 2, 2, 2)\) 1.797e-11 2.880e-09
10 32 \((0.216, 0.065, 0.134, 0.079, 0.169, 0.338)\) \((12, 8, 3, 3, 4)\) 1.698e-15 1.620e-09
11 27 \((0.065, 0.259, 0.286, 0.067, 0.135, 0.189)\) \((4, 5, 5, 3, 8)\) 2.136e-10 1.803e-08
12 20 \((0.298, 0.123, 0.258, 0.110, 0.111, 0.100)\) \((3, 3, 3, 2, 5)\) 8.123e-09 9.418e-09
13 22 \((0.140, 0.100, 0.130, 0.432, 0.129, 0.069)\) \((2, 7, 3, 4, 5)\) 3.396e-11 1.494e-09
14 16 \((0.324, 0.227, 0.212, 0.034, 0.168, 0.036)\) \((2, 3, 5, 2, 2)\) 4.567e-08 5.939e-08
15 22 \((0.038, 0.391, 0.151, 0.067, 0.043, 0.310)\) \((3, 4, 2, 7, 3)\) 1.917e-11 4.032e-09
16 18 \((0.332, 0.110, 0.179, 0.116, 0.216, 0.046)\) \((3, 3, 3, 2, 3)\) 2.900e-08 3.595e-08
17 15 \((0.073, 0.143, 0.456, 0.051, 0.126, 0.150)\) \((2, 2, 3, 3, 4)\) 2.792e-10 3.288e-09
18 22 \((0.131, 0.180, 0.152, 0.412, 0.071, 0.053)\) \((5, 3, 3, 3, 5)\) 6.566e-10 6.533e-09
19 26 \((0.255, 0.169, 0.147, 0.219, 0.176, 0.034)\) \((12, 3, 2, 2, 6)\) 4.444e-11 4.066e-09
20 23 \((0.084, 0.109, 0.152, 0.129, 0.152, 0.375)\) \((2, 7, 6, 2, 4)\) 6.135e-14 5.597e-10
21 15 \((0.093, 0.156, 0.172, 0.048, 0.292, 0.239)\) \((2, 2, 2, 3, 3)\) 9.867e-11 9.526e-10
22 13 \((0.060, 0.224, 0.113, 0.314, 0.244, 0.044)\) \((4, 2, 2, 2, 2)\) 2.120e-11 9.828e-09
23 19 \((0.212, 0.109, 0.099, 0.112, 0.285, 0.182)\) \((7, 2, 4, 2, 2)\) 5.476e-11 1.329e-08
24 24 \((0.237, 0.108, 0.085, 0.255, 0.187, 0.128)\) \((2, 2, 2, 11, 6)\) 6.958e-09 5.101e-08
25 19 \((0.196, 0.331, 0.170, 0.196, 0.043, 0.064)\) \((2, 2, 3, 8, 3)\) 1.018e-08 1.691e-08
26 23 \((0.054, 0.140, 0.041, 0.183, 0.077, 0.505)\) \((2, 6, 4, 2, 6)\) 3.087e-13 9.298e-09
27 16 \((0.048, 0.150, 0.223, 0.295, 0.162, 0.124)\) \((2, 3, 4, 3, 3)\) 1.320e-11 5.857e-10
28 12 \((0.142, 0.095, 0.145, 0.259, 0.093, 0.266)\) \((2, 2, 3, 2, 2)\) 8.607e-11 2.671e-09
29 22 \((0.245, 0.362, 0.144, 0.124, 0.092, 0.033)\) \((2, 7, 5, 2, 4)\) 2.188e-09 2.483e-09
30 16 \((0.198, 0.155, 0.174, 0.202, 0.169, 0.102)\) \((4, 4, 3, 2, 2)\) 2.013e-10 4.250e-09
31 34 \((0.204, 0.293, 0.108, 0.152, 0.187, 0.055)\) \((13, 2, 6, 3, 4)\) 5.261e-10 4.489e-08
32 21 \((0.178, 0.238, 0.333, 0.053, 0.083, 0.115)\) \((2, 6, 5, 2, 2)\) 6.204e-10 1.012e-09
33 15 \((0.377, 0.170, 0.110, 0.104, 0.092, 0.148)\) \((2, 3, 5, 2, 2)\) 1.353e-09 5.478e-09
34 19 \((0.049, 0.271, 0.038, 0.254, 0.085, 0.302)\) \((4, 2, 4, 5, 3)\) 1.864e-12 7.881e-09
35 13 \((0.125, 0.165, 0.216, 0.079, 0.205, 0.210)\) \((2, 3, 3, 2, 2)\) 6.518e-10 1.518e-08
36 33 \((0.115, 0.095, 0.090, 0.113, 0.256, 0.331)\) \((2, 10, 2, 3, 2)\) 1.130e-10 5.570e-09
37 32 \((0.144, 0.088, 0.388, 0.146, 0.150, 0.084)\) \((6, 10, 4, 4, 5)\) 4.174e-11 5.535e-08
38 20 \((0.097, 0.063, 0.125, 0.175, 0.137, 0.403)\) \((4, 2, 2, 3, 5)\) 2.592e-11 3.506e-09
39 29 \((0.111, 0.064, 0.438, 0.037, 0.093, 0.257)\) \((2, 3, 4, 9, 6)\) 7.617e-09 5.288e-08
40 21 \((0.162, 0.212, 0.399, 0.039, 0.107, 0.081)\) \((5, 4, 4, 2, 3)\) 3.311e-09 2.001e-08

Errors larger than \(10^{-6}\)

are highlighted in bold.

Disclosure statement↩︎

The author declares no conflicts of interest.

Funding↩︎

F.Ouimet is supported by the Natural Sciences and Engineering Research Council of Canada through Discovery Grant RGPIN-2026-04471 and Discovery Launch Supplement DGECR-2026-00449.

References↩︎

References↩︎

[1]
A. V. Carter and D. Pollard, MR2154001Tusnády’s inequality revisited,” Ann. Statist., vol. 32, no. 6, pp. 2731–2741, 2004, doi: 10.1214/009053604000000733.
[2]
U. Einmahl, MR996984“Extensions of results of Komlós, Major, and Tusnády to the multivariate case,” J. Multivariate Anal., vol. 28, no. 1, pp. 20–68, 1989, doi: 10.1016/0047-259X(89)90097-3.
[3]
J. Komlós, P. Major, and G. Tusnády, MR0375412“An approximation of partial sums of independent RV’-s, and the sample DF. I,” Z. Wahrscheinlichkeitstheorie und Verw. Gebiete, vol. 32, no. 1–2, pp. 111–131, 1975, doi: 10.1007/BF00533093.
[4]
J. Komlós, P. Major, and G. Tusnády, MR0402883“An approximation of partial sums of independent RV’s, and the sample DF. II,” Z. Wahrscheinlichkeitstheorie und Verw. Gebiete, vol. 34, no. 1, pp. 33–58, 1976, doi: 10.1007/BF00532688.
[5]
M. Csörgő and P. Révész, MR0666546Strong approximations in probability and statistics. New York–London: Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], 1981, p. 284.
[6]
J. Bretagnolle and P. Massart, MR0972783“Hungarian constructions from the nonasymptotic viewpoint,” Ann. Probab., vol. 17, no. 1, pp. 239–256, 1989, doi: 10.1214/aop/1176991506.
[7]
P. Major, Manuscript/notes“The approximation of the normalized empirical distribution function by a Brownian bridge.” 2000, [Online]. Available: https://users.renyi.hu/~major/probability/empir.pdf.
[8]
G. Tusnády, In Hungarian. English title: A study of statistical hypotheses“Statisztikai hipotézisek vizsgálata,” Candidatus dissertation, Hungarian Academy of Sciences, Budapest, 1977.
[9]
P. Massart, MR1955348“Tusnady’s lemma, 24 years later,” Ann. Inst. H. Poincaré Probab. Statist., vol. 38, no. 6, pp. 991–1007, 2002, doi: 10.1016/S0246-0203(02)01130-5.
[10]
F. Ouimet, MR4340237“An improvement of Tusnády’s inequality in the bulk,” Adv. in Appl. Math., vol. 133, pp. Paper No. 102270, 24 pp., 2022, doi: 10.1016/j.aam.2021.102270.
[11]
A. Yu. Zaitsev, MR1616527“Multidimensional version of the results of Komlós, Major and Tusnády for vectors with finite exponential moments,” ESAIM Probab. Statist., vol. 2, pp. 41–108, 1998, doi: 10.1051/ps:1998103.
[12]
T. A. Severini, MR2168237Elements of Distribution Theory, vol. 17. Cambridge University Press, Cambridge, 2005, p. xii+515.
[13]
K. Tanabe and M. Sagae, MR1157720“An exact Cholesky decomposition and the generalized inverse of the variance–covariance matrix of the multinomial distribution, with applications,” J. Roy. Statist. Soc. Ser. B, vol. 54, no. 1, pp. 211–219, 1992, doi: 10.1111/j.2517-6161.1992.tb01875.x.
[14]
T. H. Gronwall, MR1502544“The gamma function in the integral calculus,” Ann. of Math. (2), vol. 20, no. 2, pp. 35–124, 1918, doi: 10.2307/1967180.
[15]
M. Siotani and Y. Fujikoshi, MR750392“Asymptotic approximations for the distributions of multinomial goodness-of-fit statistics,” Hiroshima Math. J., vol. 14, no. 1, pp. 115–124, 1984, doi: 10.32917/hmj/1206133150.
[16]
F. Ouimet, MR4249129“A precise local limit theorem for the multinomial distribution and some applications,” J. Statist. Plann. Inference, vol. 215, pp. 218–233, 2021, doi: 10.1016/j.jspi.2021.03.006.
[17]
F. Ouimet, Accessed 2026-05-28MultinomialIntegralRepresentation.” GitHub repository, 2026, [Online]. Available: https://github.com/FredericOuimetMcGill/MultinomialIntegralRepresentation.