U-statistics and random subgraph counts: Multivariate normal approximation via exchangeable pairs and embedding

TL;DR

This paper introduces an embedding method using exchangeable pairs for normal approximation, applied to U-statistics and subgraph counts in random graphs, overcoming previous linearity condition limitations.

Contribution

The authors extend the embedding approach to U-statistics and subgraph counts, enabling multivariate normal approximation without linearity constraints.

Findings

Successful application to U-statistics and subgraph counts

Improved normal approximation accuracy

Broader applicability of exchangeable pairs method

Abstract

In a recent paper by the authors, a new approach--called the "embedding method"--was introduced, which allows to make use of exchangeable pairs for normal and multivariate normal approximation with Stein's method in cases where the corresponding couplings do not satisfy a certain linearity condition. The key idea is to embed the problem into a higher dimensional space in such a way that the linearity condition is then satisfied. Here we apply the embedding to U-statistics as well as to subgraph counts in random graphs.

Equations

$$\text{mathbbm}{E}^W(W^{\prime }-W)=-\Lambda W+R$$

$$\text{mathbbm}EW=0,\text{mathbbm}EW{W}^t=\Sigma ,$$

$$\bigl{|}\mathbbm{E}h(W)-\mathbbm{E}h(\Sigma^{1/2}Z)\bigr{|}\leqslant\frac{|h|_{2}}{4}A+\frac{|h|_{3}}{12}B+\Bigl{(}|h|_{1}+{\textstyle\frac{1}{2}}d\|\Sigma\|^{1/2}|h|_{2}\Bigr{)}C,$$

$$A$$

$$B$$

$$C$$

$$\text{mathbbm}E{Y}^t\nabla f(Y)=\text{mathbbm}E{\nabla }^t\Sigma \nabla f(Y),\text{for all smooth f:mathbbmRd→mathbbmR}.$$

$$\bigl{|}\mathbbm{E}h(X)-\mathbbm{E}h(Y)\bigr{|}\leqslant\frac{1}{2}|h|_{2}\sum_{i,j=1}^{d}|\sigma_{i,j}-\sigma^{0}_{i,j}|.$$

$$h_{\delta }^{+}(x)$$

$$\text{mbox}and\tilde{h}(x,\delta )$$

$${\hat{\lambda }}^{(i)}=\sum _{m=1}^d|({\Sigma }^{-1/2}{\Lambda }^{-1}{\Sigma }^{1/2}{)}_{m,i}|$$

$$A^{\prime }$$

$$B^{\prime }$$

$$\mathcal{l}abel{9}C^{\prime}=\sum_{i}\hat{\lambda}^{(i)}\sqrt{\mathbbm{E}\Bigl{(}\sum_{k}\Sigma^{-1/2}_{i,k}{R_{k}}\Bigr{)}^{2}},$$

$$\sup_{h\in{\mathcal{H}}}|\mathbbm{E}h(W)-\mathbbm{E}h(Z)|\leqslant\gamma^{2}\Bigl{(}-D^{\prime}\log(T^{\prime})+\frac{B^{\prime}}{2\sqrt{T^{\prime}}}+C^{\prime}+a\sqrt{T^{\prime}}\Bigr{)},$$

$$T^{\prime}=\frac{1}{a^{2}}\biggl{(}D^{\prime}+\sqrt{\frac{aB^{\prime}}{2}+D^{\prime 2}}\biggr{)}^{2}\qquad\text{ and }\qquad D^{\prime}=\frac{A^{\prime}}{2}+C^{\prime}d.$$

$$\begin{split}&\mathop{\mathrm{Var}}\mathbbm{E}^{Y}\sum_{k,\ell}\Sigma^{-1/2}_{i,k}\Sigma^{-1/2}_{j,\ell}(W_{k}^{\prime}-W_{k})(W_{\ell}^{\prime}-W_{\ell})\\ &\qquad\leqslant d^{4}\|\Sigma^{-1/2}\|^{4}\sup_{k,\ell}\mathop{\mathrm{Var}}\mathbbm{E}^{W}\bigl{\{}(W_{k}^{\prime}-W_{k})(W_{\ell}^{\prime}-W_{\ell})\bigr{\}}\end{split}$$

