Conditional expectation

Question

Suppose I have a bivariate normal random vector $(Y,Z)'$ where $Y\sim N(μ1,\sigma_1^2)$, $Z\sim N(μ2,\sigma_2^2)$, and $Cov[Y,Z]=\sigma_{12}$. If I want Mathematica to compute the conditional expectation of $Y|Z=z$, then the following command

Expectation[Y \[Conditioned] Z == z, {Y, Z} \[Distributed] MultinormalDistribution[{μ1, μ2}, {{Subscript[\[Sigma],1]^2, Subscript[\[Sigma], 12]}, {Subscript[\[Sigma], 12], Subscript[\[Sigma], 2]^2}}]]

returns the correct answer:

$$E[Y|Z] = \mu1+ \frac{ \sigma_{12} }{\sigma_2^2}(z- \mu2).$$

Is there a way to have Mathematica compute $E[Y| |Z|=z]$? I have tried

Expectation[Y \[Conditioned] Abs[Z] == z, {Y, Z} \[Distributed] MultinormalDistribution[{μ1, μ2}, {{Subscript[\[Sigma],1]^2, Subscript[\[Sigma], 12]}, {Subscript[\[Sigma], 12], Subscript[\[Sigma], 2]^2}}]]

But as output Mathematica returns what I have written as input. More generally, I would like to know whether it is possible to have Mathematica compute $E[Y| |Z|=z, X=x, W=w]$, with $(Y,Z,X,W)$ multivariate normally distributed random variables.

Answer 1

If (as stated by @wolfies in a comment) the joint distribution is a bivariate normal, then here's a hand-waving approach (which I think could be made more concrete).

Because you can get the conditional mean for $Z=z$, then the conditional mean for $|Z|=z$ should be a weighted average of the conditional means for $Z=z$ and $Z=-z$ with the weights being the marginal densities of $Z$ at $z$ and $-z$:

dist = MultinormalDistribution[{μ1, μ2}, {{Subscript[σ, 1]^2, Subscript[σ, 12]}, {Subscript[σ, 12], Subscript[σ, 2]^2}}];

μNegative = Expectation[Y \[Conditioned] Z == -z, {Y, Z} \[Distributed] dist] // FullSimplify;
μPositive = Expectation[Y \[Conditioned] Z == z, {Y, Z} \[Distributed] dist] //    FullSimplify;
pdfNegative = PDF[NormalDistribution[μ2, Subscript[σ, 2]], -z];
pdfPositive = PDF[NormalDistribution[μ2, Subscript[σ, 2]], z];
μ = (μNegative pdfNegative + μPositive pdfPositive)/(pdfNegative + pdfPositive) // FullSimplify

As a somewhat feeble check, consider taking a large random sample from the bivariate distribution and then keep only those values where the realized values of the absolute value of $Z$ are "close" to $z$:

parms = {Subscript[σ, 12] -> 2/3, μ1 -> 2, μ2 -> 2, 
   z -> 1, Subscript[σ, 1] -> 1, Subscript[σ, 2] -> 1};
SeedRandom[12345];
yz = RandomVariate[dist /. parms, 1000000];
yz = Select[yz, (Abs[Abs[#[[2]]] - Abs[z]] /. parms) < 0.01 &];
Length[yz]
(* 5013 *)
Mean[yz[[All, 1]]]
(* 1.31002 *)
μ /. parms // N 
(* 1.30935 *)

Here is a generalization to find

$$E(X_1 |\,\,\, |X_2|=x_2, X_3=x_3, X_4=x_4, \ldots, X_n=x_n)$$

conditionalMean[n_] := Module[{mean, pdf, cm},
  mean = Expectation[X[1] \[Conditioned] (Table[X[i] == x[i], {i, 2, n}] /. List -> And),
    Table[X[i], {i, n}] \[Distributed] MultinormalDistribution[Table[μ[i], {i, n}], 
      Table[σ[Min[i, j], Max[i, j]], {i, n}, {j, n}]]];
  pdf = PDF[MultinormalDistribution[Table[μ[i], {i, 2, n}], 
     Table[σ[Min[i, j], Max[i, j]], {i, 2, n}, {j, 2, n}]], Table[x[i], {i, 2, n}]];
  cm = (mean*pdf + (mean /. x[2] -> -x[2])*(pdf /. x[2] -> -x[2]))/(pdf + (pdf /. x[2] -> -x[2]));
  cm /. σ[i_, i_] -> Subscript[σ, i]^2 /. σ[i_, j_] -> Subscript[σ, i, j] /.
    μ[i_] -> Subscript[μ, i] /. x[i_] -> Subscript[x, i]]

conditionalMean[2] // FullSimplify

conditionalMean[3] // FullSimplify

So there looks to be a pattern here. conditionalMean[4] shows up quickly but applying FullSimplify afterwards might take forever.