Covariance of transformed random variables

I am trying to solve the following:

\begin{align*} &X \sim N(1,1)\\ &\mathrm{cov}(X, X^3) = \text{?} \end{align*}

where $\mathrm{cov}$ is the covariance.

How would you do this in Mathematica?

I have tried

X = NormalDistribution[1, 1]

cov[x_, y_] := Mean[TransformedDistribution[a*b,
{a \[Distributed] x, b \[Distributed] y}]] - Mean[x] Mean[y]

cov[X, TransformedDistribution[a^3, a \[Distributed] X]]


But this doesn't seem to work.

(This is too long for a comment.)

About your comment under 0x4A4D's answer: I think you didn't make it clear enough if your $X$ is a random variable or fixed data generated from some distribution. Usually, we interpret the notation in your question with the former meaning. In that case, $Y=X^3$ merely means a certain relationship between the distributions of two independent random variables, so $\mathbb{E}(XY)\neq\mathbb{E}(X\cdot X^3)$.

• Compare the differences between the following two cases:

xdata = RandomVariate[NormalDistribution[1, 1], 10^6];
x3data = xdata^3;
Covariance[xdata, x3data]


5.96164

xdata = RandomVariate[NormalDistribution[1, 1], 10^6];
x3data = RandomVariate[NormalDistribution[1, 1], 10^6]^3;
Covariance[xdata, x3data]


0.00830199

Mathematica has the ability to deal with multivariable distributions:

multiDist = TransformedDistribution[{x1, x2^3},
{x1 \[Distributed] NormalDistribution[1, 1],
x2 \[Distributed] NormalDistribution[1, 1]}]
Covariance[multiDist, 1, 2]


0

In the latter case:

multiDist = TransformedDistribution[{x, x^3}, x \[Distributed] NormalDistribution[1, 1]];
Covariance[multiDist, 1, 2]


6

Also we can calculate the covariance with the expanded formula:

xyDist = TransformedDistribution[x1 x2^3,
{x1 \[Distributed] NormalDistribution[1, 1],
x2 \[Distributed] NormalDistribution[1, 1]}]
Exy = Expectation[xy, xy \[Distributed] xyDist]
Ex = Expectation[x, x \[Distributed] NormalDistribution[1, 1]]
Ey = Expectation[x^3, x \[Distributed] NormalDistribution[1, 1]]
COVxy = Exy - Ex Ey


0

In the latter case:

xyDist = TransformedDistribution[x x^3, x \[Distributed] NormalDistribution[1, 1]];
Exy = Expectation[xy, xy \[Distributed] xyDist];
Ex = Expectation[x, x \[Distributed] NormalDistribution[1, 1]];
Ey = Expectation[x^3, x \[Distributed] NormalDistribution[1, 1]];
COVxy = Exy - Ex Ey


6.

Please note the difference between the distributions.

• @0x4A4D Thanks for the edit.. I will be more careful about the formatting. – Silvia Jun 5 '13 at 9:37

I guess something like this:

d1 = NormalDistribution[1, 1];

xa = Mean[d1]; xa3 = Mean[TransformedDistribution[u^3, u \[Distributed] d1]];

Mean[TransformedDistribution[(x - xa) (x^3 - xa3), x \[Distributed] d1]]
6

• Here, $Y = X^3, E(XY) = E(X^4) = 10 ; E(X) = 1 ; E(Y) = E(X^3) = 4 ; Cov(X,Y) = E(XY) - E(X)E(Y) = 10 - (1)(4) = 6$, not 0. There is a good chance I made some mistake, but I don't know where. – klaufir Jun 5 '13 at 5:02
• Are you sure? I used the more traditional definition of covariance, as opposed to the "expanded" version... – J. M. is away Jun 5 '13 at 5:04
• I am not sure about which result is correct (although your solution is very convincing), just looking for an explanation about the difference. What bothers me, is that the "expanded" and traditional definition should give the same result. – klaufir Jun 5 '13 at 5:08
• The main reason this yields the wrong answer is because ... although you understand that d2 is the distribution of $X^3$, in Mma syntax, d2 is just the distribution of some slot variable. Then, in TransformedDistribution[(x-xa)(y-ya), {x~d1, y~d2}]] ... the x and the y are seen as different variables ... not as $x$ and $x^3$. But there are more problems: Your APPROACH is to instruct mma to derive the pdf of (x-xa)(y-ya), and thus the mean ...but: (i) you don't need the pdf to perform the expectation and (ii) Mma cannot derive the pdf anyway (try), so it is not even doing what you think. – wolfies Jun 5 '13 at 12:15
• @wolfies, okay, I've removed the assumption of independence... – J. M. is away Jun 5 '13 at 12:18

Given $X$ ~ $N(\mu, \sigma^2)$ with pdf $f(x)$:

$$f=\frac{1}{\sqrt{2 \pi } \sigma } {\text{Exp} \left[-\frac{(x-\mu )^2}{2 \sigma ^2}\right]}; \text{ domain}[f]=\{x,-\infty ,\infty \}\land \{\mu \in \text{Reals},\sigma >0\};$$

Then, using the mathStatica package for Mathematica, the solution is simply:

 Cov[{x, x^3}, f]


$3 \sigma ^2 \left(\mu ^2+\sigma ^2\right)$

In your specific case, with $\mu =1$ and $\sigma^2=1$, the answer is thus 6.