# 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 Jun 23 '17 at 13:43
• Thank you. It looks like it will take me a while to digest the linked answer. – Christoph Hanck Jun 23 '17 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. – Bogdan Jun 23 '17 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? – Christoph Hanck Jun 24 '17 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? – Christoph Hanck Jun 24 '17 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. – wolfies Jun 24 '17 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. – Christoph Hanck Jun 24 '17 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. – wolfies Jun 25 '17 at 20:10
• Also, the Mma implementation of OwenT appears to have some numerical instability problems - may require sending in a bug report – wolfies Jun 25 '17 at 20:15
• @wolfies Where did you observe this numerical instability ? – Stéphane Laurent Jun 25 '17 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] – wolfies Jun 30 '17 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). – Stéphane Laurent Jul 8 '17 at 18:47