Expectation function vs Integral function

Question

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:

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:

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

The second way is based on Expectation function:

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:

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

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 *)