Sum of Independent Half-Normal Distributions with unequal variance

Question

$\DeclareMathOperator{\erf}{erf}$

I want to find the cumulative distribution function (CDF) of the random variable $Z = |X| + |Y|$ where $X$ and $Y$ are two independent random variables which are normally distributed with mean 0 but different variance $\sigma_X^2$ and $\sigma^2_{Y}$.

$|X|$ and $|Y|$ follow the so-called Half-Normal distribution, which is a special case of the Folded-Normal distribution, when $\mu = 0$. The probability density functions (PDF) of $|X|$ and $|Y|$ for $x > 0$ are defined as follows (wikipedia):

\begin{align} \label{pdf} f_{|X|}(x, \sigma_{X}) = \frac{\sqrt{2}}{\sigma_{X} \sqrt{\pi}} \exp \left(- \frac{x^2}{2\sigma_{X}^2}\right) && f_{|Y|}(y, \sigma_{Y}) = \frac{\sqrt{2}}{\sigma_{Y} \sqrt{\pi}} \exp\left(- \frac{y^2}{2\sigma_{Y}^2}\right) \end{align}

And their CDF:

\begin{align} \label{cdf} F_{|X|}(x, \sigma_{X}) = \erf\left(\frac{x}{\sigma_{X} \sqrt{2}}\right) && F_{|Y|}(y, \sigma_{Y}) = \erf\left(\frac{y}{\sigma_{Y} \sqrt{2}}\right) \end{align}

@DilipSarwate gave a very elegant solution for the case $\sigma_{Y} = \sigma_{X}$ here, but I am interested in the general case when this is not true, i.e. one of the variable has a higher weight in the final value (bigger scale).

So far, I have determined the PDF of $Z$ by doing the convolution of $f_{|X|}$ and $f_{|Y|}$. Since $f_{|X|}$ and $f_{|Y|}$ are defined only for $x > 0$, we can truncate the convolution to the interval $[0,z]$:

\begin{align} f_Z(z) &= \int_0^z f_{|Y|}(z-x) f_{|X|}(x) \, dx \\ f_Z(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \int_0^z \exp\left(\frac{-(z-x)^2}{2\sigma_Y^2}\right) \exp\left(\frac{-x^2}{2\sigma_X^2}\right) \, dx \\ f_Z(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \int_0^z \exp\left(\frac{-(x^2 - 2xz + z^2)}{2\sigma_Y^2}\right) \exp\left(\frac{-x^2}{2\sigma_X^2}\right) \, dx \\ f_Z(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2\sigma_Y^2}\right) \int_0^z \exp\left(\frac{-x^2+2xz}{2\sigma_Y^2}\right) \exp\left(\frac{-x^2}{2\sigma_X^2}\right) \, dx \\ f_Z(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2\sigma_Y^2}\right) \int_0^z \exp\left(\frac{-(\sigma_X^2 + \sigma_Y^2)x^2 + 2z\sigma_X^2 x}{2 \sigma_X^2 \sigma_Y^2}\right) \, dx \end{align}

We define $\sigma_Z = \sqrt{\left(\sigma_X^2 + \sigma_Y^2\right)}$:

\begin{align} f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2\sigma_Y^2}\right) \int_{0}^{z} \exp\left(\frac{-\sigma_Z^2 x^2 + 2z\sigma_X^2 x}{2 \sigma_X^2 \sigma_Y^2}\right) dx \end{align}

We complete the square $ax^2 + bx$ to the form $a(x-h)^2 + k$ with $h = -\frac{b}{2a}$ and $k = - \frac{b^2}{4a}$, i.e., $h = \frac{z\sigma_X^2}{\sigma_Z^2}$ and $k = \frac{(z \sigma_X^2)^2}{\sigma_Z^2}$. This results in:

