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?


1 Answer 1


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

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:


## 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
  • $\begingroup$ Unfortunately it won't be nearly as nice when $\mu \neq 0$. You lose symmetry which makes it quite a bit sloppier. $\endgroup$ Oct 13, 2016 at 15:10
  • $\begingroup$ True, but at least this is still a principled approach. $\endgroup$
    – jld
    Oct 13, 2016 at 15:13
  • $\begingroup$ 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. $\endgroup$
    – Peter
    Oct 13, 2016 at 15:37
  • $\begingroup$ @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$. $\endgroup$
    – jld
    Oct 13, 2016 at 16:15
  • $\begingroup$ @Peter I believe that I've completed extended this to the multivariate case at this point. $\endgroup$
    – jld
    Oct 13, 2016 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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