1
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ Can't we argue that $E[y]=0$ from first principles? $\endgroup$
    – David G. Stork
    Feb 27, 2017 at 18:44
  • $\begingroup$ yes, E[y] looks like 0 $\endgroup$
    – Yaroslav Bulatov
    Feb 27, 2017 at 18:58
  • $\begingroup$ 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. $\endgroup$
    – David G. Stork
    Feb 27, 2017 at 19:01
  • $\begingroup$ 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. $\endgroup$
    – mikado
    Feb 27, 2017 at 19:33
  • $\begingroup$ 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. $\endgroup$
    – David G. Stork
    Feb 27, 2017 at 21:29

3 Answers 3

Reset to default
5
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks, that's helpful. That expression seems weird though. If you plug in values of k close to 1 it gives implausibly large results $\endgroup$
    – Yaroslav Bulatov
    Feb 27, 2017 at 22:09
  • $\begingroup$ I think Mathematica had problems solving for both the cases when $k>1$ and $k<1$. The new answer is remarkably simple! $\endgroup$
    – David G. Stork
    Feb 27, 2017 at 22:15
  • $\begingroup$ thanks, that works. Actually it works with multinormal too and seems a bit cleaner, posted that example for posterity $\endgroup$
    – Yaroslav Bulatov
    Feb 27, 2017 at 22:36
  • $\begingroup$ it seems they still haven't made it possible to use subscripts in variable names... $\endgroup$
    – Yaroslav Bulatov
    Feb 27, 2017 at 22:38
  • $\begingroup$ 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). $\endgroup$
    – JimB
    Feb 27, 2017 at 22:45
4
$\begingroup$ 
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}$$

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ I tried your code Assuming[0<k, ...] and got different answers for $k<1$ and $k>1$. $\endgroup$
    – David G. Stork
    Feb 27, 2017 at 23:56
1
$\begingroup$

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]]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.