Mean of a function of $K$ Gaussian random variables

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.

• Q3 would be more readable were you to adopt a more congenial and abstract notation. There's no need to keep reproducing the details of the normal pdf; simple symbols like $f$ and $F$ for the pdf and cdf will do just fine and, by reducing much of the extraneous clutter, help reveal what is really going on. Note, too, that it is clear a prior that the ratio depends only on $\mu/\sigma$, so you don't need both parameters. You might also find it a helpful guide first to obtain the means in the limiting cases $\mu/\sigma=0$ and $\mu/\sigma\to\infty$.
– whuber
Mar 13, 2013 at 19:00
• As you sample, you can tell whether you are in case 1, 2, or 3. So, you should be able to check which of your estimates for $E_i$ and $P_i$ match the simulation. In case you have only made a sign error, this should narrow your search considerably. Mar 13, 2013 at 21:02
• I have a solution involving nine numerical integrations but haven't time to write it up. Using the theory of order statistics it's straightforward to compute the bivariate distribution of the largest two of $K-1$ iid Half-Gaussians. Multiplying that by the pdf of the absolute value of a Gaussian of mean $\mu$ gives the trivariate pdf. The integral can be broken into just three regions of integration in each of which the ratio has a simple expression, thus reducing the problem to three $3$d integrations. Calculations take about $2$ seconds in Mathematica 8 and agree with simulation.
– whuber
Mar 14, 2013 at 0:01
• @whuber I don't understand why you would need the bivariate distribution of the largest two of $K-1$ iid Half-Gaussians. Mar 14, 2013 at 1:41
• @Sneha Because it can occur that two of the Half-Gaussians are the two largest values.
– whuber
Mar 14, 2013 at 3:22

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"))


• Thanks so much for your comprehensive solution! I arrived at a similar-looking result using a rigorous proof. I shall post my solution when I find the time tomorrow. Mar 14, 2013 at 5:21
• @Sneha I hope you're not suggesting my approach is anything less than rigorous: please do not mistake economy of notation for lack of rigor. (If you're referring to my parenthetical remark about "this can be rigorously derived," that's because the formulae for joint distributions of order statistics are so routine and well-known that I felt there was no need to go fully into the details here, nor is there any reason for you to do so either.) I do wonder, though, what might constitute a "rigorous proof" of a numerical integration: after all, it has no closed form.
– whuber
Mar 14, 2013 at 5:35
• Oh no, I did not mean it that way. I meant a different way of looking at the same thing. I apologize for using the wrong word. Your help is much appreciated! Mar 14, 2013 at 9:13

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)

• Ratio R is the ratio of second largest to largest magnitude. So, for your example it would be 2/3. I was able to obtain the sample average as shown in your second pseudo-code sample. I am trying to figure out if numerical integration is feasible is or not. Thanks! Mar 13, 2013 at 17:14
• Armin, the question is clear that absolute values ("magnitudes") are meant. Incidentally, simulations like those shown here don't help much in concluding whether means exist or not. They can suggest a mean does not exist when the variation in sample mean from one simulation to another is huge, or when that variation does not get lower as a function of larger sample size: but you haven't conducted either of those studies.
– whuber
Mar 13, 2013 at 18:53
• @whuber Thanks for re-clarifying what Sneha already did about the "magnitudes". Indeed, suggesting is not a proof, as you have noticed. Sneha: The second code (not pseudocode) does a numerical integration via Monte Carlo integration for any $k$, you just have to tune an appropriate number of simulations ns.
– user21965
Mar 13, 2013 at 20:52
• A monte-carlo simulation is not a bad idea (because it can also obtain an estimate of the entire distribution), but to estimate the mean ratio to better than three significant figures requires over $10^6$ iterations and unfortunately your code won't scale well at those sizes. Consider rewriting it to take advantage of R's built-in vectorization of calculations: that will speed it up by an order of magnitude at ns=10^6 and progressively more for larger numbers of iterations.
– whuber
Mar 13, 2013 at 21:47