$$A^{\prime}\leqslant d^{3}\|\Sigma^{-1/2}\|^{2}\sum_{i}\hat{\lambda}^{(i)}\sup_{k,\ell}\sqrt{\mathop{\mathrm{Var}}\mathbbm{E}^{W}\bigl{\{}(W_{k}^{\prime}-W_{k})(W_{\ell}^{\prime}-W_{\ell})\bigr{\}}};$$

$$\begin{array}{rl}\text{mathbbm}E(W^{\prime }-W)(W^{\prime }-W{)}^t& =2\text{mathbbm}EW(\Lambda W{)}^t=2\Sigma {\Lambda }^t\end{array}$$

$$\text{mathbbm}E(W^{\prime }-W)(W^{\prime }-W{)}^t=\Lambda \Sigma +\Sigma {\Lambda }^t.$$

$$F$$

$$E(J_n^{\prime }(f_n)-J_n(f_n)|W)$$

$$\psi_{k}(X_{1},\dots,X_{k}):=\mathbbm{E}\bigl{(}\psi(X_{1},\dots,X_{d})\bigm{|}X_{1},\dots,X_{k}\bigr{)}.$$

$$U_k:=\sum _{|\alpha |=k}{\psi }_k(\alpha ),$$

$$W_k:=n^{1/2}{\left(\genfrac{}{}{0ex}{}{n}{k}\right)}^{-1}{U}_k.$$

$$|\text{mathbbm}Eh(W)-\text{mathbbm}Eh({\Sigma }^{1/2}Z)|⩽n^{-\frac{1}{2}}\left(4{\rho }^{1/2}{d}^6|h{|}_2+{\rho }^{3/4}{d}^7|h{|}_3\right).$$

$$\begin{split}\mathbbm{E}^{X}(U^{\prime}_{k}-U_{k})&=\frac{1}{n}\sum_{j=1}^{n}\sum_{|\alpha|=k,\@@atop\alpha\ni j}\mathbbm{E}^{X}\bigl{\{}\psi^{\prime}_{j,k}(\alpha)-\psi_{k}(\alpha)\bigr{\}}\\ &=-\frac{k}{n}U_{k}+\frac{1}{n}\sum_{j=1}^{n}\sum_{|\alpha|=k,\@@atop\alpha\ni j}\psi_{k-1}(\alpha\setminus\{j\})\\ &=-\frac{k}{n}U_{k}+\frac{n-k+1}{n}\sum_{|\beta|=k-1}\psi_{k-1}(\beta)\\ &=-\frac{k}{n}U_{k}+\frac{n-k+1}{n}U_{k-1},\end{split}$$

$$\sum _{|\alpha |=k,\text{@@atop}\alpha \ni j}{\psi }_{k-1}(\alpha \setminus \{j\})=\sum _{|\beta |=k-1,\text{@@atop}\beta \not\ni j}{\psi }_{k-1}(\beta ),$$

$$\text{mathbbm}{E}^X(W_k^{\prime }-W_k)=-\frac{k}{n}(W_k-W_{k-1}).$$

$$\Lambda =\frac{1}{n}\left[\begin{array}{ccccc}1& & & & \\ -2& 2& & & \\ & -3& 3& & \\ & & \ddots & \ddots & \\ & & & -d& d\end{array}\right],$$

$$({\Lambda }^{-1}{)}_{k,l}=n/l,$$

$${\lambda }^{(l)}⩽dn.$$

Full text

University of Oxford and National University of Singapore

1. Introduction

