In several applications in surveys, it would be helpful to be able to generate a set of $R$ $n$-dimensional variates with the following properties:

  1. Has mean vector $1$
  2. Has a specified variance-covariance matrix ($\Sigma$, a positive semi-definite matrix, with rank $\leq R$)
  3. Is nonnegative

Some First Thoughts about Solutions

The first two requirements are easy to satisfy. We can just draw from a multivariate Normal distribution, as discussed here: Generating data with a given sample covariance matrix

But if we add the additional requirement that all of the variates are nonnegative, then a multivariate Normal distribution won't work. In the applications I have in mind (described below), typically the diagonal of $\Sigma$ has a few entries which are $1$ or even as large as $1.5$, so a multivariate Normal will easily generate negative values.

In this application, it doesn't matter at all what the skew or kurtosis are, and it doesn't matter whether the random variates come from a particular distribution. All that matters is that they're nonnegative, have mean vector $1$, and have the specified covariance matrix (ideally, it would have the exact specified sample covariance matrix, but it would be OK if it just had the specified covariance matrix in expectation).

The multivariate lognormal and Gamma distributions are nonnegative and seem like fairly natural options, except that there are constraints on their precise shape (due to their density functions) which mean that they will often not be able to attain the desired variance-covariance. So these parametric distributions, at least, seem unnecessarily limiting.

In low dimensions, one can generate random variates and then "fix them up" to satisfy constraints, at least approximately. This StackExchange gets at this kind of approach, but it doesn't really cover the constraint of nonnegativity and is really focused on a univariate case.

How to simulate data that satisfy specific constraints such as having specific mean and standard deviation?

Motivating Application

In survey statistics, resampling methods such as the jackknife or bootstrap are typically implemented by generating replicate weights, which are random variates with a specific mean vector and covariance matrix. To be concrete, if one has a dataset of size $n$, then one generates $R$ sets of replicate weights, which can be represented as a matrix of dimension $R \times n$. In general, we want these replicate weights to be nonnegative, since the weights will frequently be used in statistical procedures that require nonnegative weights and usually even for many procedures that allow negative weights, the software implementations won't allow negative weights because they're just unexpected.

This R package vignette describes the generation of replicate weights for the "generalized survey bootstrap", and describes how multivariate normal distributions are used to generate replicate weights that have mean vector $1$ and a specified variance-covariance matrix, but which can sometimes be negative. https://cran.r-project.org/web/packages/svrep/vignettes/bootstrap-replicates.html#forming-adjustment-factors It also describes a rescaling adjustments that makes the random variates nonnegative and with mean vector $1$, but unfortunately that rescaling adjustment increases the variance-covariance by a constant (which is a problem).

  • $\begingroup$ Thanks- I added that tag to the question. $\endgroup$
    – bschneidr
    Mar 3 at 21:45
  • $\begingroup$ What I meant is that you should investigate copulas as a direct manner of imposing a given covariance matrix on rv's with given margins. $\endgroup$
    – Xi'an
    Mar 4 at 8:55
  • 1
    $\begingroup$ Copulae could be a real pain (after having to use them for multivariate Gamma), however here's one good beginners' tutorial: twiecki.io/blog/2018/05/03/copulas $\endgroup$
    – Spätzle
    Mar 14 at 6:49
  • 1
    $\begingroup$ Another option you should check out is the folded normal distribution, more specifically the multivariate folded normal: jstor.org/stable/42003783 $\endgroup$
    – Spätzle
    Mar 14 at 6:50
  • 1
    $\begingroup$ What were typical/interesting sizes for $R$ and $n$? I ask because attempts at direct computational solutions may not work well for those being large. $\endgroup$
    – g g
    Mar 16 at 11:00

1 Answer 1


First notice that the requirements of non-negativity and mean equals 1 cannot be simultaneously satisfied for an arbitrary covariance matrix: if $x_i$ is the $i$-th variate ($i=1...n$) with $E[x_i]=1$, then

