# Mean and Variance of a function of a multivariate normal

I have a multivariate normally-distributed random variable:

$X \sim \mathcal{N}(\mu, \Sigma)$

And another RV which is a deterministic elementwise function of this variable (producing another random vector of the same dimensionality):

$Y=f(X)$ .

In my case $f(X) = \max(0, X)$, which truncates negative parts of the distribution for each element of $X$ to $X=0$, so that $p(Y_i=0)=p(X_i<0)$. But the ideal solution would work for any function.

I'd like to find $\text{E}(Y)$ and $\text{Var}(Y)$. Is there an analytical way to do this?

Suppose that $$Z \sim \mathcal N(0,1)$$ and $$W = \max(0,Z)$$. Note that $$P(W = 0) = \frac{1}{2}$$ so $$W$$ is like a continuous random variable plus a point mass at 0. This means that the CDF $$F_W$$ of $$W$$ is $$F_W(w) = \begin{cases} 0 \hspace{15mm} w < 0 \\ \frac{1}{2} \hspace{14mm} w = 0 \\ \Phi(w) \hspace{8mm} w > 0 \\ \end{cases}$$ with a discontinuity at 0.

We want the pdf of this. Let $$\delta_0$$ be the probability measure defined by $$\delta_0(A) = I(0 \in A)$$ where $$I$$ is the indicator function. Let the Lebesgue measure be $$m$$.

Consider some measurable set $$A \subset \mathbb R$$. Note that $$(m + \delta_0)(A) = 0 \implies P(W \in A) = 0$$ therefore $$P_W$$ is absolutely continuous with respect to $$m + \delta_0$$. This means we can obtain the Radon-Nikodym pdf $$f_W$$ of $$W$$, denoted by $$f_W = \frac{\text d P_W}{\text d(m + \delta_0)}.$$

Once we have this we can compute the expectation as follows: $$E(W) = \int_\mathbb R w \,\text dP_W = \int_\mathbb R w f_W(w) \,\text d(m + \delta_0)$$

$$=\int_{\{0\}} w f_W(w) \,\text d\delta_0 + \int_{(0, \infty)} w f_W(w) \,\text dm$$

$$= 0 + \int_{(0, \infty)} w f_W(w) \,\text dm$$

Because $$F_W(w) = \Phi(w)$$ on $$(0, \infty)$$ we know that $$f_W(w) = \phi(w)$$ on this interval therefore $$E(W) = \int \limits_0^\infty w \phi(w) \,\text dw = \frac{1}{\sqrt{2\pi}}.$$

We can confirm this via simulation:

z <- rnorm(1e5)
w <- sapply(z, function(v) max(c(0,v)))
mean(w)
1 / sqrt(2 * pi)


The second moment $$E(X^2)$$ can be found in a similar fashion. I was not rigorous with much of this, and your particular case is a little more complicated, but I think this shows the overall way to get an analytic solution out of this.

Update: $$X \sim \mathcal N(\mu, \sigma^2)$$

Let $$X \sim \mathcal N(\mu, \sigma^2)$$. Then

$$F_W(w) = \begin{cases} 0 \hspace{15mm} w < 0 \\ F_X(0) \hspace{6mm} w = 0 \\ F_X(w) \hspace{6mm} w > 0 \\ \end{cases}$$

so we find that $$E(W) = \frac{1}{\sqrt{2\pi\sigma^2}}\int \limits_0^\infty w \exp \left( -\frac{(w - \mu)^2}{2 \sigma^2}\right) \,\text dw.$$

We can write this in terms of the standard normal CDF $$\Phi$$: letting $$u = \frac{w-\mu}{\sigma}$$, we obtain $$E(W) = \frac{1}{\sqrt{2\pi\sigma^2}}\int \limits_0^\infty w \exp \left( -\frac{(w - \mu)^2}{2 \sigma^2}\right) \,\text dw = \frac{1}{\sqrt{2\pi}}\int \limits_{-\mu/\sigma}^\infty (\sigma u + \mu) \exp \left( -u^2/2\right) \,\text du$$

$$= \frac1{\sqrt{2\pi}}\left[ \sigma \int \limits_{-\mu/\sigma}^\infty u e^{-u^2/2} \,\text du + \mu \int \limits_{-\mu/\sigma}^\infty e^{-u^2/2} \,\text du \right].$$ Substituting $$z = u^2/2$$ in the left integral we find $$E(W) = \frac\sigma{\sqrt{2\pi}} e^{-\frac{\mu^2}{2\sigma^2}} + \frac\mu{\sqrt{2\pi}} \left[ \int \limits_{-\infty}^\infty e^{-u^2/2} \,\text du - \int \limits_{-\infty}^{-\mu/\sigma} e^{-u^2/2} \,\text du\right]$$

