# Defining symbolic expectation function

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.

• Does the issue of $n_i = 0$ or $m_i = 0$ actually result in something that you would need to handle differently? Do those instances get changed to 1 such that by the time that the expect function would operate, randX[1]^5 * randX[2]^0 * randX[3]^4 would be changed to randX[1]^5 * randX[3]^4 ? – JimB Aug 16 '15 at 0:20
• @JimBaldwin No, not in particular. – user32416 Aug 16 '15 at 0:45
• Not sure if your response is for my first question or the second question. Maybe if you showed what rules you've constructed so far (and even the ones that don't do what you want them to do), that would be helpful. – JimB Aug 16 '15 at 0:57
• @JimBaldwin Thanks for the suggestions! I'm currently now working off of solution of Jens (see below) at the moment. In particular, I'm following the usual rules of algebra that $X_1^5 X_2^0 X_3^4 = X_1^5 X_3^4$. – user32416 Aug 16 '15 at 1:02
• Good. At point you'll need to convert things like expect[randX[1]^2] to the associated moments of the random variables that you're considering. This example (assuming that the moments exist) would result in $\sigma^2+\mu^2$ (if $\sigma^2$ is the variance of randX[1] and $\mu$ is the mean of randX[1]. Do you have a convention for naming the moments of X[i] and Y[i] and the expectations of their products? – JimB Aug 16 '15 at 1:07

## 2 Answers

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]

• Thanks for the answer. Is it possible to expand your solution? In particular, there are often cases when one needs to consider, say $Z := c + a X_1^{n_1} X_2^{n_2} X_3^{n_3} X_4^{n_4}$ and consider the computation $E[ (Z - E[Z])^2 ]$ (i.e. the variance of $Z$); and likewise, consider the other central moments of $Z$. In this case, your current code would generate something like expect[ expect[randZ] ] and does not pull out the deterministic expect[randZ] inside. Thanks (in particular, since I'm a novice to Mathematica)! – user32416 Aug 16 '15 at 1:00
• What @Jens gave you will do the job if you use the Expand function. z = c + a*randX[1]^n1*randX[2]^n2*randX[3]^n3*randX[4]^n4; Simplify[expect[Expand[z^2]] - Expand[expect[z]^2]] results in a^2 (-expect[randX[1]^n1 randX[2]^n2 randX[3]^n3 randX[4]^n4]^2 + expect[randX[1]^(2 n1) randX[2]^(2 n2) randX[3]^(2 n3) randX[4]^( 2 n4)]). – JimB Aug 16 '15 at 1:24
• Thanks again! This does indeed work --- however, it seems to require the user to know the algebraic identity that $Var(Z) := E[ (Z - E[Z])^2 ] = E[Z^2] - ( E Z )^2$, which is fine for the variance. But I'm not sure how does one think about using Expand when one considers say, the 3rd centralized moment (for skewness) $E[ (Z - E[Z])^3 ]$, as expanding that would still involve a computation of a term like $E [ E[Z] Z ]$ and the code Jens currently has doesn't allow for the computation $E [ E[Z] Z ] = E[Z] E[Z] = ( E[Z] )^2$. Same comment applies for other higher centralized moments. – user32416 Aug 16 '15 at 1:32
• Might this work: Expand[expect[Expand[(z - ez)^3]] /. ez -> expect[z]]. And, of course, instead of 3 one could put in any positive integer. – JimB Aug 16 '15 at 1:51

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.