Consider a length $n$ vector $\mathbf{x}$ containing $n$ i.i.d. observations $\{x_i\}_{i=1}^n$ of a standard normal random variable $X$. Let $\mathbf{z}$ be a length $n-1$ vector whose entries are $z_i = x_{i+1}-x_i$. I have a statistic $s(\mathbf{z}) = \frac{1}{n-1}\sum_{i=1}^{n-1}(z_i)^2$ that is meant to measure how close to sorted a given permutation of random samples is, and I want to know its properties under the null hypothesis where the order is random, for the purpose of calculating a $p$-value.

I have two questions:

  1. Is $s(\mathbf{z})$ asymptotically normal?
  2. What is the variance of $s(\mathbf{z})$ ?

Both of these questions would be easy if the entries of $\mathbf{z}$ were i.i.d., but I have been getting stuck with the fact that they are dependent. For instance, the value of $x_{i+1}-x_i$ depends on the value of $x_{j+1}-x_j$ since $i \neq j$.

For 1., the answer seems to be "yes" based on simulations, but I'm having trouble proving this.

For 2., I thought of estimating the variance by calculating the sample variance of the squared pairwise Euclidean distances of the entries of $\mathbf{x}$, which would work if each $z_i$ were a random sample drawn from $\{x_i - x_j\}_{i,j \in \{1,\ldots,n\}}$, but this is not the case due to the dependence mentioned above. I have found empirically from simulations that dividing this pairwise distance-based estimate by 2 gives the correct result within simulation error, but this might be a lucky coincidence.

  • $\begingroup$ I notice you characterize $s$ as "a measure of how close to sorted a given permutation of random samples is." Maybe so, but it's not a good one because it depends too much on the distribution of those samples. Why aren't you using a more direct measure based on the actual permutation? For instance, if $\sigma_n$ is the permutation that would order first $n$ values, you could use some measure of correlation between $\sigma_n$ and the vector $(1,2,\ldots,n).$ $\endgroup$
    – whuber
    Commented Jun 9, 2020 at 17:47

2 Answers 2


TLDR; $s(z)$ is asymptotically normal, and its variance is $\frac {12} {n-1}$ according to CLT for Markov chains. It can be shown that the distribution is a special case of generalized $\chi^2$ distribution.

Markov Chain approach, asymptotics and variance

The sequence $z_i$ is Markov chain because once you know $z_i$ the value of $z_{i+1}$ doesn't depend on $z_k$ where $k<i$. Therefore, Markov chain CLT is applicable. Here's how we apply it.

The sum, or any linear combination of normal r.v.s, is a r.v. itself. Knowing that $z_i\sim\mathcal N(0,2)$ or $z_i\sim\sqrt 2\space \mathcal N(0,1)$, we know that $z_i^2\sim 2\space\chi^2_1$, see the definition of $\chi^2$ distribution. Thus, $\sigma_z^2=\operatorname{var}[z_i^2]=2^2\times 2=8$.

Markov chain CLT states: $$\sqrt{n-1}(s(z)-\mu)\sim\mathcal N(0,\sigma^2),$$ where $\mu=E[z_i^2]$ and $\sigma^2 = \sigma_z^2 + 2\sum_{k=1}^\infty \operatorname{cov}( z_{1}^2, z_{1+k}^2)=8+2\times 2=12$. Hence $\operatorname{var}[s(z)]=\frac{12}{n-1}$

Here's the proof by simulation (Python):

import numpy as np
n = 51
s = np.mean(np.diff(np.random.randn(10000,n))**2,axis=1)
vars = np.var(s)



Note, that if $z_i^2$ weren't correlated then $s(z)$ would have been from scaled $\chi^2$ distribution with variance $\frac 8 {n-1}$. However, due to overlapping terms $x_i$ in $z_i$ and $z_{i+1}$ we had to apply modified CLT to obtain the asymptotic distribution of $s(z)$.

Acknowledgements: My initial answer, which I updated a few times already, did not account for correlation, which was pointed out by @Sextus Empiricus. Also, I used this answer for $\operatorname{cov}( z_{1}^2, z_{1+k}^2)$, where the correlation $\rho=\operatorname{corr}[z_i,z_{i+1}]=-1/2$ and we know that correlation disappears between $z_i$ and $z_j$ when $|i-j|>1$.

The distribution

Let's start with independent row vector of randoms $X'=(x_1,\dots,x_n)$. We get the row vector of differences $Z'$ aby pplying Toeplitz matrix $B'$ as follows $Z'=X'B'$, where $$B' = \begin{bmatrix} {-1} & 0&\dots & 0 &0\\ 1 & -1 & \dots&0 & 0 \\ 0 & {1}& \dots & 0& 0\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0& 0 & \dots & 1 & -1 \\ 0& 0 & \dots & 0 & {1} \end{bmatrix}$$ Your quantity then is a quadratic form $$s(z)=\frac 1 {n-1} X'B'BX$$ where $B'B$ has a form of a tridiagonal Toeplitz matrix:

enter image description here

