Issues computing expectation

Question

I am interested in computing a conditional expectation, but not sure what is the best way to do it. All the methods I have used seem to have issues. Let $Y=X_1+X_2$ where $X_1$ and $X_2$ are independent random variables, both normally distributed. I am interested in computing $E[X_1 \mid Y \geq c]$ where $c$ is an arbitrary constant. This is $$ E[X_1 \mid Y>c] = \frac{1}{1-\Phi_Y(c)} \int_c^\infty \int_{-\infty}^{\infty} x_1 \cdot f_{X_1}(x_1) \cdot f_{X_2}(y-x_1) \quad dx_1 ~dy $$ In Mathematica, with $c=0.5$, I am comparing three ways of computing this expectation:

Method 1 - Using the symbolic Expectation function

z[m1_, s1_, m2_, s2_, c_] := 
  Expectation[
   x1 \[Conditioned] 
    c < x1 + x2, {x1 \[Distributed] NormalDistribution[m1, s1], 
    x2 \[Distributed] NormalDistribution[m2, s2]}];

Method 2 - Using NExpectation function

zN[m1_, s1_, m2_, s2_, c_] := 
  NExpectation[
   x1 \[Conditioned] 
    c < x1 + x2, {x1 \[Distributed] NormalDistribution[m1, s1], 
    x2 \[Distributed] NormalDistribution[m2, s2]}];

Method 3 - Using the formula above with NIntegrate

gN[m1_, s1_, m2_, s2_, 
   c_] := (1/(1 - 
       CDF[NormalDistribution[m1 + m2, Sqrt[s1^2 + s2^2]], 
        c])) NIntegrate[
    x1 PDF[NormalDistribution[m1, s1], x1] PDF[
      NormalDistribution[m2, s2], y - x1], {x1, -Infinity, 
     Infinity}, {y, c, Infinity}];

Methods (2) and (3) work, and perform almost the same in terms of computation time. However, method (1) will only work under some given parametrisation but not others, and if it works would understandably take longer. For instance, with $$ \{\mu_{X_1},\sigma_{X_1},\mu_{X_2},\sigma_{X_2},c\}=\{0,1,0,1,0.5\} $$ all methods work, but with $$ \{\mu_{X_1},\sigma_{X_1},\mu_{X_2},\sigma_{X_2},c\}=\{\mathbf{1,2,1,2},0.5\} $$ method 3 no longer works. Q1: Is there a reason why the method (1) only works under certain parametrizations?

Also, if I change the problem such that I want to compute the expectation within a region defined over a finite interval as in $$ E[X_1 \mid c_L<Y<c_H] = \frac{1}{\Phi_Y(c_H)-\Phi_Y(c_L)} \int_{c_L}^{c_H} \int_{-\infty}^{\infty} x_1 \cdot f_{X_1}(x_1) \cdot f_{X_2}(y-x_1) \quad dx_1 ~dy $$ the parametrisation also matters when comparing method 2 and method 3. For instance, with $$ \{\mu_{X_1},\sigma_{X_1},\mu_{X_2},\sigma_{X_2},c_L,c_H\}=\{0,1,0,1,0,0.5\} $$ method (2) also fails to find a solution, while method (3) works. With $$ \{\mu_{X_1},\sigma_{X_1},\mu_{X_2},\sigma_{X_2},c_L,c_H\}=\{0,1,0,1,\mathbf{-0.5},0.5\} $$ method (3) produces NIntegrate warning messages. Methods (1) and (2) don't work. Q2: Is there a reason for this? Is method (3) the best way for me to compute the expectation even though NIntegrate will struggle sometimes and give warning messages? My code with the finite interval below:

zN[m1_, s1_, m2_, s2_, cL_, cH_] := 
  NExpectation[
   x1 \[Conditioned] 
    cL < x1 + x2 <= cH, {x1 \[Distributed] NormalDistribution[m1, s1],
     x2 \[Distributed] NormalDistribution[m2, s2]}];

gN[m1_, s1_, m2_, s2_, cL_, 
   cH_] := (1/(CDF[NormalDistribution[m1 + m2, Sqrt[s1^2 + s2^2]], 
        cH] - CDF[NormalDistribution[m1 + m2, Sqrt[s1^2 + s2^2]], 
        cL])) NIntegrate[
    x1 PDF[NormalDistribution[m1, s1], x1] PDF[
      NormalDistribution[m2, s2], y - x1], {x1, -Infinity, 
     Infinity}, {y, cL, cH}];

Answer 1

Your question (and now my question, too) is "Why doesn't Method 1 work all of the time?" I don't know. But I do have a workaround.

On the hunch/wild guess that having a non-zero mean has something to do with it (maybe having two non-zero means makes things too complicated?), I modify your function z to shift the normal distributions so that one of the distributions has a zero mean and also change c to the appropriate value.

z2[m1_, s1_, m2_, s2_, c_] := 
 m1 + Expectation[x1 \[Conditioned] (c - 2 m1) < x1 + x2, 
{x1 \[Distributed] NormalDistribution[0, s1], 
 x2 \[Distributed] NormalDistribution[m2 - m1, s2]}]

z2[1, 2, 1, 2, 1/2]

$$\frac{2}{e^{9/64} \sqrt{\pi } \left(\text{erf}\left(\frac{3}{8}\right)+1\right)}+1$$

z2[1, 2, 1, 2, 1/2] // N
(* 1.6982 *)

zN[1, 2, 1, 2, 1/2]
(* 1.69819 *)