1
$\begingroup$

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}];
$\endgroup$
  • $\begingroup$ Just curious, did you try solving them while just writing out the integrals, rather than using the statistics tools in MMA? $\endgroup$ – MikeY Jan 23 at 13:39
  • $\begingroup$ Hi @MikeY, not sure I understand the question. Basically I did not know about the Expectation function. I started with Method 3 but using Integrate rather than NIntegrate. Integrate took a long time and sometimes did not find a solution at all. So started to use NIntegrate, even though it also has issues. I posted a separate question on why Integrate does not work and davidgstork suggested the Expectation function, which I tried. But then came across the problems described in this question and davidgstork suggested to post a separate question. $\endgroup$ – mariodrumblue Jan 24 at 20:05
  • $\begingroup$ Here is the link to my first post link $\endgroup$ – mariodrumblue Jan 24 at 20:06
2
$\begingroup$

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 *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.