Stein’s method, first introduced in the 70s [11], has proven a powerful tool for assessing distributional distances, such as to the normal distribution, in the presence of dependence. When considering sums $W$ of random variables, the dependence between these random variables needs to be weak in order for the distance to a normal distribution to be small. For quantifying weak dependence, Stein [12] introduced the method of exchangeable pairs: construct a sum $W^{\prime}$ such that $(W,W^{\prime})$ form an exchangeable pair, and such that $\mathbbm{E}^{W}(W^{\prime}-W)$ is (at least approximately) linear in $W$. This linearity condition arises naturally when thinking of correlated bivariate normals. The generalisation of this approach to a multivariate setting remained untackled until recently Chatterjee and Meckes [2] solved the problem in the case of exchangeable vectors $(W,W^{\prime})$ such that $\mathbbm{E}^{W}(W^{\prime}-W)=-\lambda W+R$ for a scalar $\lambda$ and a remainder vector $R$ such that $\mathbbm{E}|R|$ is small. This is a rather special case; the authors [8] tackled the general setting that

[TABLE]

for a matrix $\Lambda$ and a vector $R$ with small $\mathbbm{E}|R|$ is treated. In in a followup paper by Meckes [5] the results by Chatterjee and Meckes [2] and by the authors [8] are combined using slightly different smoothness conditions on test functions as compared to [8]; non-smooth test functions are not treated by Meckes [5], but the bounds obtained there improve on those from [8] for the example of $d$-runs with respect to smooth test functions.

A surprising finding in [8] was that it is often possible to embed a random vector $W$ into a random vector $\hat{W}$ of larger, but still finite, dimension, such that (1.1) holds with $R=0$; yet this embedding does not correspond to Hoeffding projections. Here we explore the embedding method further, by illustrating its use on two important examples. The first example is complete non-degenerate U-statistics, and the second example considers the joint count of edges and triangles in Bernoulli random graphs. In both examples the limiting covariance matrix is not of full rank; yet the bounds on the normal approximation are of the expected order.

The paper is organised as follows. In Section 2 we review the theoretical results in [8], giving bounds on the distance to normal under the linearity condition (1.1), both for smooth test functions and for non-smooth test functions. In Section 3 we discuss the embedding method, and point out a link to Rademacher integrals and chaos decompositions. Section 4 illustrates the embedding method for complete non-degenerate U-statistics; the embedding vector contains lower-order U-statistics which are obtained via fixing components. Section 5 gives a normal approximation for the joint counts of the number of edges and the number of triangles in a Bernoulli random graph; to our knowledge these are the first explicit bounds for this multivariate problem. The embedding method suggests to count the number of 2-stars as well, which makes the results not only more informative but also, surprisingly, easier to derive.

2. Theoretical bounds for a multivariate normal approximation

2.1. Notation

Denote by $W=(W_{1},W_{2},\dots,W_{d})^{t}$ random vectors in $\mathbbm{R}^{d}$, where $W_{i}$ are $\mathbbm{R}$-values random variables for $i=1,\dots,d$. We denote by $\Sigma$ symmetric, non-negative definite matrices, and hence by $\Sigma^{1/2}$ the unique symmetric square root of $\Sigma$. Denote by $\mathop{\mathrm{Id}}$ the identity matrix, where we omit the dimension $d$. Throughout this article, $Z$ denotes a random variable having standard $d$-dimensional multivariate normal distribution. We abbreviate the transpose of the inverse of a matrix $\Lambda$ as $\Lambda^{-t}:=(\Lambda^{-1})^{t}$.

For derivatives of smooth functions $h:\mathbbm{R}^{d}\to\mathbbm{R}$, we use the notation $\nabla$ for the gradient operator. Denote by $\|\cdot\|$ the supremum norm for both functions and matrices. If the corresponding derivatives exist for some function $h:\mathbbm{R}^{d}\to\mathbbm{R}$, we abbreviate |h|_{1}:=\sup_{i}\bigl{\|}\frac{\partial}{\partial x_{i}}h\bigr{\|}, |h|_{2}:=\sup_{i,j}\bigl{\|}\frac{\partial^{2}}{\partial x_{i}\partial x_{j}}h\bigr{\|}, and so on.

