SAS/IML software is often used for sampling and simulation studies. For simulating data from univariate distributions, the RANDSEED and RANDGEN subroutines suffice to sample from a wide range of distributions. (I use the terms "sampling from a distribution" and "simulating data from a distribution" interchangeably.)

For multivariate simulations, the IMLMLIB library contains a series of modules for sampling from multivariate distributions.

Included in this module library is the RANDNORMAL module which enables you to sample from a multivariate normal distribution with a given mean and covariance structure. You need to specify three parameters:

**N**: the number of observations that you want to sample**Mean**: the mean of the distribution (This is the "population mean.")**Cov**: the covariance structure of the distribution. (This is the "population covariance.") The covariance structure is specified by using a matrix whose (*i*,*j*)th element specifies the covariance between the*i*th and*j*th variables. If you want to sample from variables that are uncorrelated, use a diagonal covariance matrix. To sample from identical uncorrelated variables, use an identity matrix.

### Sample from the Multivariate Normal Distribution

The following statements generate 1,000 random observations from a multivariate normal distribution with a specified mean and covariance structure. The first five observations are displayed.

proc iml; /** specify the mean and covariance of the population **/ Mean = {1, 2, 3}; Cov = {3 2 1, 2 4 0, 1 0 5}; NumSamples = 1000; call randseed(4321); /** set seed for the RandNormal module **/ X = RandNormal(NumSamples, Mean, Cov); print (X[1:5,])[label="First 5 Obs: MV Normal"];

How can you verify that the computation is correct? That's an interesting question and a complete answer will have to wait for a future post. However, you can convince yourself that the observations probably came from the specified distribution by computing the sample mean and the sample covariance matrix. These sample statistics should be close to the parameters of the population

### Computing the Mean of the Sample

The following statements compute the mean vector of the sample:

/** check the sample mean **/ SampleMean = X[:, ]; print SampleMean;

The sample mean is, indeed, close to the mean of the population. You should expect the sample mean to approach the population mean as the number of samples gets large.
(Recall that, in the univariate case, the sample mean is distributed as σ^{2}/*N*, where σ^{2} is the population variance. The same idea holds for the multivariate mean vector.)

### Computing the Covariance of the Sample

You can compute the sample covariance matrix by using the following statements:

/** check the sample covariance **/ C = X - SampleMean; SampleCov = C` * C / (NumSamples-1); print SampleCov;

For readers who have upgraded to SAS/IML 9.22, the following statement is an alternative way to compute the covariance matrix:

/** SAS/IML 9.22 and beyond: call the COV function **/ SampleCov2 = cov(X);

The sample covariance matrix is also close to the covariance parameters for the population. The sample covariance matrix should approach the population covariance as the number of samples gets large. (The distribution of the sample covariance matrix is called the Wishart distribution.)

### Sampling When You Have a Correlation Matrix

A correlation matrix is just a covariance matrix for data that have been standardized. Consequently, it is acceptable to use a correlation matrix as the third argument to the RANDNORMAL function. Suppose that S is a covariance matrix and R is the corresponding correlation matrix. If you use R as the third argument to the RANDNORMAL function, the resulting sample will have the same general shape as if you had used S, but the shape will be linearly scaled in certain directions.

If you have correlation matrix instead of covariance, and you also have a set of variances, you can transform the correlation matrix into a covariance matrix.

## 4 Comments

How could I sampling from the bivariate exponetial distribution and bivariate weibull distribution?

That's a difficult question, and there are several parameterizations for those distributions. I recommend the book by Mark Johnson (1987)

Multivariate Statistical Simulation. In it, he shows how to simulate from the Morganstern distribution with correlated exponential marginals. Unfortunately, it only models weakly correlated variables. For the bivariate Weibull, start with X1 and X2 that are bivarate exponential. Then define Y1 = X1^{1/alpha) and Y2 = X2^{1/alpha). Y1 and Y2 are bivaraite Weibull and correlated, although I'm not sure how to express the correlation analytically in terms of the corr(X1,X2).What is the most efficient way to create multiple data sets (in longitudinal format) from the multivariate normal distribution? This is especially useful for running simulations.

Let's say you want to simulate k MVN samples, where each sample contains N observations. You want an ID variable with values 1-k that identifies each sample, and variables x1-x_p that contain the MVN data.

## 6 Trackbacks

[...] are different for each column of x. As I've discussed previously, when you want correlated vectors you can sample from a multivariate normal distribution by using the RANDNORMAL module. Each of the resulting vectors is univariate normal and they are correlated with each other. You [...]

[...] Sampling from the multivariate normal distribution: This article describes how to obtain correlated multivariate normal vectors in SAS. [...]

[...] blogged several times about multivariate normality, including how to generate random values from a multivariate normal distribution. But given a set of multivariate data, how can you determine if it is likely to have come from a [...]

[...] about how to simulate various quantities. Random integers? Check! Random univariate samples? Check! Random multivariate samples? [...]

[...] statistical task varies, but one place where this problem occurs is in simulating multivariate normal data. I have previously written about why an estimated matrix of pairwise correlations is not always a [...]

[...] have previously shown how to use the RANDNORMAL function in SAS/IML to simulate multivariate normal data. Now suppose that you want to generate 10 samples, where each sample contains five observations [...]