# Expectation of function of independent random variables

I see that Mathematica added some random variable features in the last 5 years, are there some nice tutorials for finding expectations of functions of random variables?

In particular, I have x1 and x2 are sampled independently from $\mathcal{N}(0,\Sigma)$ where

$$\Sigma = \left( \begin{array}{cc} 1 & 0 \\ 0 & k \\ \end{array} \right)$$

Then I have a random variable $y$ defined as

$$y = \frac{<x1, x2>}{\|x1\|\|x2\|}$$

and I need to find $E[y]$ and $E[y^2]$. Basically it's finding how the angle between two random vectors depends on their covariance matrix.

• Can't we argue that $E[y]=0$ from first principles? Feb 27, 2017 at 18:44
• yes, E[y] looks like 0 Feb 27, 2017 at 18:58
• For every condition where there is a positive value of ${{\bf x}_1 \cdot {\bf x}_2 \over \| {\bf x}_1 \| \ \| {\bf x}_2 \|}$, there is a condition with the corresponding negative value. Feb 27, 2017 at 19:01
• You might be better asking on the Mathematics site, as this is not really specific to Mathematica. I think that this might have a solution that is well known to the right audience. Feb 27, 2017 at 19:33
• Here is how you'd express this in Mathematica: Expectation[x1.x2/(Norm[x1] Norm[x2]), {x1, x2} \[Distributed] MultinormalDistribution[{0, 0}, {{1, 0}, {0, k}}]] and likewise for the $E[y^2]$ case. Unfortunately, Mathematica cannot solve these in this form. Feb 27, 2017 at 21:29

For some reason, Mathematica can handle the problem when represented as one-dimensional distributions, which is allowed given your covariance matrix is diagonal.

$E[y]$

Assuming[k > 0,
Expectation[
(x1x x2x + x1y x2y)/((x1x^2 + x1y^2) (x2x^2 + x2y^2))^(1/2),
{x1x \[Distributed] NormalDistribution[0, 1],
x1y \[Distributed] NormalDistribution[0, k],
x2x \[Distributed] NormalDistribution[0, 1],
x2y \[Distributed] NormalDistribution[0, k]}]
]


$0$ (as it should, based on symmetry arguments)

$E[y^2]$

Mathematica has difficulty solving for the general case when $k>1$ and $k<1$, but can do each individually (and they give the same answer):

Assuming[k > 1,
Expectation[(x1x x2x + x1y x2y)^2/((x1x^2 + x1y^2) (x2x^2 + x2y^2)),
{x1x \[Distributed] NormalDistribution[0, 1],
x1y \[Distributed] NormalDistribution[0, Sqrt[k]],
x2x \[Distributed] NormalDistribution[0, 1],
x2y \[Distributed] NormalDistribution[0, Sqrt[k]]}]]


$\frac{k+1}{(\sqrt{k}+1)^2}$

and

Assuming[0 < k < 1,
Expectation[(x1x x2x + x1y x2y)^2/((x1x^2 + x1y^2) (x2x^2 + x2y^2)),
{x1x \[Distributed] NormalDistribution[0, 1],
x1y \[Distributed] NormalDistribution[0, Sqrt[k]],
x2x \[Distributed] NormalDistribution[0, 1],
x2y \[Distributed] NormalDistribution[0, Sqrt[k]]}]]


$\frac{k+1}{(\sqrt{k}+1)^2}$.

An explicit check for the case $k=1$:

Expectation[(x1x x2x + x1y x2y)^2/((x1x^2 + x1y^2) (x2x^2 + x2y^2)),
{x1x \[Distributed] NormalDistribution[0, 1],
x1y \[Distributed] NormalDistribution[0, 1],
x2x \[Distributed] NormalDistribution[0, 1],
x2y \[Distributed] NormalDistribution[0, 1]}]


${1 \over 2}$, which agrees with the general case.

• Thanks, that's helpful. That expression seems weird though. If you plug in values of k close to 1 it gives implausibly large results Feb 27, 2017 at 22:09
• I think Mathematica had problems solving for both the cases when $k>1$ and $k<1$. The new answer is remarkably simple! Feb 27, 2017 at 22:15
• thanks, that works. Actually it works with multinormal too and seems a bit cleaner, posted that example for posterity Feb 27, 2017 at 22:36
• it seems they still haven't made it possible to use subscripts in variable names... Feb 27, 2017 at 22:38
• I think you need to use Sqrt[k] in place of k to match the original question (which has the variance being k rather than the standard deviation being k).
– JimB
Feb 27, 2017 at 22:45
normal := MultinormalDistribution[{0, 0}, ( {{1, 0}, {0, k} } )];
x := {x1, x2};
y := {y1, y2};
vars := {x \[Distributed] normal, y \[Distributed] normal};
Assuming[0 < k < 1, Expectation[(x.y/(Norm[x] Norm[y]))^2, vars]]


Result

$$\frac{k+1}{\left(\sqrt{k}+1\right)^2}$$

• I tried your code Assuming[0<k, ...] and got different answers for $k<1$ and $k>1$. Feb 27, 2017 at 23:56

Just for fun:

f[k_] := MultinormalDistribution[{0, 0}, ({{1, 0}, {0, k}})];
rv[k_, n_] :=
Mean[(#1.#2)^2/(#1.#1 #2.#2) & @@@ RandomVariate[f[k], {n, 2}]]
vis[n_] :=
Show[Plot[(1 + k)/(1 + Sqrt[k])^2, {k, 0, 1}, PlotStyle -> Red,
PlotRange -> All],
ListPlot[Table[{j, rv[j, n]}, {j, 0.01, 1, 0.01}]],
PlotRange -> {0, 1}, AxesOrigin -> {0, 0}, Frame -> True,
GridLines -> {None, {1/2}}, PlotLabel -> Row[{"n= ", n}],
ImageSize -> 200]
Grid[Partition[vis /@ {100, 1000, 10000, 100000}, 2]]