Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
up vote 8 down vote accepted

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


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


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]


In the latter case:

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


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


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


Please note the difference between the distributions.

share|improve this answer
@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]]
share|improve this answer
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. 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. 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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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