calculate variance of conditional sum of Normals

Question

Let $\quad W \sim N(\mu_0, \sigma_0^2), \quad X_1 \sim N(\mu_1, \sigma_1^2), \quad X_2 \sim N(\mu_2, \sigma_2^2)$

denote independent random variables, and let $c$ denote a constant.

I want to find the cdf and variance of $Z$, where:

$$ Z = \begin{cases}W + X_1 & \text{if } W \leq c \\ W + X_2 & \text{if } W > c \end{cases}$$

I know from a previous question (see statsSE) that the pdf of $Z$ is:

$\frac{\frac{e^{-\frac{\left(\mu _0+\mu _1-z\right){}^2}{2 \left(\sigma _0^2+\sigma _1^2\right)}} \left(\text{erf}\left(\frac{\sigma _1^2 \left(g-\mu _0\right)+\sigma _0^2 \left(g+\mu _1-z\right)}{\sqrt{2} \sigma _0 \sigma _1 \sqrt{\sigma _0^2+\sigma _1^2}}\right)+1\right)}{\sqrt{\sigma _0^2+\sigma _1^2}}-\frac{e^{-\frac{\left(\mu _0+\mu _2-z\right){}^2}{2 \left(\sigma _0^2+\sigma _2^2\right)}} \left(\text{erf}\left(\frac{\sigma _2^2 \left(g-\mu _0\right)+\sigma _0^2 \left(g+\mu _2-z\right)}{\sqrt{2} \sigma _0 \sigma _2 \sqrt{\sigma _0^2+\sigma _2^2}}\right)-1\right)}{\sqrt{\sigma _0^2+\sigma _2^2}}}{2 \sqrt{2 \pi }}$

My code is as follows:

hsol = (((1/
     Sqrt[Subscript[σ, 0]^2 + 
       Subscript[σ, 
         1]^2]) (Exp[(-(-z + Subscript[μ, 0] + 
            Subscript[μ, 1])^2)/(2*(Subscript[σ, 0]^2 +
            Subscript[σ, 1]^2))]*(1 + 
       Erf[((g - z + Subscript[μ, 1])*
            Subscript[σ, 0]^2 + (g - 
              Subscript[μ, 0]) Subscript[σ, 
              1]^2)/(Sqrt[2]*Subscript[σ, 0]*
           Subscript[σ, 1]*
           Sqrt[Subscript[σ, 0]^2 + 
             Subscript[σ, 1]^2])]))) - ((1/
     Sqrt[Subscript[σ, 0]^2 + 
       Subscript[σ, 
         2]^2]) (Exp[(-(-z + Subscript[μ, 0] + 
            Subscript[μ, 2])^2)/(2*(Subscript[σ, 0]^2 +
            Subscript[σ, 2]^2))]*(-1 + 
       Erf[((g - z + Subscript[μ, 2])*
            Subscript[σ, 0]^2 + (g - 
              Subscript[μ, 0]) Subscript[σ, 
              2]^2)/(Sqrt[2]*Subscript[σ, 0]*
           Subscript[σ, 2]*
           Sqrt[Subscript[σ, 0]^2 + 
             Subscript[σ, 2]^2])]))))/(2 Sqrt[2 Pi]);
domain[hsol] = {z, -Infinity, Infinity};
Prob[hsol]
Var[z, hsol]

But after some time mathematica (using the mathStatica package) returns an expression containing an integral or stops calculating. I assume it is too difficult to calculate these expressions. Is there anyway make to circumvent this behaviour?

Answer 1

The Variance Problem

Find $Var(Z)$, where:

$$Z = \begin{cases}W + X_1 & \text{if } W \leq c \\ W + X_2 & \text{if } W > c \end{cases}$$

where $\quad W \sim N(\mu_0, \sigma_0^2), \quad X_1 \sim N(\mu_1, \sigma_1^2), \quad X_2 \sim N(\mu_2, \sigma_2^2)$ are independent, and $c$ is a constant.

Solution

Let $f(w, x_1,x_2)$ denote the joint pdf of $(W, X_1, X_2)$ which, by independence, is simply the product of the 3 individual Normal pdf's:

Define:

Then $Var(Z)$ is:

where I am also using the Var function from the mathStatica package for Mathematica which the OP is using. All done.

Monte Carlo check

It is worthwhile doing a Monte Carlo check. To illustrate, let:

Wdata  = RandomReal[NormalDistribution[16, 6], 1000000];
X1data = RandomReal[NormalDistribution[3, .1], 1000000];
X2data = RandomReal[NormalDistribution[2, 2],  1000000];

and

Zdata = MapThread[ If[#1 < 12, #1 + #2, #1 + #3] &, {Wdata, X1data, X2data}]; 

The sample variance of the $Z$ data is:

Variance[Zdata]

35.3591

The exact symbolic $Var(Z)$, given the same parameter values, is:

which 'matches' neatly. Looks fine :)