Let's apply eigen decomposition $B'B=P'\Lambda P$ then we have: $$s(z)=\frac 1 {n-1} X'P'\Lambda PX=\frac 1 {n-1} Y'\Lambda Y$$ where $Y=PX\sim\mathcal N(0,I_{n-1})$, i.e. each $Y_i$ (principal component) is an independent normal.

Hence, $$s(z)=\frac 1 {n-1} \sum_{i=1}^{n-1}\lambda_i Y_i^2$$ where $Y_i^2\sim\chi^2_1$ and $\lambda_i$ are eigenvalues. The eigenvalues of Toeplitz tridiagonal matrices are known to form a sine wave and easy to find, see "The Eigenproblem of a Tridiagonal P-Toeplitz Matrix" by Gover .

So the distribution can be seen as a linear combination of $\chi^2$ variables or a generalized $\chi^2$ distribution.


We can define a row vector $V'=(x_1,z_1,\dots,z_{n-1})$, then it can be obtained by applying a matrix $D'$ to the original observations $V'=X'D'$, the matrix $B'$ above is the subset of columns of $D'$:

enter image description here

The matrix $D'D$ looks like this:

enter image description here

We can get the matrix $U'$ that recovers the original vector from $V$ as follows: $X'=V'U'$, and $U'=D'^{-1}$. The matrix $U'$ is upper unit triangular, meaning $u_{ij}=1_{i\ge j}$:

enter image description here

Matrix $A=U'U$, which appears in the quadratic form, has a very interesting form: $a_{ij} = n+1-min(i,j)$, e.g. $n=5$: enter image description here

  • $\begingroup$ (+1) To be explicit, you might add that $s$ can be shown to be asymptotically Normal in many ways, of which perhaps the most appealing is to note that it has the distribution of twice the mean of $n-1$ iid squared Normal variates--so just apply the CLT to them. $\endgroup$
    – whuber
    Commented Jun 8, 2020 at 21:26

Small pedantic note: Below I changed the coefficient into $1/\sqrt{n-1}$ otherwise the limiting distribution will be a degenerate distribution (zero variance). In that case one would also need to subtract the mean of the $z_i^2$. That means, only a scaled and shifted sum like $\sum_{i=1}^{n-1} \frac{(x_{i+1}-x_i)^2-1}{\sqrt{n-1}}$ will approach a normal distribution.

Distribution of $s(\mathbf{z})$ as a linear sum of chi-squared variables

The sum $s(\mathbf{z}) = \frac{1}{\sqrt{n-1}}\sum_{i=1}^{n-1}(z_i)^2$ is similarly distributed as the sum $s(\mathbf{y}) = \frac{1}{\sqrt{n-1}}\sum_{i=1}^{n-1}(y_i)^2$ where the $y_i$ are $n-1$ independent normally distributed variables with variance $\lambda_i = 2 + 2 \cos(\frac{i}{n}\pi)$ $$s(\mathbf{z}) \sim s(\mathbf{y}) \quad \text{where} \quad y_i \sim N\left(0,\lambda_i \right)$$


  • The variance of $s(\mathbf{z})^2$ is equal to $1/\sqrt{n-1}$ times the sum of the variances of the individual terms $y_i^2$ (which relate to scaled $\chi_{(1)}^2$ distributions or to gamma distributions).

    For the individual terms which are the squares of normal distributed variables we have $$\begin{array}{rcl} \text{var}(y_i^2) &=& 2 \text{var}(y_i)^2 \\ &=&2\left( 2 + 2 \cos\left(\frac{i}{n}\pi\right)\right)^2 \\ &=& 12 + 4\cos\left(2\frac{i}{n}\pi\right) + 16 \cos\left(\frac{i}{n}\pi\right) \end{array}$$ and for the sum $$\begin{array}{rcl} \text{var}[s(\mathbf{z})^2] &=& \frac{1}{n-1} \sum_{i=1}^{n-1}12 + { 4\cos\left(2\frac{i}{n}\pi\right)} \overbrace{ + 16 \cos\left(\frac{i}{n}\pi\right)}^{\substack{\text{these terms cancel}\\\text{ due to symmetry}}} \\ &=& \frac{1}{n-1} ( \sum_{i=1}^{n-1}12 + 4 \underbrace{\sum_{i=1}^{n-1}\cos\left(2\frac{i}{n}\pi\right))}_{=-1} \\& =& \frac{12(n-1) -4}{n-1} \\&\approx& 12 \end{array}$$

    where we used this to derive that the sum of cosines is equal to -1.

  • We can not express the probability density function in a closed-form but we can express the cumulants of the distribution $\kappa_k(s(\mathbf{z}))$ in terms of the cumulants of a single chi-squared variable $\kappa_k(\chi_{(1)}^2)$. For increasing $n$ the 1st order cumulant will go to infinity (so to make the limiting distribution a normal distribution you should not only change the factor $1/(n-1)$ but also subtract the mean), the 2nd cumulant will approach to $12$ and the other, higher order cumulants, will approach zero, which means that you approach a normal distribution.

    $$\kappa_k(s(\mathbf{z})) = \kappa_k(\chi_{(1)}^2)\frac{1}{\sqrt{n-1}^k} \sum_{i=1}^{n-1} \lambda_i^k \approx \kappa_k(\chi_{(1)}^2)\frac{n-1}{\sqrt{n-1}^k} \int_0^1 ( 2 + 2 \cos(x\pi))^k dx$$

    Maybe there is a more direct way to use some version of the CLT for a sum of independent variables that only differ by a scaling constant instead of manually computing the cumulants. But I couldn't find one.