We start by considering smooth test functions.

Theorem 2.1** (c.f. Theoem 2.1 [8]).**

Assume that $(W,W^{\prime})$ is an exchangeable pair of $\mathbbm{R}^{d}$-valued random variables such that

[TABLE]

with $\Sigma\in\mathbbm{R}^{d\times d}$ symmetric and positive definite. Suppose further that (1.1) is satisfied for an invertible matrix $\Lambda$ and a $\sigma(W)$-measurable random variable $R$. Then, if $Z$ has $d$-dimensional standard normal distribution, we have for every three times differentiable function $h$,

[TABLE]

where, with $\lambda^{(i)}=\sum_{m=1}^{d}|(\Lambda^{-1})_{m,i}|$,

[TABLE]

The proof of Theorem 2.1 is based on the Stein characterization of the normal distribution that ${Y}\in\mathbbm{R}^{d}$ is a multivariate normal $\mathop{\mathrm{MVN}}(0,\Sigma)$ if and only if

[TABLE]

In the paper by Meckes [5] a different norm for functions and for operators is used, to obtain a similar result, and the difference in the bounds depending on the chosen norm is illustrated for the example of runs on the line.

Theorem 2.1 can be extended to allow for covariance matrices which are not full rank, using the triangle inequality in conjunction with the following proposition.

Proposition 2.2** (c.f. Proposition 2.9 [8]).**

Let $X$ and $Y$ be $\mathbbm{R}^{d}$-valued normal variables with distributions $X\sim\mathop{\mathrm{MVN}}(0,\Sigma)$ and $Y\sim\mathop{\mathrm{MVN}}(0,\Sigma_{0})$, where $\Sigma=(\sigma_{i,j})_{i,j=1,\ldots,d}$ has full rank, and $\Sigma_{0}=(\sigma^{0}_{i,j})_{i,j=1,\ldots,d}$ is non-negative definite. Let $h:\mathbbm{R}^{d}\to\mathbbm{R}$ have 3 bounded derivatives. Then

[TABLE]

For non-smooth test functions, following Rinnot and Rotar [9], let $\Phi$ denote the standard normal distribution in $\mathbbm{R}^{d}$, and $\phi$ the corresponding density function. For $h:\mathbbm{R}^{d}\rightarrow R$ set

[TABLE]

Let ${\mathcal{H}}$ be a class of measurable functions $\mathbbm{R}^{d}\rightarrow R$ which are uniformly bounded by 1. Suppose that for any $h\in{\mathcal{H}}$, for any $\delta>0$, $h_{\delta}^{+}(x)$ and $h_{\delta}^{-}(x)$ are in ${\mathcal{H}}$; for any $d\times d$ matrix $A$ and any vector $b\in\mathbbm{R}^{d}$, $h(Ax+b)\in{\mathcal{H}}$; and for some constant $a=a({\mathcal{H}},\delta)$, $\sup_{h\in{\mathcal{H}}}\left\{\int_{\mathbbm{R}^{d}}{\tilde{h}}(x,\delta)\Phi(dx)\right\}\leqslant a\delta$. Obviously we may assume $a\geq 1$. The class of indicators of measurable convex sets is such a class where $a\leqslant 2\sqrt{d}$; see the paper by Bolthausen and Götze [1].

Let $W$ have mean vector [math] and covariance matrix $\Sigma$. If $\Lambda$ and $R$ are such that (1.1) is satisfied for $W$, then $Y=\Sigma^{-1/2}W$ satisfies (1.1) with $\hat{\Lambda}=\Sigma^{-1/2}\Lambda\Sigma^{1/2}\text{ and }R^{\prime}=\Sigma^{-1/2}R.$ With

[TABLE]

as well as

[TABLE]

and

[TABLE]

we have the following result [8].

