Generating data with a given sample covariance matrix

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?

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.)

• +1. The other day, I needed to generate some data with a given sample covariance matrix, didn't know how to do it, and for some reason it took me a lot of time to find your answer. To increase the visibility of this thread and to illustrate your suggestions, I posted another answer here with some Matlab code. – amoeba Feb 3 '15 at 23:21
• @amoeba I wonder if there's a possibility of putting one of the search terms you used which are not already present here into the question tags (or possibly inserting several in a small edit to the text of the question, which should still help it be found). I'm now wondering whether I should do the same thing in R ... but then does it go better in my answer, or as an addition to yours? – Glen_b Feb 3 '15 at 23:27
• I have already taken the liberty to edit the question, and also tried to formulate my answer such that it includes as many keywords as possible. Hope this will help. I was surprised, by the way, that this simple tip (whitening the generated data before transforming to the required covariance) was so difficult to google; could not find anything (on CV or elsewhere), until I finally found your answer. – amoeba Feb 3 '15 at 23:32
• @amoeba Oh, okay, thanks. Yeah, actually, I can't say I recall ever seeing it mentioned anywhere for the multivariate case (no doubt it has been, since it's a fairly obvious idea, especially if you've already thought of it for the univariate case, or have already seen it in the univariate case). – Glen_b Feb 3 '15 at 23:36
• @Glen_b As you say, the resulting distribution of these "cleansed" samples cannot be normal. Do you have any idea what the resulting distribution could be? Or maybe whether it is equal/not equal to the conditional distribution $Z\mid \bar{z}=\mu, Cov(z)=\Sigma$ – g g Apr 30 at 15:05

@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
• +1. Somehow this question is in various forms in various locations in CV. Is there a way to proceed if we are aware that the multivariate distribution is non-gaussian ? – rgk Apr 12 at 20:44
• If you know the form of the multivariate distribution you want it to look like, perhaps in some cases. – Glen_b May 2 at 0:35