# 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. Commented Jun 5, 2013 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. Commented Jun 5, 2013 at 5:02
• Are you sure? I used the more traditional definition of covariance, as opposed to the "expanded" version... Commented Jun 5, 2013 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. Commented Jun 5, 2013 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. Commented Jun 5, 2013 at 12:15
• @wolfies, okay, I've removed the assumption of independence... Commented Jun 5, 2013 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.