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.

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

$\endgroup$
  • $\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 '17 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 '17 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 '17 at 22:36
  • $\begingroup$ it seems they still haven't made it possible to use subscripts in variable names... $\endgroup$ – Yaroslav Bulatov Feb 27 '17 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 '17 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}$$

$\endgroup$
  • 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 '17 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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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