Generating data with a given sample covariance matrix

Question

Given a covariance matrix $\boldsymbol \Sigma_s$, how to generate data such that it would have the sample covariance matrix $\hat{\boldsymbol \Sigma} = \boldsymbol \Sigma_s$?

---

More generally: we are often interested in generating data from a density $ f(x \vert \boldsymbol\theta) $, with data $x$ given some parameter vector $\boldsymbol\theta$. This results in a sample, from which we may then again estimate a value $\boldsymbol{\hat\theta}$. What I am interested in is the reverse problem: What if we are given a set of parameters $\boldsymbol\theta_{s}$, and we would like to generate a sample $x$ such, that $ \boldsymbol{\hat\theta} = \boldsymbol\theta_{s}$.

Is this a known problem? Is such a method useful? Are algorithms available?

Answer 1

There are two different typical situations for these kind of problems:

i) you want to generate a sample from a given distribution whose population characteristics match the ones specified (but due to sampling variation, you don't have the sample characteristics exactly matching).

ii) you want to generate a sample whose sample characteristics match the ones specified (but, due to the constraints of exactly matching sample quantities to a prespecified set of values, don't really come from the distribution you want).

You want the second case -- but you get it by following the same approach as the first case, with an extra standardization step.

So for multivariate normals, either can be done in a fairly straightforward manner:

With first case you could use random normals without the population structure (such as iid standard normal which have expectation 0 and identity covariance matrix) and then impose it - transform to get the covariance matrix and mean you want. If $\mu$ and $\Sigma$ are the population mean and covariance you need and $z$ are iid standard normal, you calculate $y=Lz+\mu$, for some $L$ where $LL'=\Sigma$ (e.g. a suitable $L$ could be obtained via Cholesky decomposition). Then $y$ has the desired population characteristics.

With the second, you have to first transform your random normals to remove even the random variation away from the zero mean and identity covariance (making the sample mean zero and sample covariance $I_n$), then proceed as before. But that initial step of removing the sample deviation from exact mean $0$, variance $I$ interferes with the distribution. (In small samples it can be quite severe.)

This can be done by subtracting the sample mean of $z$ ($z^*=z-\bar z$) and calculating the Cholesky decomposition of $z^*$. If $L^*$ is the left Cholesky factor, then $z^{(0)}=(L^*)^{-1}z^*$ should have sample mean 0 and identity sample covariance. You can then calculate $y=Lz^{(0)}+\mu$ and have a sample with the desired sample moments. (Depending on how your sample quantities are defined, there may be an extra small fiddle involved with multiplying/dividing by factors like $\sqrt{\frac{n-1}{n}}$, but it's easy enough to identify that need.)

Answer 2

@Glen_b gave a good answer (+1), which I want to illustrate with some code.

How to generate $n$ samples from a $d$-dimensional multivariate Gaussian distribution with a given covariance matrix $\boldsymbol \Sigma$? This is easy to do by generating samples from a standard Gaussian and multiplying them by a square root of the covariance matrix, e.g. by $\mathrm{chol}(\boldsymbol \Sigma)$. This is covered in many threads on CV, e.g. here: How can I generate data with a prespecified correlation matrix? Here is a simple Matlab implementation:

n = 100;
d = 2;
Sigma = [ 1    0.7  ; ...
          0.7   1   ];
rng(42)
X = randn(n, d) * chol(Sigma);

The sample covariance matrix of the resulting data will of course not be exactly $\boldsymbol \Sigma$; e.g. in the example above cov(X) returns

1.0690    0.7296
0.7296    1.0720

How to generate data with a pre-specified sample correlation or covariance matrix?

As @Glen_b wrote, after generating data from a standard Gaussian, center, whiten, and standardize it, so that it has sample covariance matrix $\mathbf I$; only then multiply it with $\mathrm{chol}(\boldsymbol \Sigma)$.

Here is the continuation of my Matlab example:

X = randn(n, d);
X = bsxfun(@minus, X, mean(X));
X = X * inv(chol(cov(X)));
X = X * chol(Sigma);

Now cov(X), as required, returns

1.0000    0.7000
0.7000    1.0000