6

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
Christoph Hanck
  • 25,948
  • 3
  • 57
  • 106
  • 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 https://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 – 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 07:57

2 Answers2

4

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

enter image description here

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:

enter image description here

... 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$:

enter image description here

... 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$:

enter image description here


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 :

enter image description here

... 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$:

enter image description here

wolfies
  • 6,963
  • 1
  • 22
  • 27
  • 1
    Thanks a lot! May I suggest that the notation $(X|X – Christoph Hanck Jun 24 '17 at 07:55
  • 1
    I don't think the notation "joint pdf of $( X|X – wolfies Jun 24 '17 at 14:03
  • 1
    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
4

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)

enter image description here

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

Stéphane Laurent
  • 17,425
  • 5
  • 59
  • 101
  • 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 08 '17 at 18:47