\begin{align} f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2\sigma_Y^2}\right) \int_{0}^{z} \exp\left(\frac{-\sigma_Z^2\left(x - \frac{z \sigma_X^2}{\sigma_Z^2}\right)^2 + \frac{(z \sigma_X^2)^2}{\sigma_Z^2}}{2\sigma_X^2 \sigma_Y^2}\right) dx \\ f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2\sigma_Y^2}\right) \exp\left(\frac{(z\sigma_X^2)^2}{2\sigma_X^2 \sigma_Y^2 \sigma_Z^2}\right) \int_{0}^{z} \exp\left(\frac{-\sigma_Z^2 \left(x - \frac{z \sigma_X^2}{\sigma_Z^2}\right)^2}{2\sigma_X^2 \sigma_Y^2}\right) dx \\ f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2(\sigma_Z^2 - \sigma_X^2)}{2\sigma_Y^2 \sigma_Z^2}\right) \int_{0}^{z} \exp\left(\frac{-\sigma_Z^2 \left(x - \frac{z \sigma_X^2}{\sigma_Z^2}\right)^2}{2\sigma_X^2 \sigma_Y^2}\right) dx \\ f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2 \sigma_Z^2}\right) \int_{0}^{z} \exp\left(\frac{-\sigma_Z^2 \left(x - \frac{z \sigma_X^2}{\sigma_Z^2} \right)^2}{2\sigma_X^2 \sigma_Y^2}\right) dx \end{align}

\begin{align} f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2 \sigma_Z^2}\right) \int_{0}^{z} \exp\left(\frac{-\left(x - \frac{z \sigma_X^2}{\sigma_Z^2}\right)^2}{2\frac{\sigma_X^2 \sigma_Y^2}{\sigma_Z^2}}\right) dx \\ \end{align}

We add the term $\frac{\sqrt{2 \pi (\frac{\sigma_X \sigma_Y}{\sigma_Z})^2}}{\sqrt{2 \pi (\frac{\sigma_X \sigma_Y}{\sigma_Z})^2}}$:

\begin{align} f_{Z}(z) &= \frac{2}{\sigma_X \sigma_Y \pi} \exp\left(\frac{-z^2}{2 \sigma_Z^2}\right) \sqrt{2 \pi \left(\frac{\sigma_X \sigma_Y}{\sigma_Z}\right)^2} \int_{0}^{z} \frac{1}{\sqrt{2 \pi \left(\frac{\sigma_X \sigma_Y}{\sigma_Z}\right)^2}} \exp\left(\frac{-\left(x - \frac{z \sigma_X^2}{\sigma_Z^2}\right)^2}{2\left(\frac{\sigma_X \sigma_Y}{\sigma_Z}\right)^2}\right) dx \\ f_{Z}(z) &= \frac{2 \sqrt{2}}{\sigma_Z \sqrt{\pi}} \exp\left(\frac{-z^2}{2 \sigma_Z^2}\right) \int_{0}^{z} \frac{1}{\sqrt{2 \pi \left(\frac{\sigma_X \sigma_Y}{\sigma_Z}\right)^2}} \exp\left(\frac{-\left(x - \frac{z \sigma_X^2}{\sigma_Z^2}\right)^2}{2\left(\frac{\sigma_X \sigma_Y}{\sigma_Z}\right)^2}\right) dx \end{align}

As we can see, the expression in the integral represents the density of a Gaussian distribution with mean $\frac{z \sigma_X^2}{\sigma_Z^2}$ and standard deviation $\frac{\sigma_X \sigma_Y}{\sigma_Z}$. As a results we can evaluate it as follows:

\begin{align} f_{Z}(z) &= \frac{2 \sqrt{2}}{\sigma_Z \sqrt{\pi}} \exp\left(\frac{-z^2}{2 \sigma_Z^2}\right) \left[\Phi\left(z, \frac{z \sigma_X^2}{\sigma_Z^2}, \frac{\sigma_X \sigma_Y}{\sigma_Z}\right) - \Phi\left(0, \frac{z \sigma_X^2}{\sigma_Z^2}, \frac{\sigma_X \sigma_Y}{\sigma_Z}\right)\right] \end{align}