Corollary 2.3**.**

Let $W$ be as in Theorem 2.1. Then, for all $h\in{\mathcal{H}}$ with $|h|\leqslant 1$, there exist $\gamma=\gamma(d)$ and $a>1$ such that

[TABLE]

with

[TABLE]

Remark 2.4**.**

We can simplify the above bound further. Using Minkowski’s inequality we have that $\mathop{\mathrm{Var}}\sum_{i=1}^{k}X_{i}\leqslant k^{2}\sup_{i}{\mathop{\mathrm{Var}}X_{i}},$ and thus obtain the simple estimate

[TABLE]

and hence

[TABLE]

in $B^{\prime}$ and $C^{\prime}$ we could similarly bound $\Sigma^{-1/2}_{i,k}$ by $\|\Sigma^{-1/2}\|$ to obtain a simpler bound. There are however examples, such as the random graph example in Section 5, where $\|\Sigma^{-1/2}\|$ provides a non-informative bound.

Remark 2.5**.**

Note that, if $(W,W^{\prime})$ is exchangeable and (1.1) is satisfied we have

[TABLE]

On the other hand, if we only have $\mathscr{L}(W)=\mathscr{L}(W^{\prime})$, we obtain

[TABLE]

Hence, to check in an application whether the often tedious calculation of $\Sigma$ and $\Lambda$ has been carried out correctly, we can combine Equations (2.6) and (2.7), to conclude that, under the conditions of Theorem 2.1, we must have $\Lambda\Sigma=\Sigma\Lambda^{t}$.

3. The embedding method

Assume that an $\ell$-dimensional random vector $W_{(\ell)}$ of interest is given. Often, the construction of an exchangeable pair $(W_{(\ell)},W_{(\ell)}^{\prime})$ is straightforward. If, say, $W_{(\ell)}=W_{(\ell)}(\mathbbm{X})$ is a function of i.i.d. random variables $\mathbbm{X}=(X_{1},\dots,X_{n})$, one can choose uniformly an index $I$ from $1$ to $n$, replace $X_{I}$ by an independent copy $X^{\prime}_{I}$, and define $W_{(\ell)}^{\prime}:=W_{(\ell)}(\mathbbm{X}^{\prime})$, where $\mathbbm{X}^{\prime}$ is now the vector $\mathbbm{X}$ but with $X_{I}$ replaced by $X^{\prime}_{I}$.

In general there is no hope that $(W_{(\ell)},W_{(\ell)}^{\prime})$ will satisfy Condition (1.1) with $R$ being of the required smaller order or even equal to zero, so that in this case Theorem 2.1 would not yield useful bounds.

Surprisingly often it is possible, though, to extend $W_{(\ell)}$ to a vector $W\in\mathbbm{R}^{d}$ such that we can construct an exchangeable pair $(W,W^{\prime})$ which satisfies Condition (1.1) with $R=0$. If we can bound the distance of the distribution ${\mathcal{L}}(W)$ to a $d$-dimensional multivariate normal distribution, then a bound on the distance of the distribution ${\mathcal{L}}(W_{(\ell)})$ to an $\ell-$dimentional multivariate normal distribution follows immediately.

In order to obtain useful bounds in Theorem 2.1, the embedding dimension $d$ should not be too large. In the examples below it will be obvious how to choose $W^{(d-\ell)}$ to make the construction work.

As a first illustration of the method, it was observed in [6] that for functions which depend on the first $d$ coordinates of an infinite Rademacher sequence, that is, a sequence of symmetric $\{-1,1\}$ random variables, the natural embedding vector is a vector of Rademacher integrals of lower order. A similar construction works fairly generally, as follows. Assume that $F=F(X_{1},...,X_{d})$ is a random variable that depends uniquely on the first $d$ coordinates of a sequence $X$ of i.i.d. mean zero random variables, with $E(F)=0$ and $E(F^{2})=1$, of the form

[TABLE]

