I am having trouble understanding how Mathematica computes a conditional expectation.

Consider random variables $v\sim N(0,\tau_v^{-1})$, $s_0=v+\epsilon_0$, $s_1=v+\eta+\epsilon_1$, where $\epsilon_0\sim N(0,\tau_0^{-1})$, $\eta \sim N(0,\tau_{\eta}^{-1})$, $\epsilon_1\sim N(0,\tau_1^{-1})$. For simplicity, $(v,\epsilon_0,\eta,\epsilon_1)$ are mutually independent, so $(v,s_0,s_1)$ are jointly normal.

It is easy to show that

$E[s_1|s_0]=\frac{\tau_0}{\tau_v+\tau_0}s_0$, $Var[s_1|s_0]=\frac{1}{\tau_v+\tau_0}+\frac{1}{\tau_{\eta}}+\frac{1}{\tau_1}$, and $E[v|s_0,s_1]=\frac{\tau_0(\tau_{\eta}+\tau_1)s_0+\tau_1\tau_{\eta}s1}{\tau_1(\tau_0+\tau_v)+(\tau_0+\tau_1+\tau_v)\tau_{\eta}}$.

Also, $s_1|s_0\sim N\big(E[s_1|s_0],Var[s_1|s_0]\big)$.

My goal is to compute $E\Big[E[v|s_0,s_1]\times {\bf\large 1}\{E[v|s_0,s_1]>0\}|s_0\Big]$. Below is my code in Mathematica.

First, I enter the expressions above:

enter image description here

Es1[s0_] := τ0/(τv + τ0) s0; (* E[s1|s0] *)
σs1s0 := Sqrt[1/(τv + τ0) + 1/τη + 1/τ1];  (* σ[s1|s0] *)
fs1[s0_] := PDF[NormalDistribution[Es1[s0], σs1s0], s1];  (* pdf of s1|s0 *)
Ev[s0_, s1_] := (τ0 (τ1 + τη) s0 + τ1 τη s1)/(τ1 (τ0 + τv) + (τ0 + τ1 + τv) τη);   (* E[v|s0,s1] *)

Next, I attempt to compute the conditional expectation in two ways. The first way is based on Integrate function:

enter image description here

\*SubsuperscriptBox[\(∫\), \(\(-τ0\) \((
\*FractionBox[\(1\), \(τ1\)] + 
\*FractionBox[\(1\), \(τη\)])\) s0\), \(+∞\)]\(Ev[
   s0, s1] fs1[s0] \[DifferentialD]s1\)\)

The second way is based on Expectation function:

enter image description here

Expectation[Ev[s0, s1] Boole[Ev[s0, s1] > 0] \[Conditioned] s0, 
 s1 \[Distributed] NormalDistribution[Es1[s0], σs1s0]]

The results are somewhat surprising to me. The first way (after FullSimplify) returns a very long expression: enter image description here

whereas the second way returns 0. Could you help me understand what is going on here?


1 Answer 1


They work out the same for me. I'm using Version 13.1.0 for Mac. It takes about 7 seconds on my MacBook Pro M1 Max.

assumptions = \[Tau]0 > 0 && \[Tau]1 > 0 && \[Tau]v > 0 && \[Tau]\[Eta] > 0 && s0 > 0;
bound = s0 \[Tau]0 (-(1/\[Tau]1) - 1/\[Tau]\[Eta]);

one = Integrate[
        Ev[s0, s1] fs1[s0], 
        {s1, bound, \[Infinity]}, 
        GenerateConditions -> False,
        Assumptions -> assumptions

two = Expectation[
        Ev[s0, s1] Boole[s1 > bound], 
        s1 \[Distributed] NormalDistribution[Es1[s0], \[Sigma]s1s0], 
        GenerateConditions -> False,
        Assumptions -> assumptions

Simplify[one == two, assumptions]
(* True *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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