$$= \frac\sigma{\sqrt{2\pi}} e^{-\frac{\mu^2}{2\sigma^2}} + \mu \left( 1- \Phi(-\mu/\sigma) \right).$$

This expression still has an integral in it but at least it is a standard integral that we are comfortable with.

In R:

x <- rnorm(1e6, mean=2, sd=1.5)
w <- sapply(x, function(v) max(c(0,v)))

## checking CDF
mean(w <= 1)
pnorm(1, 2, 1.5)

## checking expectation
mean(w)

expectation <- function(y, mu, sigma) {
y * exp(-.5 * (y - mu)^2 / sigma^2) / sqrt(2 * pi) / sigma
}
integrate(expectation, 0, Inf, mu=2, sigma=1.5)

## using simpler expression
E_W <- function(mu, sigma)
sigma / sqrt(2 * pi) * exp(-.5 * (mu/sigma)^2) + mu * (1 - pnorm(-mu/sigma))
E_W(mu=2, sigma=1.5)


Update 2

Here's the complete multivariate case. Let $$X \sim \mathcal N (\mu, \Sigma)$$ and $$Y = (Y_1, \dots, Y_n)$$ for $$Y_i = \max(0, X_i)$$. What I have done above is sufficient for $$E(Y) = (E(Y_1), \dots, E(Y_n))$$ so what remains is the variance $$\Sigma_Y := Var(Y)$$. The diagonal of $$\Sigma_Y$$ we can do already: it is just $$E(Y_i^2) - E(Y_i)^2$$ and we have (albeit ugly) expressions for those.

So in order to fill in the rest of $$\Sigma_Y$$ we need $$Cov(Y_i, Y_j) = E(Y_iY_j) - E(Y_i)E(Y_j)$$, so the only missing part at this point is $$E(Y_iY_j)$$.

The CDF of $$(Y_i, Y_j)$$ has a form analogous to the univariate CDFs so it turns out that $$E(Y_i Y_j) = \int \limits_0^\infty \int \limits_0^\infty y_1 y_2 f_{X_i, X_j}(y_1, y_2)\,\text dy_1 \,\text dy_2.$$

Confirming this:

library(cubature)
library(mvtnorm)

## generating data
sx = 1; sy = 2; rho = .5
Sigma <- matrix(c(sx^2, sx * sy * rho, sx * sy * rho, sy^2), 2)
mu = c(-1, 2)
x <- MASS::mvrnorm(5e5, mu, Sigma)
y <- t(apply(x, 1, function(v) sapply(v, function(v2) max(c(0,v2)))))

## checking CDF
mean(y[,1] <= .5 & y[,2] <= 1.25)
mean(x[,1] <= .5 & x[,2] <= 1.25)
mvtnorm::pmvnorm(upper=c(.5,1.25), mean=mu, sigma=Sigma)

## checking E(Y_i Y_j)
mean(y[,1] * y[,2])

expectation <- function(vec, mu, sigma){
prod(vec) * mvtnorm::dmvnorm(vec, mu, sigma)
}
adaptIntegrate(expectation, lowerLimit = c(0,0), upperLimit = c(100, 100), mu=mu, sigma=Sigma)$integral  • Unfortunately it won't be nearly as nice when$\mu \neq 0$. You lose symmetry which makes it quite a bit sloppier. Oct 13, 2016 at 15:10 • True, but at least this is still a principled approach. – jld Oct 13, 2016 at 15:13 • Thanks for your response. However, I have trouble seeing how this could be extended to answer the original problem. The main difficulty would be in solving the integral when (a)$\mu$is not zero, and (b)$\Sigma$is not diagonal. Oct 13, 2016 at 15:37 • @Peter I've added a bit on$X \sim \mathcal N(\mu, \sigma)$. Also note that this does address your general case for$E( \bf Y)$since the expectation of a random vector is the vector of expectations of each element, and we've got$E(Y_i)$for each$i\$.
– jld
Oct 13, 2016 at 16:15
• @Peter I believe that I've completed extended this to the multivariate case at this point.
– jld
Oct 13, 2016 at 17:07