Generate nonnegative variates with mean 1 and specified variance-covariance

Question

Problem

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

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
S=torch.Tensor([-2])
while S.min() < -1:
    A = torch.randn(n,n)/2
    S = A.matmul(A.T)


def loss(W,S):
    [R,n]=W.shape
    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):
    optimizer.zero_grad()
    L = loss(W,S)
    L.backward()
    optimizer.step()

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

output:

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>)
mean:
tensor([1.0026, 1.0003, 1.0093], grad_fn=<MeanBackward1>)