Why is it equivalent to a sum of chi-squared variables?

(see also here Show that the distribution of $x'Ax$ is linear combination of chi-squared)

Example for $n=3$

In the case of $n=3$ then the $z_1$ and $z_2$ are distributed like a multivariate normal distribution with a negative correlation. Geometrically it looks like an elongated shape.


n <- 10^4
x <- matrix(rnorm(3*n),3)
z1 <- x[2,]-x[1,]
z2 <- x[3,]-x[2,]
plot(z1,z2, xlab = expression(z[1]), ylab = expression(z[2]))

We can express the square in terms of alternative variables $Y_1 = \sqrt{0.5}(Z_1-Z_2) \sim N(0,3)$ and $Y_2 = \sqrt{0.5}(Z_1+Z_2) \sim N(0,1)$

$$Z_1^2 + Z_2^2 = 0.5(Z_1-Z_2)^2 + 0.5(Z_1+Z_2)^2 = Y_1^2 + Y_2^2$$

Note that the $Y_i$ are independent. So the distribution is similar to the distribution of a sum of independent squared normal distributed variables, but with different variance.

Generalized for all $n$

More generally $z_i$ is a multivariate normal distribution (any linear combination of the $z_i$ is a linear combination of the $x_i$ which is a normal distributed variable).

The variance of each $z_i$, being the sum of two standard normal variables, is $2$. The covariance of two neighbouring variables is $-1$ (which you can find with covariance of sums). So the covariance matrix is like:

$$\Sigma = \begin{bmatrix} {2} & -1 & 0 & \dots & 0 &0\\ -1 & 2 & -1 & 0& \dots & 0 \\ 0 & -1 & {2}& \dots & 0& 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0& 0 & \dots & {2} & -1 \\ 0 & 0& 0 & \dots & -1 & {2} \end{bmatrix}$$

In this general case we can do the same as for $n=3$ and restate the dependent $z_i^2$ as a sum of independent squared normal variables $y_i^2$. We use the same geometrical interpretation and rotate the distribution (keeping the radial distance invariant) and the distribution of $z_i$ is equivalent to rotated $y_i$ which have a variance that relate to the eigenvalues of the covariance matrix $\Sigma$. These eigenvalues will be in between 0 and 4 (see more about that below).

These eigenvalues follow a cosine function

$$\lambda_i = 2 + 2 \cos(\frac{i}{n}\pi)$$

for $1\leq i\leq n+1$. Which can be derived from the general description of eigenvalues of triadagonal Toeplitz matrices (as mentioned by Aksakal in the comments, you can see previous edits of this post for an alternative derivation of that relation with cosines)

  • $\begingroup$ (1) @Aksakal's answer already addresses all the issues you raise at the end. (2) This quadratic form has a rich history: it is the Schläfli matrix $A_{n-1}.$ $\endgroup$
    – whuber
    Commented Jun 8, 2020 at 22:20
  • $\begingroup$ For Tridiagonal Toeplitz matrix like yours the eigenproblem solution is known to be cosine, but it seems a bit different than yours $\endgroup$
    – Aksakal
    Commented Jun 9, 2020 at 13:31
  • $\begingroup$ I think using a vector $v=(x_1,z_1,\dots,z_{n-1})$ should simplify the rotation (PCA) argument. When you apply PCA on $z_i$ only the eigen vectors are complicated somewhat. Converting between $v$ and $x$ vectors is trivial, and the variance is preserved. I'm experimenting with this approach $\endgroup$
    – Aksakal
    Commented Jun 9, 2020 at 18:00
  • $\begingroup$ @Aksakal, I believe that I am getting to the eigenvalues without much troubles now (I do not see what your change with the vector $v$ improves). In your previous comment you said that there was a difference but maybe it is because I used in a previous version of my post $(2 - 2\cos(x))^2= 16 \sin^4(x/2)$ thus not a cosine and a sine instead, but effectively the same. $\endgroup$ Commented Jun 9, 2020 at 18:11
  • $\begingroup$ Great answer."Maybe there is a more direct way to use some version of the CLT for a sum of independent variables that only differ by a scaling constant(...) But I couldn't find one." Here it is: en.wikipedia.org/wiki/Lindeberg%27s_condition ; and in this case the maximum of the variances of any of the independent random variable in the sum is negligible compared to the sum of the variances, so Lindenberg's condition holds if and only if the CLT holds. To show Linderberg's condition is satisfied here, it seems Cauchy-Schwarz inequality, a bound on the 4th moment of $\chi^2_1$ is enough. $\endgroup$
    – jlewk
    Commented Jun 15, 2020 at 12:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.