where $\Phi(x,m,s)$ is the value of the cumulative distribution function at $x$ of a Gaussian distribution with mean $m$ and standard deviation $s$.

Now we can compute the cumulative distribution function $F_Z$:

\begin{align} F_{Z}(z) &= \int_{0}^{z} f_{Z}(x) dx \\ F_{Z}(z) &= \frac{2 \sqrt{2}}{\sigma_Z \sqrt{\pi}} \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \left[\Phi\left(x, \frac{x \sigma_X^2}{\sigma_Z^2}, \frac{\sigma_X \sigma_Y}{\sigma_Z}\right) - \Phi\left(0, \frac{x \sigma_X^2}{\sigma_Z^2}, \frac{\sigma_X \sigma_Y}{\sigma_Z}\right)\right] dx \\ F_{Z}(z) &= \frac{2 \sqrt{2}}{\sigma_Z \sqrt{\pi}} \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \left[ \frac{1}{2} \left( 1 + \erf \left(\frac{x-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) \right) - \frac{1}{2} \left( 1 + \erf \left( \frac{-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) \right) \right] dx \\ F_{Z}(z) &= \frac{\sqrt{2}}{\sigma_Z \sqrt{\pi}} \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \left[ \erf \left( \frac{x-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) - \erf \left( \frac{-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) \right] dx \end{align}

Since the integral of the sum is equal to the sum of the integrals:

\begin{align} F_{Z}(z) &= \frac{\sqrt{2}}{\sigma_Z \sqrt{\pi}} \left[ \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( \frac{x-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) dx - \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( \frac{-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) dx \right] \\ F_{Z}(z) &= \frac{\sqrt{2}}{\sigma_Z \sqrt{\pi}} \left[ \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( \frac{x \left(1-\frac{\sigma_X^2}{\sigma_Z^2} \right)}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) dx - \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( \frac{-\frac{x \sigma_X^2}{\sigma_Z^2}}{\frac{\sigma_X \sigma_Y}{\sigma_Z} \sqrt{2}} \right) dx \right] \\ F_{Z}(z) &= \frac{\sqrt{2}}{\sigma_Z \sqrt{\pi}} \left[ \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( x \frac{\sigma_Y}{\sqrt{2} \sigma_Z \sigma_X } \right) dx - \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left(- x \frac{\sigma_X}{\sqrt{2} \sigma_Z \sigma_Y } \right) dx \right] \end{align}

We know that $\erf(x)$ is an odd function and that $\exp(x)$ is an even function, as a result, their product is odd, and thus:

\begin{align} F_{Z}(z) &= \frac{\sqrt{2}}{\sigma_Z \sqrt{\pi}} \left[ \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( x \frac{\sigma_Y}{\sqrt{2} \sigma_Z \sigma_X } \right) dx + \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left(x \frac{\sigma_X}{\sqrt{2} \sigma_Z \sigma_Y } \right) dx \right] \end{align}

This is where I am stuck. It seems that the integrals are not elementary. What can I do?

Just to verify my results in the case of equal variance, I am going to assume that $\sigma_X = \sigma_Y$. From 3.3 (41) in this paper, We know that:

$$ \int \exp \left( -a^2 x^2 \right) \left[ \erf(ax)\right]^n \, dx = \frac{\sqrt{\pi}}{2a(n+1)} \left[\erf(ax)\right]^{n+1} $$

Thus:

\begin{align} F^{\sigma_X = \sigma_Y}_{Z}(z) &= \frac{2 \sqrt{2}}{\sigma_Z \sqrt{\pi}} \int_{0}^{z} \exp\left(\frac{-x^2}{2 \sigma_Z^2}\right) \erf \left( \frac{x}{\sqrt{2} \sigma_Z } \right) dx \\ F^{\sigma_X = \sigma_Y}_{Z}(z) &= \frac{2 \sqrt{2}}{\sigma_Z \sqrt{\pi}} \frac{\sqrt{\pi}}{4 \frac{1}{\sqrt{2} \sigma_Z}} \erf(\frac{z}{\sqrt{2} \sigma_Z})^2 \\ F^{\sigma_X = \sigma_Y}_{Z}(z) &= \erf(\frac{z}{\sqrt{2} \sigma_Z})^2 \end{align}

Which is consistent with the solution from here.

However, as I said, I don't know how to solve the problem when $\sigma_X \neq \sigma_Y$. Any suggestions? I think I have no choice but to use a numerical approximation method, such as Simpson's rule.

Answer 1

$\DeclareMathOperator{\erf}{erf}$In the last step (so the one where you got stuck) you found two integrals of the form $$\int_{0}^{z} e^{-a^{2} x^{2}} {\rm erf}(bx)\, dx \quad.$$ I asked about this integral on Mathoverflow via this question. It turns out it can be written by means of a special function called Owen's T: $$\int_{0}^{z} e^{-a^{2} x^{2}} {\rm erf}(bx)\, dx = \frac{\tan^{-1}(b/a)}{a \sqrt{\pi}} - 2 \frac{\sqrt{\pi}}{a} T(\sqrt{2}az,b/a) \quad.$$ When you plug in the values $a = \sqrt{\frac{1}{2\sigma_{Z}^{2}}} = \frac{1}{\sqrt{2}\sigma_{Z}}, b_{1} = \frac{\sigma_{Y}}{\sqrt{2}\sigma_{Z}\sigma_{X}}$ for the first integral, and $b_{2} = \frac{\sigma_{X}}{\sqrt{2}\sigma_{Z}\sigma_{Y}} $ for the second integral (with the same value for $a$ as in the first one), you will obtain the cumulative distribution function of the convolution of two independent half-normal distributions.

When $\sigma_{X} = \sigma_{Y} = \sigma$, you can use the property $T(h,1) = \frac{1}{2} \Phi(h) (1- \Phi(h)) $ to find $F_{Z} (z) = {\rm erf} \big{(} \frac{z}{\sqrt{2}\sigma} \big{)}^{2}$, which I think is similar but not quite the same as the solution you found already, because you have $\sigma_{Z}$ in stead of $\sigma$, and $\sigma_{Z} = \sqrt{2} \sigma$ when $\sigma_{X} = \sigma_{Y} = \sigma$. So I'm not sure what's going on over here. Perhaps you know?

---

Edit: I made a small mistake. The cumulative distribution function and the error function are related as follows: $\Phi(x) = 1/2 + 1/2 \erf(x/\sqrt{2})$. I did not take into account the $\frac{1}{\sqrt{2}}$ factor. So the expression you found in case $\sigma_{1} = \sigma_{2} = \sigma$ is correct.

Answer 2

If $X \sim N(0,\sigma^2)$ and $Y \sim N(0,s^2)$, then given independence, the pdf of $Z = |X| + |Y|$ is given by the current developmental build of mathStatica as say $h(z)$:

$$h(z) = \frac{1}{\sqrt{s^2+\sigma ^2}} \sqrt{\frac{2}{\pi }} \exp \left( -\frac{z^2}{2 \left(s^2+\sigma ^2\right)} \right) \left(\text{erf}\left(\frac{\sigma z}{\sqrt{2} s \sqrt{s^2+\sigma ^2}}\right)+\text{erf}\left(\frac{s z}{\sqrt{2} \sigma \sqrt{s^2+\sigma ^2}}\right)\right)$$

This is just a one-liner to evaluate, and may be helpful to verify your manually derived solution $f(z)$.

The following is a quick Monte Carlo check of the theoretical solution above when $\sigma = \frac12$ and $s = 4$.

Looks fine. That just leaves integrating the pdf :)