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

Question

This question is motivated by my question on meta-analysis. But I imagine that it would also be useful in teaching contexts where you want to create a dataset that exactly mirrors an existing published dataset.

I know how to generate random data from a given distribution. So for example, if I read about the results of a study that had:

• a mean of 102,
• a standard deviation of 5.2 , and
• a sample size of 72.

I could generate similar data using rnorm in R. For example,

set.seed(1234)
x <- rnorm(n=72, mean=102, sd=5.2)

Of course the mean and SD would not be exactly equal to 102 and 5.2 respectively:

round(c(n=length(x), mean=mean(x), sd=sd(x)), 2)
##     n   mean     sd 
## 72.00 100.58   5.25 

In general I'm interested in how to simulate data that satisfies a set of constraints. In the above case, the constaints are sample size, mean, and standard deviation. In other cases, there might be additional constraints. For example,

• a minimum and a maximum in either the data or the underlying variable might be known.
• the variable might be known to take on only integer values or only non-negative values.
• the data might include multiple variables with known inter-correlations.

Questions

• In general, how can I simulate data that exactly satisfies a set of constraints?
• Are there articles written about this? Are there any programs in R that do this?
• For the sake of example, how could and should I simulate a variable so that it has a specific mean and sd?

Answer 1

The general technique is the 'Rejection Method', where you just reject results that don't meet your constraints. Unless you have some sort of guidance (like MCMC), then you could be generating a lot of cases (depending on your scenario) which are rejected!

Where you're looking for something like a mean and standard deviation and you can create a distance metric of some kind to say how far you are away from your goal, you can use optimisation to search for the input variables which give you the desired output values.

As an ugly example where we will look for a random uniform vector with length 100 which has mean=0 and standard deviation=1.

# simplistic optimisation example
# I am looking for a mean of zero and a standard deviation of one
# but starting from a plain uniform(0,1) distribution :-)
# create a function to optimise
fun <- function(xvec, N=100) {
  xmin <- xvec[1]
  xmax <- xvec[2]
  x <- runif(N, xmin, xmax)
  xdist <- (mean(x) - 0)^2 + (sd(x) - 1)^2
  xdist
}
xr <- optim(c(0,1), fun)

# now lets test those results
X <- runif(100, xr$par[1], xr$par[2])
mean(X) # approx 0
sd(X)   # approx 1

Answer 2

In general, to make your sample mean and variance exactly equal to a pre-specified value, you can appropriately shift and scale the variable. Specifically, if $X_1, X_2, ..., X_n$ is a sample, then the new variables

$$ Z_i = \sqrt{c_{1}} \left( \frac{X_i-\overline{X}}{s_{X}} \right) + c_{2} $$

where $\overline{X} = \frac{1}{n} \sum_{i=1}^{n} X_i$ is the sample mean and $ s^{2}_{X} = \frac{1}{n-1} \sum_{i=1}^{n} (X_i - \overline{X})^2$ is the sample variance are such that the sample mean of the $Z_{i}$'s is exactly $c_2$ and their sample variance is exactly $c_1$. A similarly constructed example can restrict the range -

$$ B_i = a + (b-a) \left( \frac{ X_i - \min (\{X_1, ..., X_n\}) }{\max (\{X_1, ..., X_n\}) - \min (\{X_1, ..., X_n\}) } \right) $$

will produce a data set $B_1, ..., B_n$ that is restricted to the interval $(a,b)$.

Note: These types of shifting/scaling will, in general, change the distributional family of the data, even if the original data comes from a location-scale family.

Within the context of the normal distribution the mvrnorm function in R allows you to simulate normal (or multivariate normal) data with a pre-specified sample mean/covariance by setting empirical=TRUE. Specifically, this function simulates data from the conditional distribution of a normally distributed variable, given the sample mean and (co)variance is equal to a pre-specified value. Note that the resulting marginal distributions are not normal, as pointed out by @whuber in a comment to the main question.

Here is a simple univariate example where the sample mean (from a sample of $n=4$) is constrained to be 0 and the sample standard deviation is 1. We can see that the first element is far more similar to a uniform distribution than a normal distribution:

library(MASS)
 z = rep(0,10000)
for(i in 1:10000)
{
    x = mvrnorm(n = 4, rep(0,1), 1, tol = 1e-6, empirical = TRUE)
    z[i] = x[1]
}
hist(z, col="blue")

$ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $

Answer 3

Are there any programs in R that do this?

The Runuran R package contains many methods for generating random variates. It uses C libraries from the UNU.RAN (Universal Non-Uniform RAndom Number generator) project. My own knowledge of the field of random variate generation is limited, but the Runuran vignette provides a nice overview. Below are the available methods in the Runuran package, taken from the vignette:

Continuous distributions:

• Adaptive Rejection Sampling
• Inverse Transformed Density Rejection
• Polynomial Interpolation of Inverse CDF
• Simple Ratio-of-Uniforms Method
• Transformed Density Rejection

Discrete distributions:

• Discrete Automatic Rejection Inversion
• Alias-Urn Method
• Guide-Table Method for Discrete Inversion

Multivariate distributions:

• Hit-and-Run algorithm with Ratio-of-Uniforms Method
• Multivariate Naive Ratio-of-Uniforms Method

Example:

For a quick example, suppose you wanted to generate a Normal distribution bounded between 0 and 100:

require("Runuran")

## Normal distribution bounded between 0 and 100
d1 <- urnorm(n = 1000, mean = 50, sd = 25, lb = 0, ub = 100)

summary(d1)
sd(d1)
hist(d1)

The urnorm() function is a convenient wrapper function. I believe that behind the scenes it uses the Polynomial Interpolation of Inverse CDF method but am not sure. For something more complex, say, a discrete Normal distribution bounded between 0 and 100:

require("Runuran")

## Discrete normal distribution bounded between 0 and 100
# Create UNU.RAN discrete distribution object
discrete <- unuran.discr.new(pv = dnorm(0:100, mean = 50, sd = 25), lb = 0, ub = 100)

# Create UNU.RAN object using the Guide-Table Method for Discrete Inversion
unr <- unuran.new(distr = discrete, method = "dgt")

# Generate random variates from the UNU.RAN object
d2 <- ur(unr = unr, n = 1000)

summary(d2)
sd(d2)
head(d2)
hist(d2)