such representations occur as chaotic decompositions for functionals of Rademacher sequences. A natural exchangeable pair construction is as follows. Pick an index $I$ so that $P(I=i)=\frac{1}{d}$ for $i=1,\ldots,d$, independently of $X_{1},...,X_{d}$, and if $I=i$ replace $X_{i}$ by an independent copy $X_{i}^{*}$ in all sums in the decomposition (3.1) which involve $X_{i}$. Call the resulting expression $F^{\prime}$, and the corresponding sums $J_{n}^{\prime}(f_{n});n=1,\ldots,d$. Now choosing as embedding vector ${W}=(J_{1}(f_{1}),\ldots,J_{d}(f_{d})),$ we check that for all $n=1,\ldots,d$,

[TABLE]

Thus, with ${W^{\prime}}=(J_{1}^{\prime}(f_{1}),\ldots,J_{d}^{\prime}(f_{d})),$ the condition (1.1) is satisfied, with $\Lambda=(\lambda_{i,j})_{1\leqslant i,j\leqslant d}$ being zero off the diagonal and $\lambda_{n,n}=\frac{n}{d}$ for $n=1,\ldots,d$. Note that, although diagonal, the diagonal entries of this $\Lambda$ are not equal. It is not possible to correct this by simple coordinate-wise scaling of $W$ as this will change $\Sigma$ only and leave $\Lambda$ unaffected; see also the discussion in [8, Section 5]. Hence, again, the generality of $\eqref{1}$ is essential here.

4. Complete non-degenerate $U$-statistics

Using the exchangeable pairs coupling, Rinott and Rotar [10] proved a univariate normal approximation theorem for non-degenerate weighted $U$-statistics with symmetric weight function under fairly mild conditions on the weights. Using the typical coupling, where uniformly a random variable $X_{i}$ is choosen and replaced by an independent copy, they show that (1.1) is satisfied for the one-dimensional case and a non-trivial remainder term, being Hoeffding projections of smaller order. It should not be difficult (but nevertheless cumbersome) to generalise their result to the multivariate case, where $d$ different $U$-statistics are regarded based on the same sample of independent random variables, such that (1.1) is satisfied with $\Lambda=I$ and non-trivial remainder term, again of lower order; for multivariate approximations of several $U$-statistics see also the book by Lee [4]. However, as we want to emphasize the use of Theorem 2.1 for non-diagonal $\Lambda$, we take a different approach.

Let $X_{1},\dots,X_{n}$ be a sequence of i.i.d. random elements taking values in a measure space $\mathcal{X}$. Let $\psi$ be a measurable and symmetric function from $\mathcal{X}^{d}$ to $\mathbbm{R}$, and, for each $k=1,\dots,d$, let

[TABLE]

Assume without loss of generality that $\mathbbm{E}\psi(X_{1},\dots,X_{d})=0$. For any subset $\alpha\subset\{1,\dots,n\}$ of size $k$ write $\psi_{k}(\alpha):=\psi(X_{i_{1}},\dots,X_{i_{k}})$ where the $i_{j}$ are the elements of $\alpha$. Define the statistics

[TABLE]

where $\sum_{E(\alpha)}$ denotes summation over all subsets $\alpha\subset\{1,\dots,n\}$ which satisfy the property $E$. Then $U_{d}$ coincides with the usual $U$-statistics with kernel $\psi$. Assume that $U_{d}$ is non-degenerate, that is, $\mathbbm{P}[\psi_{1}(X_{1})=0]<1$. Put

[TABLE]

It is well known that $\mathop{\mathrm{Var}}W_{k}\asymp 1$ (e.g. [4]). Note also that, as $n\to\infty$, $\Sigma:=\mathbbm{E}(WW^{t})$ will converge to a covariance matrix with all entries equal to $\mathop{\mathrm{Var}}\psi_{1}(X_{1})$ and which is thus of rank $1$, as we assume non-degeneracy and hence $U_{1}=\sum_{i=1}^{n}\psi_{1}(X_{i})$ will dominate the behaviour of each $U_{k}$.

