# Finding mean and covariance of an arbitrary joint pdf

Let's say that I have the following joint pdf, defined on all of $\mathbb{R}^2$:

$f(x,y) = \frac{1}{54.4366}\exp(\frac{-x^2}{200}-\frac{1}{2}(y+0.05x^{2}-5)^2)$

How would I go about calculating the mean and covariance matrix associated with this pdf?

• You sure that it is $0.05x^2$ in the exponent and not $0.05x$? The problem is relatively straightforward if it is $0.05x$. That $0.05x^2$ makes it much messier. – Dilip Sarwate Aug 10 '17 at 3:23
• Unfortunately it is $0.05x^2$ in the exponent. Despite this, how would I go about finding the mean and covariance? – tattybojangler Aug 10 '17 at 4:07
• I believe that 54.4366 should be $20\pi \simeq 62.8319$ for $f(x,y)$ to be a joint pdf. – JimB Aug 10 '17 at 4:49

As Jim Baldwin mentioned in the comments, your expression for the joint PDF doesn't integrate to 1, so it isn't a valid probability distribution. I'll assume the normalizing constant is $20 \pi$ instead of 54.366, which gives a valid distribution.

Joint distribution

$$f_{XY}(x, y) = \frac{1}{20 \pi} \exp \left ( -\frac{x^2}{200} - \frac{1}{2} (y + 0.05 x^2-5)^2 \right )$$

This distribution also has a convenient form when expressed as the product of the marginal distribution of $X$ and the conditional distribution of $Y$ given $X$:

$$f_{XY}(x, y) = f_X(x) \ f_{Y \mid X}(x, y)$$

$f_{X}$ is Gaussian with mean 0 and variance 100 (as shown below), and $f_{Y \mid X}$ is Gaussian with mean $5-\frac{x^2}{20}$ and variance 1 (which can be shown by dividing $f_{XY}$ by $f_{X}$). $f_{XY}$ has a narrow, parabolic shape, so the mean and covariance matrix won't do a great job of summarizing it. But, they can be computed as follows.

Mean

The mean is $[\mu_x, \mu_y]$, where $\mu_x$ and $\mu_y$ are the expected values of the marginal distributions.

The marginal distribution of $X$ can be found by integrating out $Y$:

$$f_X(x) = \int_{=\infty}^\infty f_{XY}(x, y) dy = \frac{1}{10 \sqrt{2 \pi}} \exp \left ( -\frac{x^2}{200} \right )$$

This is a Gaussian with mean 0 and variance 100, so $\mu_X = 0$.

Similarly, $\mu_Y$ can be found by integrating out $X$, then taking the expected value:

$$f_Y(y) = \int_{-\infty}^\infty f_{XY}(x, y) dx$$

$$\mu_y = \int_{-\infty}^\infty y \ f_Y(y) dy$$

I couldn't find closed form expressions, but the integrals can be computed numerically using standard software. This typically requires specifying finite integration bounds. You can use bounds at which the density tapers to negligibly small values (check that the distribution integrates numerically to ~1 over your chosen bounds).

Covariance

The covariance matrix is:

$$\left [ \begin{array}{cc} \sigma_{X}^2 & \sigma_{XY}^2 \\ \sigma_{XY}^2 & \sigma_{Y}^2 \\ \end{array} \right ]$$

where $\sigma_X^2$ and $\sigma_Y^2$ are the variances of $X$ and $Y$, and $\sigma_{XY}^2$ is the covariance of $X$ and $Y$. As shown above, $X$ has variance 100.

$\sigma_Y^2$ can be computed (again using numerical integration) as:

$$\sigma_Y^2 = \int_{-\infty}^\infty (y-\mu_y)^2 f_Y(y) dy$$

For $\sigma_{XY}^2$, we can use the following relationship:

$$\text{cov}(X, Y) = E[(X-E[X])(Y-E[Y])] = E[XY] - E[X]E[Y]$$

Therefore:

$$\sigma_{XY}^2 = E[XY] - \mu_X \mu_Y$$

All that remains is to calculate $E[XY]$ which can be done in closed form:

$$E[XY] = \int_{-\infty}^\infty \int_{-\infty}^\infty x y \ f_{XY}(x, y) dx dy = 0$$

Therefore, $\sigma_{XY}^2 = 0$. This makes sense, looking at the symmetry of the joint distribution about the y axis.

• Thank you for your detailed explanation, this makes perfect sense! – tattybojangler Aug 10 '17 at 19:25

Just to add to @user20160 's answer, here is some Mathematica code to perform the calculations. (And it does not appear that the mean and variance of $y$ can be given in a closed form hence the numerical integration below.)

(* Check that things integrate to 1 *)
Integrate[Exp[-x^2/200 - (1/2) (y + x^2/20 - 5)^2], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]/(20 Pi)
(* 1 *)

(* Marginal distribution of x *)
fx = Integrate[f[x, y], {y, -Infinity, Infinity}]
(* \[ExponentialE]^(-(x^2/200))/(10 Sqrt[2 \[Pi]]) *)

(* Mean of x *)
xMean = Integrate[x Integrate[f[x, y], {y, -Infinity, Infinity}], {x, -Infinity, Infinity}]
(* 0 *)

(* Variance of x *)
xVariance = Integrate[x^2 Integrate[f[x, y], {y, -Infinity, Infinity}], {x, -Infinity, Infinity}] - xMean^2
(* 100 *)

(* Marginal density of y is fy1 if y < 49/10 and fy2 if y > 49/10 *)
fy1 = Integrate[E^(-(x^2/200) - 1/2 (-5 + x^2/20 + y)^2)/(20 \[Pi]),
{x, -Infinity, Infinity}, Assumptions -> y < 49/10]
fy2 = Integrate[E^(-(x^2/200) - 1/2 (-5 + x^2/20 + y)^2)/(20 \[Pi]),
{x, -Infinity, Infinity}, Assumptions -> y >= 49/10]
(* Density of y at y = 49/10 *)
Limit[fy2, y -> 49/10]
(* Gamma[1/4]/(2 2^(3/4) Sqrt \[ExponentialE]^(1/200) \[Pi]) *)

(* Mean of y *)
yMean = NIntegrate[y fy1, {y, -30, 49/10}] + NIntegrate[y fy2, {y, 49/10, 15}]
(* 0.3207576948514401 *)

(* Variance of y *)
yVariance = NIntegrate[y^2 fy1, {y, -30, 49/10}] + NIntegrate[y^2  fy2, {y, 49/10, 15}] - yMean^2
(* 37.65761182719062 *)

(* Covariance of x and y *)
cov = Integrate[x y f[x, y], {y, -Infinity, Infinity}, {x, -Infinity, Infinity}]
(* 0 *) where $I_n (x)$ and $K_n (x)$ are the modified Bessel functions of the first and second kind, respectively. • Thank you for your help! I admit it is a rather messy distribution to be working with! – tattybojangler Aug 10 '17 at 19:26