Defining symbolic expectation function

Question

I want to define the standard expectation $E$ operator in Mathematica. In particular, I want it to satisfy,

$E[ c + a \cdot X_1^{n_1} X_2^{n_2} X_3^{n_3} X_4^{n_4} + b \cdot Y_1^{m_1} Y_2^{m_2} Y_3^{m_3} Y_4^{m_4} ] = c + a \cdot E[ X_1^{n_1} X_2^{n_2} X_3^{n_3} X_4^{n_4}] + b \cdot E[ Y_1^{m_1} Y_2^{m_2} Y_3^{m_3} Y_4^{m_4} ]$

$c, a, b$ are deterministic constants, and where $n_i$ and $m_j$ are positive integers and, in particular, can take on the value $0$. I do not want to enforce a particular distribution function on these random variables. (It's even better if the code can incorporate arbitrary number of products and powers of the random variables).

In Mathematica, I will associate the random variables $X_i, Y_j$ as functions, say, of the form randX[i], randY[j], so any value that does not match the randX[i] and randY[j] form will be regarded as deterministic constants. And the resulting moments should look like, in Mathematica,

expect[c + a * randX[1]^n1 * randX[2]^n2 * randX[3]^n3 * randX[4]^n4 + 
  b * randY[1]^m1 * randY[2]^m2 * randY[3]^m3 * randY[4]^m4]

c + a * expect[randX[1]^n1 * randX[2]^n2 * randX[3]^n3 * randX[4]^n4] + 
  b * expect[randY[1]^m1 * randY[2]^m2 * randY[3]^m3 * randY[4]^m4] 

The difficulty I'm having is that I find it hard to write patterns and rules that take on the zero-value powers of those random variables.

Answer 1

Here is one way of defining the linearity of your expectation value.

Updated in response to comments

To allow for powers as they arise in the definition of the variance, I just added rules that expand such powers when possible (it is only possible when there is a sum under the power, that's why I restrict the pattern expect[Power[expr_Plus, n_]] to expressions expr with head Plus), and recognize any complete expectation value expect[...] as a constant that can be pulled out from under a surrounding expect. In the treatment of Power, the case n == 1 is also matched (because Power has a definition for a Default second argument set equal to 1).

Clear[expect]

expect[expr_Plus] := Map[expect, expr]

expect[Times[x_, y__]] /; (FreeQ[x, randX[_] | randY[_]]) := 
 x expect[Times[y]]

expect[Times[x_expect, y__]] := x expect[Times[y]]

expect[expr_?(FreeQ[#, randX[_] | randY[_]] &)] := expr

expect[Power[expr_Plus, n_]] := expect[Expand[Power[expr, n]]]

expect[Power[x_expect, n_]] := x^n

expect[
 c + a*randX[1]^n1*randX[2]^n2*randX[3]^n3*randX[4]^n4 + 
  b*randY[1]^m1*randY[2]^m2*randY[3]^m3*randY[4]^m4]

There is nothing special that needs to be done for powers of 0. Powers of 0 make the random variable disappear, so that case is covered. In deciding whether to pull out constants, I use a PatternTest of the form expr_?(...) check that the random variables don't appear in expr (unless wrapped by expect, which when raised to any power counts as a constant, too - that's the last definition).

With the rules for Power and Times involving expect itself, we also get the desired simplification of the variance:

expect[(randX[1] - expect[randX[1]])^2]

-expect[randX[1]]^2 + expect[randX[1]^2]

Answer 2

I recently ran into this same question and came up with a slightly more convenient solution based on Jens answer that I think is different enough to be worth sharing. First I define two functions, RandomVariableQ and MakeRandomVariable

(*atoms are not RVs unless explicitly marked*) 
RandomVariableQ[expr_?AtomQ] := False

(*recurse down expression tree until we reach an atom or expectation value*)
RandomVariableQ[expr_ /; (!AtomQ[expr]) && (Head[expr] =!= AngleBracket)] := Or @@ Map[RandomVariableQ,List @@ expr]

(*expectation values are not RVs*)
RandomVariableQ[expr__AngleBracket] := False

(*make a symbol into a RV*)
MakeRandomVariable[x__Symbol] := (# /: RandomVariableQ[#] := True;) & /@ List[x]

MakeRandomVariable can be used to mark a particular symbol as being a random variable that cannot be factored out of an expectation. RandomVariableQ returns true if an expression contains a random variable that is not wrapped in an expectation value.

I then define AngleBracket as the expectation value

(*Linear*)
AngleBracket[expr_Plus] := AngleBracket /@ expr
AngleBracket[Times[x_, y__]] /; !RandomVariableQ[x] := x AngleBracket[Times[y]]

(*Expectation values only operate on RVs*)
AngleBracket[expr_ /; !RandomVariableQ[expr]] := expr

(*Expand before taking expectation values*)
AngleBracket[Power[expr_Plus, n_]] := AngleBracket[Expand[Power[expr, n]]]
AngleBracket[Times[expr_Plus, x_?RandomVariableQ]] := AngleBracket[Expand[expr x]]

Note here I am using RandomVariableQ instead of FreeQ. This allows the code to correctly recognize that arbitrary functions of expectation values can be factored out of an expectation value, i.e $\left\langle f(\left\langle x \right\rangle) \right \rangle = f(\left\langle x \right\rangle)$

We can test the code by computing the first 5 central moments of a random variable

MakeRandomVariable[x];
Table[AngleBracket[(x - AngleBracket[x])^n], {n, 2, 5}] // TableForm

which produces

$$ \begin{array}{c} \left\langle x^2\right\rangle -\langle x\rangle ^2 \\ \left\langle x^3\right\rangle -3 \left\langle x^2\right\rangle \langle x\rangle +2 \langle x\rangle ^3 \\ \left\langle x^4\right\rangle -4 \left\langle x^3\right\rangle \langle x\rangle +6 \left\langle x^2\right\rangle \langle x\rangle ^2-3 \langle x\rangle ^4 \\ \left\langle x^5\right\rangle -5 \left\langle x^4\right\rangle \langle x\rangle +10 \left\langle x^3\right\rangle \langle x\rangle ^2-10 \left\langle x^2\right\rangle \langle x\rangle ^3+4 \langle x\rangle ^5 \\ \end{array} $$

as desired. This is an improvement on the current solution which cannot simplify expressions such as expect[expect[randX[1]]^2 randX[1]] because it does not recognize that expect[randX[1]]^2 can be pulled out of the expectation value. This also works as expected with multiple random variables.