$$Cov(x_i,x_j)=E[(x_i-1)(x_j-1)]=E[x_ix_j]-1 \ge -1$$

since all the $x_i$'s are nonnegative and therefore $E[x_ix_j] \ge 0$. This implies that all the entries of the covariance matrix $\Sigma_{ij}$ has to be larger than -1, which is a non-trivial constraint. (Note that the same argument applies to the sample covariance. Also note that it does not imply that if $\Sigma_{ij} \ge -1$ then necessarily a solution exists; there might be other constraints).

If you are not interested in the generating distribution but just want to find some matrix $X$ which has the specified properties (assuming it exists), a straightforward approach might be to try to minimize a sum of squared errors using a gradient descent algorithm. To ensure nonnegativity the elements of $X$ can be parametrized e.g. with $X_{in} = e^{W_{in}}$, and you can take advantage of auto-differentiation libraries to calculate the gradients for you.

Here is for example an implementation in PyTorch, which finds a (approximate) solution for $R=20, n=3$ in about 2 seconds:

import torch

R, n = 20, 3

#generate n-by-n random covariance matrix with entries larger than -1
while S.min() < -1:
    A = torch.randn(n,n)/2
    S = A.matmul(A.T)

def loss(W,S):
    X = torch.exp(W)
    sse = (torch.cov(X.T,correction = 0) - S).square().sum() + (X.mean(0)-torch.ones(n)).square().sum()
    return sse

#initialize W with random entries
W = torch.randn(R,n,requires_grad=True)
optimizer = torch.optim.Adam([W],lr=0.05)

for iter in range(1000):
    L = loss(W,S)

print(f'sse = {L}')
print('specified covariance:')
print('sample covariance:')
print(torch.cov(X.T,correction = 0))


sse = 0.007093742955476046
specified covariance:
tensor([[ 1.9820, -0.2131, -0.5573],
        [-0.2131,  0.8370, -0.0203],
        [-0.5573, -0.0203,  0.1815]])
sample covariance:
tensor([[ 1.9881, -0.2122, -0.5349],
        [-0.2122,  0.8371, -0.0169],
        [-0.5349, -0.0169,  0.2584]], grad_fn=<SqueezeBackward0>)
tensor([1.0026, 1.0003, 1.0093], grad_fn=<MeanBackward1>)
  • $\begingroup$ Thanks, I appreciate your comment about the constraints the nonnegativity and mean requirements imply for the covariance matrix. For the gradient descent algorithm, I wonder how feasible this would be in real applications, where you might have $n=10,000$ and $50 \leq R \leq 1,000$. That's a lot of parameters to try to directly optimize. $\endgroup$
    – bschneidr
    Mar 16 at 15:21
  • $\begingroup$ That number of parameters ($R \times n$) is still quite small compared to typical deep learning models which use the same optimization algorithm. I am curious however, from where do you get a $10,000$-by-$10,000$ matrix with rank $50$ ? $\endgroup$
    – J. Delaney
    Mar 16 at 17:42
  • 1
    $\begingroup$ I see that, but you want the resulting covariance matrix to be equal to $\Sigma$ which is pre-specified (and needs to have rank $R$). Where does this $\Sigma$ come from ? $\endgroup$
    – J. Delaney
    Mar 16 at 20:45
  • 1
    $\begingroup$ This sounds like you are not interested in a solution for generic $\Sigma$ but in very specific highly structured examples. Since in general a solution may not exist as @J.Delaney has demonstrated, it might be much smarter to start with one of the special cases and take it from there? $\endgroup$
    – g g
    Mar 17 at 8:36
  • 1
    $\begingroup$ @bschneidr $\Sigma$ in your example has rank $n-1$. Since you want $R \ll n$, you clearly can't get this sample covariance (or even a reasonable low rank approximation, since all non-zeros eigenvalues of $\Sigma$ are close to 1) $\endgroup$
    – J. Delaney
    Mar 17 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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