# Computing the Expectation

The following solution to this problem was found by @Wolfies:

(E^(-((-2 z + μ1 + μ2)^2/(2 (σ1^2 + σ2^2)))) Sqrt[2/π] (-1 + Erf[((z - μ2) σ1^2 + (-z + μ1) σ2^2)/
(Sqrt[2] σ1 σ2 Sqrt[σ1^2 + σ2^2])]))/((-1 + Erf[(μ1 - μ2)/(Sqrt[2] Sqrt[σ1^2 + σ2^2])]) Sqrt[σ1^2 + σ2^2])


Note that $$Z = \frac{X + Y}{2}$$ where $$X$$ and $$Y$$ are independent random variables that both follow the Normal distribution, with possibly different means and variances (denoted by subscripts $$1$$ for $$X$$ and $$2$$ for $$Y$$).

I am now trying to compute the mean of the above, after first creating a probability distribution as follows:

\[ScriptCapitalD] = ProbabilityDistribution[((E^(-((-2 z + μ1 + μ2)^2/(2 (σ1^2 + σ2^2)))) Sqrt[2/π] (-1 + Erf[((z - μ2) σ1^2 + (-z + μ1) σ2^2)/(Sqrt[2] σ1 σ2 Sqrt[σ1^2 + σ2^2])]))/((-1 + Erf[(μ1 - μ2)/(Sqrt[2] Sqrt[σ1^2 + σ2^2])]) Sqrt[σ1^2 + σ2^2]),{z, -∞, ∞}]


and then,

Expectation[a, a \[Distributed] \[ScriptCapitalD]]


but Mathematica won't give me an answer. Is this beyond Mathematica? (I wonder if the mathStatica add-on will help)?

• Sometimes (and certainly no always) it helps to look at special cases to see if some general form becomes apparent. If $\sigma_1=\sigma_2=\sigma$, then the expectation is $\frac{\sigma (\text{$\mu $1}+\text{$\mu $2}) \text{erfc}\left(\frac{\text{$\mu $1}-\text{$\mu $2}}{2 \sigma }\right)}{\sqrt{2}-2 \sigma \text{erf}\left(\frac{\text{$\mu $1}-\text{$\mu $2}}{2 \sigma }\right)}$. I'm a big fan of mathStatica (bought a copy when it first came out maybe 18 years ago) and it might help. I was wrong about your previous question having an essentially closed-form solution but maybe... – JimB Apr 22 at 19:26
• @JimB careful, the given probability distribution is not normalized. When $\sigma_1=\sigma_2$ the expectation value is $(\mu_1+\mu_2)/2$. – Roman Apr 22 at 19:32
• @Roman. Uh, oh. I thought the previous question did give a normalized density. Thanks. I'll check on that. – JimB Apr 22 at 19:35
• @Roman There's a missing set of parentheses in the OP's text. After adding those in the pdf integrates to 1 and one gets the expected value that you gave. – JimB Apr 22 at 20:37

Evaluating the Integral for the conditional distribution is not required to find the conditional mean, as you can exploit the fact that you are using a linear transformation of a bivariate distribution, and the fact that the mean of a truncated normal has a nice closed form.

The idea is to decompose your transformation into two partial transformations, which, when applied in sequence, are equivalent to your full transformation.

We'll chose the first partial transformation so that the first marginal will be our difference, but the other marginal will be uncorrelated, and by joint normality, independent. This will allow us to truncate the first marginal and not effect the other.

We can then compute the mean vector after our first partial transformation and marginal truncation and then apply the second partial linear transformation to obtain the mean vector of the conditional distribution.

This technique actually allows you to find the mean and variance of the marginal distribution of a bivariate normal conditional on some truncation of the other marginal for any bivariate normal distribution. For bivariate normal distributions with small covariance, the conditional distribution is approximately normal, and you can use this technique to approximate the joint CDF for the bivariate normal, which has no nice closed form. For more information, see https://apps.dtic.mil/dtic/tr/fulltext/u2/a125033.pdf

In[1]:= means = Array[μ, 2]

Out[1]= {μ[1], μ[2]}

In[2]:= covar = DiagonalMatrix[Array[σ, 2]]

Out[2]= {{σ[1], 0}, {0, σ[2]}}

In[3]:= trans = {{1, -1}, {1/2, 1/2}}

Out[3]= {{1, -1}, {1/2, 1/2}}

In[4]:= ptransa = {{1, -1}, {1, a}}

Out[4]= {{1, -1}, {1, a}}

In[5]:= acovar = (ptransa.covar.Transpose[ptransa])[[1, 2]]

Out[5]= σ[1] - a σ[2]

In[6]:= ptrans = ptransa /. Solve[acovar == 0, a][[1]]

Out[6]= {{1, -1}, {1, σ[1]/σ[2]}}

In[7]:= ftrans = trans.Inverse[ptrans] // FullSimplify

Out[7]= {{1,
0}, {-(1/2) + σ[1]/(σ[1] + σ[2]), σ[
2]/(σ[1] + σ[2])}}

In[8]:= pmeans = ptrans.means

Out[8]= {μ[1] - μ[2], μ[
1] + (μ[2] σ[1])/σ[2]}

In[9]:= pcovar = ptrans.covar.Transpose[ptrans] // FullSimplify

Out[9]= {{σ[1] + σ[2],
0}, {0, (σ[1] (σ[1] + σ[2]))/σ[2]}}

In[10]:= dist =
TruncatedDistribution[{-Infinity, 0},
NormalDistribution[pmeans[[1]], Sqrt[pcovar[[1, 1]]]]]

Out[10]= TruncatedDistribution[{-∞, 0},
NormalDistribution[μ[1] - μ[2],
Sqrt[σ[1] + σ[2]]]]

In[11]:= tmeans = {Mean[dist], pmeans[[2]]} // FullSimplify

Out[11]= {μ[1] - μ[2] - (
E^(-((μ[1] - μ[2])^2/(2 (σ[1] + σ[2])))) Sqrt[
2/π] Sqrt[σ[1] + σ[2]])/
Erfc[(μ[1] - μ[2])/(
Sqrt[2] Sqrt[σ[1] + σ[2]])], μ[
1] + (μ[2] σ[1])/σ[2]}

In[12]:= cmean = ftrans.tmeans // FullSimplify

Out[12]= {μ[1] - μ[2] - (
E^(-((μ[1] - μ[2])^2/(2 (σ[1] + σ[2])))) Sqrt[
2/π] Sqrt[σ[1] + σ[2]])/
Erfc[(μ[1] - μ[2])/(
Sqrt[2] Sqrt[σ[1] + σ[2]])],
1/2 (μ[1] + μ[2]) + (
E^(-((μ[1] - μ[2])^2/(
2 (σ[1] + σ[2])))) (-σ[1] + σ[2]))/(
Sqrt[2 π]
Erfc[(μ[1] - μ[2])/(
Sqrt[2] Sqrt[σ[1] + σ[2]])] Sqrt[σ[
1] + σ[2]])}

In[13]:= cmean[[2]]

Out[13]= 1/2 (μ[1] + μ[2]) + (
E^(-((μ[1] - μ[2])^2/(
2 (σ[1] + σ[2])))) (-σ[1] + σ[2]))/(
Sqrt[2 π]
Erfc[(μ[1] - μ[2])/(
Sqrt[2] Sqrt[σ[1] + σ[2]])] Sqrt[σ[
1] + σ[2]])


As @Roman pointed out the expectation is $$(\mu_1+\mu_2)/2$$ when $$\sigma_1=\sigma_2$$:

\[ScriptCapitalD] = ProbabilityDistribution[Exp[-(-2 z + μ1 + μ2)^2/(2 (σ1^2 + σ2^2))] *
Sqrt[2/π] (-1 + Erf[((z - μ2) σ1^2 + (-z + μ1) σ2^2)/(Sqrt[2] σ1 σ2 Sqrt[σ1^2 + σ2^2])])/
(Sqrt[σ1^2 + σ2^2] (-1 + Erf[(μ1 - μ2)/(Sqrt[2] Sqrt[σ1^2 + σ2^2])])),
{z, -∞, ∞}, Assumptions -> {σ1 > 0, σ2 > 0}];

Expectation[z, z \[Distributed] \[ScriptCapitalD] /. {σ1 -> σ, σ2 -> σ}]
(* (μ1 + μ2)/2 *)


But I think that when $$\sigma_1 \neq \sigma_2$$, then you'll need to use numerical methods for the expectation. (And hopefully I'm wrong about that.)

NExpectation[z, z \[Distributed] \[ScriptCapitalD] /. {μ1 -> 1, μ2 -> 3, σ1 -> 1, σ2 -> 2}]
(* 2.22026 *)