import matplotlib.pyplot as plt
import numpy as np
np.random.seed(42)
Generate points on a circle.
angles = np.random.random(10000) * 2.0 * np.pi
radius = 1.0
x = np.cos(angles) * radius
y = np.sin(angles) * radius
plt.scatter(x, y, s=1)
plt.gca().set_aspect('equal')
Compute the covariance matrix.
np.cov(np.array([x, y]))
array([[ 0.50309341, 0.00297389], [ 0.00297389, 0.49683412]])
$X$, $Y$, and $Z$ are distributed uniformly in the ranges $[0, \alpha], [0, \beta], [0, \gamma]$. This means they have variance $\alpha^2 / 12, \beta^2 / 12, \gamma^2 / 12$ which give us the diagonal terms in the covariance matrix. The other terms are calculated as follows (shown for $X$ and $Y$, but similar for other pairs of random variables):
$$ cov(X, Y) = \langle (X - \overline{X}) (Y - \overline{Y}) \rangle = \langle (\alpha \lambda - \alpha / 2) (\beta \lambda - \beta / 2) \rangle = \langle \alpha (\lambda - 1 / 2) \beta (\lambda - 1 / 2) \rangle = \alpha \beta \langle (\lambda - 1 / 2) (\lambda - 1 / 2) \rangle = \alpha \beta * var(\lambda) = \frac{\alpha \beta}{12} $$Thus the covariance matrix is:
$$ C = \: \begin{bmatrix} \alpha^2 & \alpha \beta & \alpha \gamma\\ \beta \alpha & \beta^2 & \beta \gamma \\ \gamma \alpha & \gamma \beta & \gamma^2 \end{bmatrix} $$The correlation matrix is calculated as:
$$ r_{ij} = \frac{C_{i j}}{\sqrt{C_{i i} C_{j j}}} $$which results in:
$$ r = \: \begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \\ \end{bmatrix} $$For $\alpha = 1, \beta = 2, \gamma = 3$, compute the expected covariance and correlation matrices as described. Then generate points with these parameters, compute the actual matrices, and verify they match expectations.
def cov_and_corr(a, b, c, num_samples=1000000):
"""Generate points as the homework specifies. Return covariance and correlation."""
abc = np.tile(np.array([[a, b, c]]), (num_samples, 1))
lambdas = np.tile(np.array([np.random.random(num_samples)]).T, (1, 3))
xyz = abc * lambdas
cov = np.cov(xyz, rowvar=False)
corr = np.corrcoef(xyz, rowvar=False)
return cov, corr
a, b, c = 1, 2, 3
abc = np.array([[a, b, c]])
cov_expected = np.dot(abc.T, abc) / 12.0
corr_expected = np.ones((3, 3), dtype=float)
cov, corr = cov_and_corr(a, b, c)
print 'cov:\n', cov
print 'corr:\n', corr
cov: [[ 0.08329172 0.16658343 0.24987515] [ 0.16658343 0.33316686 0.49975029] [ 0.24987515 0.49975029 0.74962544]] corr: [[ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]]
print np.allclose(cov, cov_expected, rtol=0.01)
print np.allclose(corr, corr_expected, rtol=0.01)
True True