# Distribution of sum of two independent normals conditional on one of them

Assume $X$ and $Y$ are iid $N(0,1)$. I am looking for a "neat" expression for $$P\left(\frac{X+Y}{\sqrt{2}}>c\,\Biggl|\,X<c\right).$$ Related question seem to be discussed here or here, but if they already give the answer to my question, I do not see it.

Simulation suggests it is around 3% for $c$ the 95% normal quantile:

c <- qnorm(0.95)
cprob.num <- rep(NA,50000)

for (i in 1:reps){
X <- rnorm(1)
Y <- rnorm(1)
cprob.num[i] <- (X+Y)/sqrt(2) > c & X<c
}

mean(cprob.num)/0.95 # 0.03117895

• With a simple change of variables, you can reduce this to a comparable question for a bivariate Normal distribution $(U,V)$ where $U=X-c$, $V=(X+Y)/\sqrt{2}-c$, $c=0$, and $(U,V)$ are correlated. Then by following the analysis at Using the methods described at stats.stackexchange.com/a/71303/919 or using standard formulas for bivariate normal distributions, the question is reduced to measuring an angle when $c=0$ and otherwise requires numerical integration. That might explain why you have been unable to find a formula.
– whuber
Commented Jun 23, 2017 at 13:43
• Thank you. It looks like it will take me a while to digest the linked answer. Commented Jun 23, 2017 at 14:00
• X and Y have bivariate normal distribution: $X,Y\sim (\begin{bmatrix}0\\0\end{bmatrix},\begin{bmatrix}1&0\\0&1\end{bmatrix})$ However due to constraint $X<c$ we are dealing with truncated multivariate normal distribution. We need to impose linear transformation $A=[\sqrt{2}, \sqrt{2}]$ on it in order to find the result (using R tmvnorm package). Howevere I did not find even characteristic function formula for truncated normal in order to determine linear transformation consequences so it should be not straightforward. However its application could solve the problem. Commented Jun 23, 2017 at 19:52
• @whuber, I have so far been unsuccessful to turn your help into a solution. In that vein, does Mathematica - see wolfies answer below - overlook some way to produce a closed-form result? Commented Jun 24, 2017 at 7:57

Given: $X$ and $Y$ are independent standard Normals with pdf's $\phi(.)$ and cdf's $\Phi(.)$.

Since $X$ and $Y$ are independent, the joint pdf of $\big((X \; \big|\;X<c), \; Y\big)$ is $f(x,y) = {\large\frac{\phi(x)}{\Phi(c)}} \phi(y)$:

where Erf[.] denotes the error function.

Part 1: The pdf of $Z = X+Y \; | \; X<c$

Given $f(x,y)$, consider the transformation $(Z = X+Y, V=Y)$.

If $X <c$ and $Z = X+Y$, then $Z < c + Y$. That is, $Z < c + V$. This dependency is invoked in the following line using the Boole statement. Then the joint pdf of $(Z,V)$, say $g(z,v)$ can be obtained with:

... where I am using the Transform function from mathStatica/Mathematica to automate the nitty-gritties using the Method of Transformations (Jacobian etc).

The pdf of $Z$ that we seek is simply the marginal pdf of $Z$:

... which is our desired closed form solution.

The following diagram plots the pdf of $Z$ (i.e. the sum of 2 independent Normals, conditional on one of them) for six different vales of parameter $c$:

Part 2: Find $P\left(\frac{X+Y}{\sqrt{2}}>c\,\Biggl|\,X<c\right)$

To find $P\left(\frac{X+Y}{\sqrt{2}}>c\,\Biggl|\,X<c\right)$, integrate the above pdfZ over $(\sqrt2 c, \infty)$ wrt $z$.

Alternatively, $P\left(\frac{X+Y}{\sqrt{2}}>c\,\Biggl|\,X<c\right)$ can be obtained directly from the first step by :

... where I am using the Prob function from mathStatica/Mathematica to automate the nitty-gritties. This solution can be written in conventional notation as:

$$\frac{1}{\Phi(c)} \quad \int_{-\infty}^c \phi(x) \; \Phi \left(x-\sqrt{2} c\right) \, dx$$

While the probability does not appear to have a convenient closed-form, it is nevertheless a useful and practical result that is reduced to integrating a single variable. In particular:

a) when $c = 0$, the solution simplifies to $\frac14$

b) for other $c$ values, replace Integrate with NIntegrate for a solution via numerical integration in a single variable, which works very nicely. For instance, here is a plot of the desired probability, as a function of the truncation point $c$:

• Thanks a lot! May I suggest that the notation $(X|X<c, Y)$ is potentially confusing, as it might suggest that $Y$ is part of what is conditioned on? Commented Jun 24, 2017 at 7:55
• I don't think the notation "joint pdf of $( X|X<c, Y)$" is ambiguous. (a) it refers to the joint pdf; If the | nested both $X<c$ and $Y$, it would be univariate, not joint. (b) if the alternative was intended, then it would be written: $(X | \{X<c, Y=y\} )$. Either way, if you are suggesting $( (X|X<c), Y)$ is clearer, can certainly do that. Commented Jun 24, 2017 at 14:03
• You are right that the careful reader has nothing to worry about in your notation. That said, it is maybe just me, but I stumbled over this upon first reading your +1 answer. Your last suggestion avoids this, at of course the cost of more cumbersome notation. Probably it is fine as it is together with this discussion. Commented Jun 24, 2017 at 14:08

Sorry for not delivering the details, but $$\int_{-\infty}^c \phi(x) \; \Phi(x-\sqrt{2} c) \, dx = 2T(c, \sqrt{2}-1)$$ where $T$ is the Owen $T$-function.

This function is available in Mathematica/Wolfram and in the R package OwenQ.

library(OwenQ)
pr <- function(c){
2*OwenT(c, sqrt(2)-1) / pnorm(c)
}
curve(Vectorize(pr)(x), from=-6, to=6)


Alternatively you can get the Owen $T$-function with the help of the cdf of the noncentral Student distribution:

owenT <- function(h, a) 1/2*(pt(a, 1, h*sqrt(1+a^2)) - pnorm(-h))


But this implementation is not reliable for large values of the noncentrality parameter h*sqrt(1+a^2).

• That is very cool. The only other time I have come across the OwenT function was as the CDF of a skewNormal ... which makes me think that the pdf obtained here can likely be characterised as a skew-Normal distribution, with some tweaking. Commented Jun 25, 2017 at 20:10
• Also, the Mma implementation of OwenT appears to have some numerical instability problems - may require sending in a bug report Commented Jun 25, 2017 at 20:15
• @wolfies Where did you observe this numerical instability ? Commented Jun 25, 2017 at 20:31
• MyFunc[b_] := With[{c = b}, NIntegrate[(1 + Erf[c - z/2])/ (E^(z^2/4)*(2*Sqrt[Pi]*(1 + Erf[c/Sqrt[2]]))), {z, Sqrt[2]*c, Infinity}]]; Plot[{OwenT[c, 1 - Sqrt[2]], (-2^(-1))* CDF[NormalDistribution[0, 1], c]*MyFunc[c]}, {c, -3.8, -3}, PlotRange->All] Commented Jun 30, 2017 at 18:11
• @wolfies Indeed. The numerical value of OwenT[3.3, Sqrt[2]-1] returned by Mathematica 11.1 is not correct. But the one returned by Wolfram online is the right one (0.0002073998). Commented Jul 8, 2017 at 18:47