Mean of a function of $K$ Gaussian random variables

Question

Consider $K$ Gaussian random variables $Z_1, Z_2,\ldots,Z_K$ that are mutually independent with $Z_1 \sim \mathcal{N}(\mu,\sigma^2)$ and the remaining random variables $\sim \mathcal{N}(0,\sigma^2)$. Let $m_1({\bf Z})$ and $m_2({\bf Z})$ be functions of ${\bf Z} = (Z_1,\ldots,Z_K)$ such that $m_1({\bf Z})$ is the random variable with the largest magnitude and $m_2({\bf Z})$ is the random variable with the second largest magnitude among the $K$ random variables. I am trying to find the mean of the ratio $R = \frac{|m_2({\bf Z})|}{|m_1({\bf Z})|}$.

Q1: The event $m_1({\bf Z})=0$ is equivalent to the vector ${\bf Z}= {\bf 0}$ and hence it occurs with zero probability. So, I have assumed that $R$ is well-defined. The range of $R$ is $(0,1]$. So, the mean of $R$ must exist. Am I correct in assuming that the mean exists?

Q2: If the mean exists, the expression for the mean is given by $$ E[R]= \int\limits_{-\infty}^{\infty}\ldots\int\limits_{-\infty}^{\infty} \frac{|m_2({\bf z})|}{|m_1({\bf z})|} \frac{1}{(2\pi \sigma^2)^{K/2}} \exp \{-\frac{(z_{1}-\mu)^2+\sum_{i=2}^{K}z_{i}^2}{2\sigma^2}\} dz_{1}\ldots dz_{K}. $$ In general, for $ K>2$ this integral cannot be evaluated in closed form and numerical integration is required. Is numerical integration even feasible for values of $K >2 $?

Q3: In my attempt to answer Q2, I expressed the $K$-dimensional integral in terms of repeated one-dimensional integrals. My approach was to consider the following three cases.

1. $Z_1$ has the largest magnitude i.e., $m_1({\bf Z}) = Z_1$. Let $P_1$ be the conditional probability of $m_1({\bf Z}) = Z_1$ and $E_1$ be the conditional mean of symbol ratio when $m_1({\bf Z}) = Z_1$.

2. $Z_1$ has the second largest magnitude i.e., $m_2({\bf Z}) = Z_1$. Let $P_2$ be the conditional probability of $m_2({\bf Z}) = Z_1$ and $E_1$ be the conditional mean of symbol ratio when $m_2({\bf Z}) = Z_1$.

3. $Z_1$ has neither the largest nor the second largest magnitude. Let $P_3$ be the conditional probability that this is the case and $E_3$ be the conditional mean when this occurs.

Then, $E[R] = E_1P_1+E_2P_2+E_3P_3$. I was able to express each of these terms as repeated one-dimensional integrals with appropriate limits. Is this the right approach? The expressions for each of the terms are given below.

$E_1 = \int\limits_{-\infty}^{\infty} \frac{1}{|z_1|} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{(z_1-\mu)^2}{2 \sigma^2}} (K-1) \int\limits_{-|z_1|}^{|z_1|} \frac{|z_2|}{\sqrt{2\pi \sigma}} e^{-\frac{z_2^2}{2 \sigma^2}} \left[ \int\limits_{-|z_2|}^{|z_2|} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{z^2}{2\sigma^2}} dz\right]^{K-2} dz_{2} dz_{1}$, where $\Phi(.)$ is the standard Gaussian distribution function.

Using change of variables $x=\frac{z_1-\mu}{\sigma}$ and $y=\frac{z_2}{\sigma}$, we have $E_1 =\int\limits_{-\infty}^{\infty} \frac{1}{|x+\frac{\mu}{\sigma}|} \phi(x) (K-1) \int\limits_{-|x+\frac{\mu}{\sigma}|}^{|x+\frac{\mu}{\sigma}|} |y| \phi(y) \left[ 2\Phi( |y|) -1 \right]^{K-2} dy dx$

Similarly,

$P_1 = (K-1) \int\limits_{-\infty}^{\infty} \phi(x) \int\limits_{-|x+\frac{\mu}{\sigma}|}^{|x+\frac{\mu}{\sigma}|} \phi(y) \left[ 2\Phi( |y|) -1 \right]^{K-2} dy dx $,

$E_2=\int\limits_{-\infty}^{\infty} \frac{1}{|x|} (K-1) \phi(x) \int\limits_{-|x|-\frac{\mu}{\sigma}}^{|x|-\frac{\mu}{\sigma}} |y+\frac{\mu}{\sigma}| \phi(y) \left[ 2\Phi\left( \Bigg|y+\frac{\mu}{\sigma} \Bigg|\right) -1\right]^{K-2} dy dx$,