Using Stein’s method and the approach of decomposable random variables, Raič [7] proved rates of convergence for vectors of $U$-statistics, where the coordinates are assumed to be uncorrelated (but nevertheless based upon the same sample $X_{1},\dots,X_{n}$). The next theorem can be seen as a complement to Raič’s results, as in our case from above, a normalization is not appropriate.

Theorem 4.1**.**

With the above notation, and if $\rho:=\mathbbm{E}\psi(X_{1},\dots,X_{d})^{4}<\infty$ we have for every three times differentiable function $h$

[TABLE]

Proof.

Let $X_{1}^{\prime},\dots,X_{n}^{\prime}$ be independent copies of $X_{1},\dots,X_{n}$. Define the random variables $\psi^{\prime}_{j,k}(\alpha)$ analogously to $\psi_{k}(\alpha)$ but based on the sequence $X_{1},\dots,X_{j-1},X^{\prime}_{j},X_{j+1},\dots,X_{n}$. Define the coupling as in [10], that is, pick uniformly an index $J$ from $\{1,\dots,n\}$ and replace $X_{J}$ by $X^{\prime}_{J}$, so that $U^{\prime}_{k}=\sum_{|\alpha|=k}\psi^{\prime}_{J,k}(\alpha)$; it is easy to see that $(U^{\prime},U)$ is exchangeable. Note now that, if $j\not\in\alpha$, $\psi^{\prime}_{j,k}(\alpha)=\psi_{k}(\alpha)$, and, with $X=(X_{1},\dots,X_{n})$, that $\mathbbm{E}^{X}\psi^{\prime}_{j,k}(\alpha)=\psi_{k-1}(\alpha\setminus\{j\})$ if $j\in\alpha$. Thus

[TABLE]

where the second last equality follows from the observation that

[TABLE]

and thus, in the corresponding double sum of (4.1), every set $\beta$ of size $k-1$ appears exactly $n-(k-1)$ times. Thus

[TABLE]

Hence, (1.1) is satisfied for $R=0$ and

[TABLE]

with lower triangular $\Lambda^{-1}$ such that, if $l\leqslant k$,

[TABLE]

thus, for $l=1,\dots,d$,

[TABLE]

Define now $\eta_{j,k}(\alpha):=\psi^{\prime}_{j,k}(\alpha)-\psi_{k}(\alpha)$. Then we have for every $k,l=1,\dots,d$,

[TABLE]

and

[TABLE]

Note now that, if the sets $\alpha\cup\beta$ and $\gamma\cup\delta$ are disjoint (which can only happen if $i\neq j$),

[TABLE]

due to independence. The variance of (4.3), that is (LABEL:17) minus the square of the expectation of (4.3), contains only summands where $\alpha\cup\beta$ and $\gamma\cup\delta$ are not disjoint. Recall now that $\rho=\mathbbm{E}\psi(X_{1},\dots,X_{d})^{4}$. Bounding all the non-vanishing terms simply by $32\rho$, it only remains to count the number of non-vanishing terms. Thus,

[TABLE]

where the equality is just a split of the sum over $j$ into the cases whether or not $j\in\alpha\cup\beta$. In the former case we automatically have $(\alpha\cup\beta)\cap(\gamma\cup\delta)\neq\emptyset$. It is now not difficult to see that

[TABLE]

Noting that, for fixed $j$, $k$, $l$, $\alpha$ and $\beta$,

[TABLE]

we further have

[TABLE]

where we also used that ${n-|\alpha\cup\beta|\choose k-1}\geqslant{n-k-l+1\choose k-1}$. The following statements are straightforward to prove:

[TABLE]

Thus, from (4.6),

[TABLE]

From (4.6) and (4.7),

[TABLE]

Thus, for all $k$ and $l$,

[TABLE]

Notice further that for any $m=1,\dots,d$