Polynomial expectations of generic distributions

Question

Summary

Polynomial expectations depend only moments and cross moments of a multivariate distribution. I would like to use Expectation to compute polynomial expectations for generic distribution of which only the necessary moments are specified.

Some details

I would like to compute polynomial exectations with the function Expectation for variables distributed according to a generic distribution. Here is a toy example:

In[139]:= ClearAll[myDist]

myDist /: Moment[myDist[mu_, var_], 1] := mu
myDist /: Moment[myDist[mu_, var_], 2] := Moment[myDist[mu, var], 1]^2 + var
myDist /: Moment[myDist[0, var_], 4] := 
 Cumulant[myDist[0, var], 4] + 2 Moment[myDist[0, var], 2]
myDist /: Moment[myDist[0, var_], 3] := 0

This few lines provide enough information about the distribution myDist so that Mathematica can compute

In[144]:= Expectation[2 x^2, x \[Distributed] myDist[0, s^2], Method -> "Moment"]
Expectation[2 x^4, x \[Distributed] myDist[0, s^2], Method -> "Moment"]

Out[144]= 2 s^2

Out[145]= 2 (2 s^2 + Cumulant[myDist[0, s^2], 4])

I would like to carry out similar calculations with several random variables. For example.

In[147]:= ClearAll[a, b, x, y]
Expectation[
 a x + b y, {x \[Distributed] myDist[0, s^2], 
  y \[Distributed] myDist[0, s^2]}, Method -> "Moment"]
Expectation[(a x + b y)^2, {x \[Distributed] myDist[0, s^2], 
  y \[Distributed] myDist[0, s^2]}, Method -> "Moment"]

Out[148]= Expectation[
 a x + b y, {x \[Distributed] myDist[0, s^2], 
  y \[Distributed] myDist[0, s^2]}, Method -> "Moment"]

Out[149]= Expectation[(a x + b y)^2, {x \[Distributed] myDist[0, s^2], 
  y \[Distributed] myDist[0, s^2]}, Method -> "Moment"]

In this case Mathematica doesn't carry out the computation, even though it has all the information to do so. Namely

1. x and y are independent (this is implied in the specification {x\[Distributed]myDist[0,s^2],y\[Distributed]myDist[0,s^2]}, according to the definition of Expectation);

2. the result is completely determined by the first and second moment of x and y (Provided).

One of the problems is that Expectation doesn't know (or it is pre-emped by some other rule or evaluation)

In[153]:= Expectation[
  a x + b y, {x \[Distributed] myDist[0, s^2], 
   y \[Distributed] myDist[0, s^2]}, Method -> "Moment"] ===

 a Expectation[x, x \[Distributed] myDist[0, s^2], Method -> "Moment"] + 
  b Expectation[y, y \[Distributed] myDist[0, s^2], Method -> "Moment"]

Out[153]= False

Notice that all-is-well when the distribution is one that Mathematica knows

In[166]:= Expectation[
 a x + b y, {x \[Distributed] NormalDistribution[0, s^2], 
  y \[Distributed] NormalDistribution[0, s^2]}, Method -> "Moment"]

Out[166]= 0

One possible solution is to define an operator such as myExpectation that knows/uses linearity properties of the mathematical expectation, but I would rather not reinvent the wheel and leverage the power of the built in symbol Expectation (for example specializing a result to one of the built in distributions) as well as the ease with it deals with higher momenta and cumulants. I tried (with no success)

1. using ProductDistribution[{myDist[], 2}] instead of {x\[Distributed]myDist[0,s^2],y\[Distributed]myDist[0,s^2]};

2. definiting a generic distribution as suggested (here) with the symbol ProbabilityDistribution. In this case the delayed UpValues for the moments cannot be set.

Indeed

 myDist2[mu_, var_] =
 ProbabilityDistribution[myDistPDF[x, mu, var], {x, -Infinity, Infinity}]
myDist2 /: Moment[myDist[mu_, var_], 1] := mu

Out[164]= ProbabilityDistribution[
 myDistPDF[\[FormalX], mu, var], {\[FormalX], -\[Infinity], \[Infinity]}]

During evaluation of In[164]:= TagSetDelayed::tagnf: Tag myDist2 not found in Moment[ProbabilityDistribution[myDistPDF[\[FormalX],mu_,var_],{\[FormalX],-\[Infinity],\[Infinity]}],1]. >>

Out[165]= $Failed

If a HoldPattern the left-hand-side of the moment delayed assignment, the command executes with no error, but it the information won't be used by Expectation. The reason for this is that Expectation first evaluates its arguments and myDist2 evaluated to

ProbabilityDistribution[myDistPDF[\[FormalX],mu,var],{\[FormalX],-\[Infinity],\[Infinity]}]

and the upvalues of myDist are from then on "invisible" to Expectation.

Some related question on how to define an arbitrary distributions are here and here.

Answer 1

Provided you add

myDist /: Moment[myDist[mu_, var_], 0] := 1

The following (which is risky as it messes up with built in definitions) would work

freeQ[a_, b_] := FreeQ[a, b] /; Length[b] <= 1
freeQ[a_, b_] := And @@ Map[FreeQ[a, #] &, b] /; Length[b] > 1

Unprotect[Expectation];
Expectation[a__, b__, c__] := 
       Map[Expectation[#, b, c] &, a // Expand] /; Head[Expand[a]] == Plus
Expectation[a_ z_, b_, c__] := 
       a Expectation[z, b, c] /; freeQ[a, b /. u_ \[Distributed] _ -> u]
Expectation[a_ , {b__, c__}, d__] := 
       Expectation[a, b, d] /; freeQ[a, c /. u_ \[Distributed] _ -> u]
Expectation[a1_ a2_ , {b1__, b2__}, c__] := 
Expectation[a1, b1, c] Expectation[a2, b2, c] /; 
      (freeQ[a2, b1 /. u_ \[Distributed] _ -> u] && 
       freeQ[a1, b2 /. u_ \[Distributed] _ -> u] )
Expectation[a_ , {b_, c_}, d__] := 
       Expectation[a, c, d] /; freeQ[a, b /. u_ \[Distributed] _ -> u]
Protect[Expectation];

e.g.

 Expectation[(a x + b y)^2, {x \[Distributed] myDist[0, s^2], 
 y \[Distributed] myDist[0, s^2]}, Method -> "Moment"]

(* a^2 s^2+b^2 s^2 *)

 Expectation[(a x + b y)^3, {x \[Distributed] myDist[0, s^2], 
 y \[Distributed] myDist[0, s^2]}, Method -> "Moment"]

(* 0 *)

Once again this is risky and I have not checked all possible edge effects ! :-)