$P_2 = (K-1)\int\limits_{-\infty}^{\infty} \phi(x) \int\limits_{-|x|-\frac{\mu}{\sigma}}^{|x|-\frac{\mu}{\sigma}} \phi(y) \left[ 2\Phi\left( \Bigg|y+\frac{\mu}{\sigma} \Bigg|\right) -1 \right]^{K-2} dy dx$,

$E_3=(K-1)(K-2)\int\limits_{-\infty}^{\infty} \frac{1}{|x|} \phi(x) \int\limits_{-|x|}^{|x|} |y| \phi(y) \left[ \Phi\left( |y|-\frac{\mu}{\sigma} \right) -\Phi \left( -|y|-\frac{\mu}{\sigma} \right)\right] \left[ 2\Phi\left( |y|\right) -1\right]^{K-3}dydx$,

and

$P_3=(K-1)(K-2)\int\limits_{-\infty}^{\infty} \phi(x) \int\limits_{-|x|-\frac{\mu}{\sigma}}^{|x|-\frac{\mu}{\sigma}} \phi(y) \left[ \Phi\left( |y|-\frac{\mu}{\sigma} \right) -\Phi \left( -|y|-\frac{\mu}{\sigma} \right)\right] \left[ 2\Phi\left( |y|\right) -1\right]^{K-3} dy dx$.

I was able to find the sample average by generating large number of random samples of $\bf Z$ according to the distribution. I am trying to figure out if numerical integration is feasible at all. On using matlab to evaluate the expressions for $E_1$ through $E_3$ and $P_1$ through $P_3$, the answers do not match with the sample average.

This is the last puzzle in my research and your inputs will be highly appreciated.

Answer 1

The approach shown in the question basically looks right, but I think the notation might be getting in the way.

It is tempting to simplify things by comparing the largest of the $K-1$ standard Normal variates to $Z_1$. In order to know which is the largest and which is the second largest, we need to know what the two largest of the $K-1$ standard Normals are, because in some cases they will be the two largest. This reduces the question to finding the joint distribution of those two largest, because $Z_1$ is independent of them.

In general, let $F$ be any continuous distribution function and $f$ its derivative (the PDF). The PDF of the joint distribution of the two largest of $N$ iid draws from this distribution can be written

$$\binom{N}{N-2, 1, 1} F(y)^{N-2} f(y) f(x) dxdy ,\quad y \le x.$$

(This can rigorously be derived, but intuitively we need to consider the chance that one of the $N$ values lies in the interval $[x, x+dx)$, another lies in $[y, y+dy)$, and all the rest are less than or equal to $y$. The multinomial coefficient $\binom{N}{N-2, 1, 1} = N(N-1)$ counts the distinct ways this can occur while $F(y)^{N-2}$ is the chance that $N-2$ independent values lie in the interval $(-\infty, y)$, $f(y)dy$ is the chance that one value (independent of the first $N-2$ values) lies in $[y, y+dy)$, and $f(x)dx$ similarly is the chance that the last value lies in $[x, x+dx)$.)

Let, then, $F$ be the distribution function for the magnitude of the standard Normal (otherwise known as a Half-Normal distribution). From the symmetry around $0$ it is evident that $F(x) = 2\Phi(x) - 1$ when $x \ge 0$ (and otherwise is zero), where $\Phi$ is the CDF of the Standard Normal distribution.

Note that the ratio is invariant under a change of scale in all variables simultaneously. This implies that its distribution depends only on $\mu/\sigma$ rather than $\mu$ and $\sigma$ separately. For convenience, we may choose $\sigma=1$ throughout; the results will depend only on $\mu$. So, let $G_\mu$ be the distribution of the magnitude of a Normal$(\mu, 1)$ distribution. Writing $g_\mu$ for its derivative (the pdf), it is just as evident that $g_\mu(x) = \phi(x-\mu) + \phi(x+\mu)$ ($x\ge 0$) where $\phi$ is the standard Normal pdf. From the independence of $Z_1$ with respect to the other variates, we thereby obtain the joint (trivariate) distribution governing $(X,Y,Z)$ where $X$ is the largest of $N=K-1$ standard Normals, $Y$ is the second largest, and $Z=Z_1$:

$$\text{PDF}(x,y,z) = \binom{K-1}{K-3, 1, 1} F(y)^{K-3} f(y) f(x) g_\mu(z) dx dy dz,\quad y \le x.$$

