Suppose I have generated a collection of correlated sequences of samples $(S_i)_{i=1}^{n}$ from random variables $\mathbf{\underline{x}} = x_i$.

Let's fix a sequence of reals $(\sigma_i)_{i=0}^{n}$. Suppose that I want to generate a sample $S_{i+1}$ such that $Correlation(S_{k},S_{i+1}) = \sigma_k$ for any k.

So basically I want to append one result on to collection of correlated samples and dictate what its correlation to each of the pre-existing ones should be. What are good ways of doing this?


3 Answers 3


As Richard says, this is really hard to do in a general setting. But if we make extra assumptions about the distribution of the variables, it might be doable.

Assume for instance that your variables are following a multivariate normal distribution. This is interesting because

  • The distribution is characterized exclusively by the means, variances and covariances.
  • The marginals of the distribution are again multivariate normal.
  • There is a simple formula for the conditional distributions of variables that are multivariate normal.

Using the conditional distribution, you could easily generate a new sample $S_{n+1}$ assuming that the sample is from a random variable, part of a $n+1$-size multivariate normal set of random variables.

  • $\begingroup$ Ah I see what you mean. I think I can make do with MVN for a good part of it so we can use the conditional distribution. Thanks! $\endgroup$
    – Djoker
    Jan 17, 2018 at 9:23

I think it is hard to add a random variable $X$ with a predefined correaltion to a whole sample $(X_1, \ldots, X_n)$ because this would mean that you have to define relations to each of the $n$ existing rvs which could be infeasible.

A partial answer is the following: For a random variables $X$ and $Y$ uncorrelated with variance $1$ you can do the following. Choose a correlation $\rho$ and define $$ Z := \sqrt{1-\rho^2} X + \rho Y. $$ Note that $var(Z)=1$.

Then $cor(Z,Y) = cov(Z,Y) = \rho$ (cov = cor as variance is 1). To see this calculate $$ cov(Z,Y) = cov(\sqrt{1-\rho^2} X + \rho Y,Y) = cov(\sqrt{1-\rho^2} X,Y) + cov(\rho Y,Y) = \rho, $$
because $cov(X,Y) = 0$ (they are uncorrelated) and $cov(\rho Y,Y) =\rho\cdot cov(Y,Y) = \rho $.

Thus you can sample $X$ uncorrelated and define $Z$ above. Then you get a correlation to one of the others.


As other answers have pointed out, you cannot in general impose any arbitrary correlation structure on your samples. But you can try to rearrange your new sample in such a way that you get as close as possible to the desired correlation.

The idea is this: given your existing samples, generate a new one with the desired marginal distribution. Suppose you reordered this new sample, e.g. reshuffled it. The marginal distribution would be unchanged, but the correlations with the existing samples change. You could now try to search for an reordering that gets you correlations close to your desired ones.

Here is a sketch how you could do this in R. Suppose you have three samples, drawn from a standard Gaussian distribution with the following correlation matrix:

     [,1] [,2] [,3]
[1,]  1.0  0.1  0.5
[2,]  0.1  1.0  0.7
[3,]  0.5  0.7  1.0

Let's create the dataset.

C <- c(1.0, 0.1, 0.5,
       0.1, 1.0, 0.7, 
       0.5, 0.7, 1.0)
dim(C) <- c(3,3)

n <- 100  ## number of observations
A <- array(rnorm(n*3), dim = c(n, 3))
A <- A %*% chol(C)

##         [,1]    [,2]  [,3]
## [1,] 1.00000 0.08617 0.592
## [2,] 0.08617 1.00000 0.635
## [3,] 0.59197 0.63503 1.000

Create the new sample and define target correlations.

new <- runif(n)
target <- c(0.2,0.3,0.6)

The new sample comes from a uniform distribution, but the method does not depend on the distribution.

Now we need to find an ordering of the values in new that gets us close to the desired correlation. I do this with a local search, as implemented in the R package NMOF (of which I am the maintainer).

For this, I first collect all input to the search in a list Data.

Data <- list()
Data$new <- new
Data$target <- target
Data$A <- A

I define an objective function: it computes the average absolute difference between realised and target correlation for a given reordering x:

mean_diff <- function(x, Data)
    sum(abs(cor(Data$new[x], Data$A) - Data$target))

The final ingredient is a neighbourhood function. It takes a solution and slightly changes it by exchanging two elements.

neighbour <- function(x, ...) {
    ans <- x
    i <- sample(length(x), 2)
    ans[i] <- rev(x[i])

neighbour(1:5)  ## example
## 3 2 1 4 5

It remains to run the search:

sol <- LSopt(mean_diff,
             list(neighbour = neighbour, nS = 1000, x0 = 1:n),
             Data = Data)

cor(new[sol$xbest], A)
##        [,1]   [,2]   [,3]
## [1,] 0.2004 0.2996 0.5986

Which is pretty close to the defined target.


Your Answer

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

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