To find the expectation of the ratio, we look at three cases (rather than $3!=6$, because it is guaranteed that $Y \le X$):

1. $Z \gt X$. The second largest is $X$ and the largest is $Z$, whence the ratio is $X/Z$. The integration is limited to $0\le y\le x \lt z \lt \infty$.

2. $Y \lt Z \le X$. Now the ratio is $Z/X$ and the integration is limited to $0\le y\lt z \le x\lt \infty$.

3. $Z \le Y$. The ratio is $Y/X$ and the integration is limited to $0\le z \le y \le x\lt \infty$.

Here is Mathematica code to compute these integrals numerically. First, the PDFs and CDF needed for the Half-Normals:

f[x_] := 2 PDF[NormalDistribution[0, 1]][x];
ff[x_] := 2 CDF[NormalDistribution[0, 1], x] - 1;
g[x_, m_] := Evaluate[Sum[PDF[NormalDistribution[u, 1]][x], {u, {-m, m}}]];

Now the integration, expressed in terms of $K$ and the parameter $\mu$:

ratio[k_Integer, m_] := Multinomial[k-3, 1, 1]  Plus @@ {
  NIntegrate[ff[y]^(k-3) f[x] f[y] g[z, m] x/z, {x,0,Infinity}, {y,0,x}, {z,x,Infinity}],
  NIntegrate[ff[y]^(k-3) f[x] f[y] g[z, m] z/x, {x,0,Infinity}, {y,0,x}, {z,y,x}],
  NIntegrate[ff[y]^(k-3) f[x] f[y] g[z, m] y/x, {x,0,Infinity}, {y,0,x}, {z,0,y}]};

(By comparing these expressions to the foregoing formulas and to the calculations in the question, it should be evident that they are all expressing the same thing.)

By evaluating these results at a sequence of values of $\mu$ we can interpolate the rest. Here are values for $K=3$ (red) and $K=20$ (blue), with the numerically integrated values shown as points:

With about the same amount of computation (around one minute), we can simulate the same results in R. Here is the simulation function:

sim <- function(n.iter, k, mu, sigma=1) {
  x <- matrix(abs(rnorm(n.iter * k, mean=c(rep(0, k-1), mu), sd=sigma)), nrow=k);
  z <- apply(x, 2, function(y) y[order(y, decreasing=TRUE)[1:2]]);
  z[2,] / z[1,]
}

With it we may estimate the values to approximately three decimal figures of precision using n.iter=10^5 and plot the results.

set.seed(17)
mean3 = sapply(0:30, function(m) mean(sim(10^5, k=3, mu=m)))
mean20 = sapply(0:30, function(m) mean(sim(10^5, k=20, mu=m)))

blue <- rgb(.25, .3, .5)
red <- rgb(.6, .2, .25)
par(mfrow=c(1,1))
plot(c(0,30), c(0,1), type="n", xlab="Mu",ylab="Mean ratio")
points(0:30, mean20, pch=19, col=blue)
lines(0:30, mean20, lwd=2, col=blue)
points(0:30, mean3, pch=19, col=red)
lines(0:30, mean3, lwd=2, col=red)

It agrees with the numerical integrals, of course. With the simulation we may recover the full distribution of ratios for any $K$ and $\mu$. For high accuracy, let's use more iterations. This series varies $\mu$ from $0$ to $30$ for $K=3$, using one million iterations:

par(mfrow=c(2,2))
tmp <- sapply(c(0,3,10,30), function(m) hist(sim(10^6, k=3, mu=m), main=m, xlab="Ratio"))

Answer 2

First of all, note that $R$ is not necessarily in $(0,1)$. For instance suppose that $k=3$ and a sorted realisation of $(Z_1,Z_2,Z_3)$ is $(-3,-2,1)$, then $R=2/1=2$. Did you mean the sorted absolute values?

A numerical simulation suggests that the mean does not exist. See

R.mean <- function(ns,k,mu,sigma){
R <- vector()
for(i in 1:ns){
sort.samp <- sort(rnorm(k,mu,sigma))
R[i] <- abs(sort.samp[k-1])/abs(sort.samp[k])
}
return(mean(R))
}

# Run it several times
R.mean(10000,3,0,1)

If you meant the sorted absolute values, then the expectation of $R$ does exist since the variable is bounded. This can be calculated in R as follows.

R.mean <- function(ns,k,mu,sigma){
R <- vector()
for(i in 1:ns){
sort.samp <- sort(abs(rnorm(k,mu,sigma)))
R[i] <- sort.samp[k-1]/sort.samp[k]
}
return(mean(R))
}

R.mean(10000,3,0,1)