Генерация случайной величины с определенной корреляцией с существующей переменной

Для исследования моделирования я должен генерировать случайные переменные , которые показывают prefined (населения) корреляцию с существующей переменной .$Y$

Я посмотрел на Rпакеты copulaи CDVineкоторые могут производить случайные многомерные распределения с заданной структурой зависимостей. Однако невозможно зафиксировать одну из результирующих переменных в существующей переменной.

Любые идеи и ссылки на существующие функции приветствуются!

Вывод: пришли два правильных ответа с разными решениями:

1. R Сценарий по каракал, который вычисляет случайную переменную с точным (образец) корреляции с предопределенной переменной
2. R Функция я очутилась, который вычисляет случайную величину с определенной популяцией корреляцией с предопределенным переменным

[@ttnphns 'добавление: я взял на себя смелость расширить заголовок вопроса с одного случая с фиксированной переменной на произвольное количество фиксированных переменных; то есть, как генерировать переменную, имеющую предопределенную корреляцию (и) с некоторой фиксированной, существующей переменной (ами)]

Вот еще один: для векторов со средним 0 их корреляция равна косинусу их угла. Итак, один из способов найти вектор с точно желаемой корреляцией , соответствующей углу :r $x$$r$$\theta$

1. получить фиксированный вектор и случайный векторх $x_1$$x_2$
2. центрировать оба вектора (в среднем 0), задавая векторы , ˙ х$\dot{x}_{1}$$\dot{x}_{2}$
3. сделать ортогональным (проекция на ортогональное подпространство), получив ˙ x 1 ˙ x ⊥ $\dot{x}_{2}$$\dot{x}_{1}$$\dot{x}_{2}^{\perp}$
4. масштабировать и до длины 1, получая и ˙ x ⊥ 2 ˉ x 1 ˉ x ⊥ $\dot{x}_{1}$$\dot{x}_{2}^{\perp}$$\bar{x}_{1}$$\bar{x}_{2}^{\perp}$
5. ˉ x 1θ ˉ x 1rx$\bar{x}_{2}^{\perp} + (1/\tan(\theta)) \cdot \bar{x}_{1}$ - это вектор, угол которого к есть , и чья корреляция с таким образом, равна . Это также корреляция с поскольку линейные преобразования оставляют корреляцию без изменений.$\bar{x}_{1}$$\theta$$\bar{x}_{1}$$r$$x_1$

Вот код:

n     <- 20                    # length of vector
rho   <- 0.6                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- rnorm(n, 1, 1)        # fixed given data
x2    <- rnorm(n, 2, 0.5)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho

Для ортогональной проекции я использовал разложение для улучшения числовой устойчивости, поскольку тогда просто .Q R P = Q Q $P$$QR$$P = Q Q'$

Я опишу наиболее общее возможное решение. Решение проблемы в этой общности позволяет нам достичь удивительно компактной программной реализации: достаточно двух коротких строк Rкода.

Выберите вектор той же длины, что и , в соответствии с любым распределением, которое вам нравится. Пусть быть остатки регрессии наименьших квадратов против : это извлекает компонент из . Добавляя назад подходящее кратное в , мы можем производить вектор , имеющий любую требуемую корреляционную с . До произвольной аддитивной константы и положительной мультипликативной константы - которую вы можете выбрать любым способом - решение -Y Y ⊥ X Y Y X Y Y ⊥ ρ $X$$Y$$Y^\perp$$X$$Y$$Y$$X$$Y$$Y^\perp$$\rho$$Y$

(« » означает любой расчет, пропорциональный стандартному отклонению.)$\operatorname{SD}$

Вот рабочий Rкод. Если вы не предоставите , код будет извлекать свои значения из многомерного стандартного нормального распределения.$X$

complement <- function(y, rho, x) {
  if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if `x` is not given
  y.perp <- residuals(lm(x ~ y))
  rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2)
}

Чтобы проиллюстрировать это , я произвел случайный с компонентами и производится , имеющую различные заданные корреляции с этим . Все они были созданы с одинаковым начальным вектором . Вот их диаграммы рассеяния. «Коврики» внизу каждой панели показывают общий вектор$Y$$50$$X_{Y;\rho}$$Y$$X=(1,2,\ldots, 50)$$Y$

Среди сюжетов есть замечательное сходство, не так ли :-).

Если вы хотите поэкспериментировать, вот код, который создал эти данные, и рисунок. (Я не удосужился использовать свободу для сдвига и масштабирования результатов, что является простой операцией.)

y <- rnorm(50, sd=10)
x <- 1:50 # Optional
rho <- seq(0, 1, length.out=6) * rep(c(-1,1), 3)
X <- data.frame(z=as.vector(sapply(rho, function(rho) complement(y, rho, x))),
                rho=ordered(rep(signif(rho, 2), each=length(y))),
                y=rep(y, length(rho)))

library(ggplot2)
ggplot(X, aes(y,z, group=rho)) + 
  geom_smooth(method="lm", color="Black") + 
  geom_rug(sides="b") + 
  geom_point(aes(fill=rho), alpha=1/2, shape=21) +
  facet_wrap(~ rho, scales="free")

Кстати, этот метод легко обобщает более чем на один : если это математически возможно, он найдет с указанными корреляциями с целым набор . Просто используйте обычные наименьшие квадраты, чтобы убрать эффекты всех из и сформировать подходящую линейную комбинацию и остатков. (Это помогает сделать это с точки зрения двойного базиса для , который получается путем вычисления псевдообратного кода. Следующий код использует SVD для для достижения этой цели.)$Y$$X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$$Y_i$$Y_i$$X$$Y_i$$Y$$Y$

Вот эскиз алгоритма, в Rкотором представлены в виде столбцов матрицы :$Y_i$y

y <- scale(y)             # Makes computations simpler
e <- residuals(lm(x ~ y)) # Take out the columns of matrix `y`
y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
return(y.dual %*% rho + sqrt(sigma2)*e)

Ниже приведена более полная реализация для тех, кто хотел бы поэкспериментировать.

complement <- function(y, rho, x) {
  #
  # Process the arguments.
  #
  if(!is.matrix(y)) y <- matrix(y, ncol=1)
  if (missing(x)) x <- rnorm(n)
  d <- ncol(y)
  n <- nrow(y)
  y <- scale(y) # Makes computations simpler
  #
  # Remove the effects of `y` on `x`.
  #
  e <- residuals(lm(x ~ y))
  #
  # Calculate the coefficient `sigma` of `e` so that the correlation of
  # `y` with the linear combination y.dual %*% rho + sigma*e is the desired
  # vector.
  #
  y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
  sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
  #
  # Return this linear combination.
  #
  if (sigma2 >= 0) {
    sigma <- sqrt(sigma2) 
    z <- y.dual %*% rho + sigma*e
  } else {
    warning("Correlations are impossible.")
    z <- rep(0, n)
  }
  return(z)
}
#
# Set up the problem.
#
d <- 3           # Number of given variables
n <- 50          # Dimension of all vectors
x <- 1:n         # Optionally: specify `x` or draw from any distribution
y <- matrix(rnorm(d*n), ncol=d) # Create `d` original variables in any way
rho <- c(0.5, -0.5, 0)          # Specify the correlations
#
# Verify the results.
#
z <- complement(y, rho, x)
cbind('Actual correlations' = cor(cbind(z, y))[1,-1],
      'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))

Вот еще один вычислительный подход (решение адаптировано из сообщения на форуме Энрико Шумана). Согласно Вольфгангу (см. Комментарии), это вычислительно идентично решению, предложенному ttnphns.

В отличие от решения Каракала, он не дает выборку с точной корреляцией , но два вектора, корреляция населения которых равна .$\rho$$\rho$

Следующая функция может вычислить двумерное распределение выборки, взятой из совокупности с заданным значением . Он либо вычисляет две случайные переменные, либо берет одну существующую переменную (переданную в качестве параметра ) и создает вторую переменную с желаемой корреляцией:$\rho$x

# returns a data frame of two variables which correlate with a population correlation of rho
# If desired, one of both variables can be fixed to an existing variable by specifying x
getBiCop <- function(n, rho, mar.fun=rnorm, x = NULL, ...) {
     if (!is.null(x)) {X1 <- x} else {X1 <- mar.fun(n, ...)}
     if (!is.null(x) & length(x) != n) warning("Variable x does not have the same length as n!")

     C <- matrix(rho, nrow = 2, ncol = 2)
     diag(C) <- 1

     C <- chol(C)

     X2 <- mar.fun(n)
     X <- cbind(X1,X2)

     # induce correlation (does not change X1)
     df <- X %*% C

     ## if desired: check results
     #all.equal(X1,X[,1])
     #cor(X)

     return(df)
}

Функция также может использовать ненормальные предельные распределения, регулируя параметр mar.fun. Обратите внимание, однако, что исправление одной переменной работает только с нормально распределенной переменной x! (что может относиться к комментарию Макроса).

Также обратите внимание, что «небольшой поправочный коэффициент» из исходного поста был удален, так как он, кажется, смещает результирующие корреляции, по крайней мере, в случае гауссовых распределений и корреляций Пирсона (также см. Комментарии).

Пусть будет вашей фиксированной переменной, и вы хотите сгенерировать переменную которая коррелирует с на величину . Если стандартизирован, то (потому что - это бета-коэффициент в простой регрессии) , где - случайная величина из нормального распределения, имеющая среднее значение и . Наблюдаемая корреляция между данными и будет приблизительно равна ; и можно рассматривать как случайные выборки из двумерной нормальной популяции (если$X$$Y$$X$$r$$X$$r$$Y= rX+E$$E$$0$$\text{sd}=\sqrt{1-r^2}$$X$$Y$$r$$X$$Y$$X$ от нормального) с .$\rho=r$

E $r$$E$$X$$E$$X$$Y$$X_1, X_2, X_3,...$

$X$$r$$Y$$Y$$r$$Y$

Обновление 11 ноября 2017 г. Сегодня я наткнулся на эту старую ветку и решил расширить свой ответ, показав алгоритм итеративной подгонки, о котором я говорил вначале.

$Y$ $X$

Отказ от ответственности: Это итеративное решение, которое я нашел, уступает превосходному, основанному на поиске двойного базиса и предложенному @whuber в этой теме сегодня. Решение @ whuber не является итеративным, и, что для меня более важно, оно, по-видимому, влияет на значения входной переменной «pig» несколько меньше, чем алгоритм «my» (тогда было бы полезно, если бы задача «исправить») существующая переменная, а не генерировать случайные изменения с нуля). Тем не менее, я публикую свою для любопытства и потому, что это работает (см. Также сноску).

$X_1, X_2,...,X_m$$Y$$Y$$r_1, r_2,...,r_m$$X$

$Y$$X$$Y$$Y$

1. $r$$\text{df}=n-1$$S_j=r_j \text{df}$$j$$X$

2. $\text{df}$$Y$$X$$\text{df}$

3. $Y$$X$$r$$\bf b=(X'X)^{-1} S$

4. $Y$$\hat{Y}=\bf Xb$

5. $E=Y-\hat{Y}$

6. $SS_S=\text{df}-SS_{\hat {Y}}$

7. $E$$X_j$$C_j= \sum_{i=1}^n E_i X_{ij}$

8. $E$$C$$0$$i$

   $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m C_j X_{ij}} {n\sum_{j=1}^m X_{ij}^2}$$

   (знаменатель не меняется на итерациях, рассчитайте его заранее)

   Или, альтернативно, более эффективная формула дополнительно гарантирует, что среднее значение станет . Сначала сделайте центр на каждой итерации перед вычислением s на шаге 7, затем на этом шаге 8 исправьте как:$E$$0$ $E$$C$

   $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m \frac{C_j X_{ij}^3}{\sum_{i=1}^n X_{ij}^2}} {\sum_{j=1}^m X_{ij}^2}$$

   (опять же знаменатели известны заранее)$^1$

9. Приведите к его целевому значению:$SS_E$$E_i[\text{corrected}]=E_i \sqrt{SS_S/SS_E}$

   Перейдите к шагу 7. (Скажем, сделайте, например, 10-20 итераций; чем больше тем больше итераций может понадобиться. Если целевой был реалистичным, положителен, а если размер выборки не слишком мал, итерации всегда прямой к конвергенции. Конец итерации.)$m$$r$$SS_S$$n$

10. Готов: все теперь почти равны нулю, что означает, что остатки обучены восстанавливать целевые значения . Вычислить фитинг : .Е Г У У [ исправлено ] = У + $C$$E$$r$$Y$$Y[\text{corrected}]=\hat{Y}+E$

11. Полученный практически стандартизирован. В качестве последнего штриха, вы можете точно стандартизировать его, опять же, как вы это делали на шаге 2.$Y$

12. Вы можете указать с любым отклонением и означать, что вам нравится. Собственно, среди четырех статистик - мин , макс , скуп , ст. DEV . - вы можете выбрать любые два значения и линейно преобразовать переменную так, чтобы она представляла их без изменения s (корреляций), которые вы достигли (это все называется линейным масштабированием).$Y$$r$

Чтобы предупредить снова, что было сказано выше. С этим потягиванием точно к , выход не должен быть нормально распределен.r $Y$$r$$Y$

Y $^1$ Формула коррекции может быть еще более сложной, например, чтобы обеспечить большую гомоскедастичность (в терминах сумм квадратов) также для каждого , одновременно с получением корреляций, - я реализовал код для этого слишком. (Я не знаю, если такая «двойная» задача разрешима с помощью более аккуратного - не итеративного - подхода, такого как Whuber .)$Y$$X$

Мне захотелось немного программировать, поэтому я взял удаленный ответ @ Adam и решил написать хорошую реализацию на R. Я сосредотачиваюсь на использовании функционально ориентированного стиля (т. Е. Цикла цикла lapply). Общая идея состоит в том, чтобы взять два вектора, случайным образом переставляя один из векторов, пока между ними не будет достигнута определенная корреляция. Этот подход очень грубый, но простой в реализации.

Сначала мы создаем функцию, которая случайным образом переставляет входной вектор:

randomly_permute = function(vec) vec[sample.int(length(vec))]
randomly_permute(1:100)
  [1]  71  34   8  98   3  86  28  37   5  47  88  35  43 100  68  58  67  82
 [19]  13   9  61  10  94  29  81  63  14  48  76   6  78  91  74  69  18  12
 [37]   1  97  49  66  44  40  65  59  31  54  90  36  41  93  24  11  77  85
 [55]  32  79  84  15  89  45  53  22  17  16  92  55  83  42  96  72  21  95
 [73]  33  20  87  60  38   7   4  52  27   2  80  99  26  70  50  75  57  19
 [91]  73  62  23  25  64  51  30  46  56  39

... и создать пример данных

vec1 = runif(100)
vec2 = runif(100)

... написать функцию, которая переставляет входной вектор и соотносит его с опорным вектором:

permute_and_correlate = function(vec, reference_vec) {
    perm_vec = randomly_permute(vec)
    cor_value = cor(perm_vec, reference_vec)
    return(list(vec = perm_vec, cor = cor_value))
  }
permute_and_correlate(vec2, vec1)
$vec
  [1] 0.79072381 0.23440845 0.35554970 0.95114398 0.77785348 0.74418811
  [7] 0.47871491 0.55981826 0.08801319 0.35698405 0.52140366 0.73996913
 [13] 0.67369873 0.85240338 0.57461506 0.14830718 0.40796732 0.67532970
 [19] 0.71901990 0.52031017 0.41357545 0.91780357 0.82437619 0.89799621
 [25] 0.07077250 0.12056045 0.46456652 0.21050067 0.30868672 0.55623242
 [31] 0.84776853 0.57217746 0.08626022 0.71740151 0.87959539 0.82931652
 [37] 0.93903143 0.74439384 0.25931398 0.99006038 0.08939812 0.69356590
 [43] 0.29254936 0.02674156 0.77182339 0.30047034 0.91790830 0.45862163
 [49] 0.27077191 0.74445997 0.34622648 0.58727094 0.92285322 0.83244284
 [55] 0.61397396 0.40616274 0.32203732 0.84003379 0.81109473 0.50573325
 [61] 0.86719899 0.45393971 0.19701975 0.63877904 0.11796154 0.26986325
 [67] 0.01581969 0.52571331 0.27087693 0.33821824 0.52590383 0.11261002
 [73] 0.89840404 0.82685046 0.83349287 0.46724807 0.15345334 0.60854785
 [79] 0.78854984 0.95770015 0.89193212 0.18885955 0.34303707 0.87332019
 [85] 0.08890968 0.22376395 0.02641979 0.43377516 0.58667068 0.22736077
 [91] 0.75948043 0.49734797 0.25235660 0.40125309 0.72147500 0.92423638
 [97] 0.27980561 0.71627101 0.07729027 0.05244047

$cor
[1] 0.1037542

... и повторять тысячу раз:

n_iterations = lapply(1:1000, function(x) permute_and_correlate(vec2, vec1))

Обратите внимание , что правила области АиР гарантировать , что vec1и vec2находятся в глобальной среде, вне анонимной функции , используемой выше. Итак, все перестановки относятся к исходным наборам тестовых данных, которые мы сгенерировали.

Далее находим максимальную корреляцию:

cor_values = sapply(n_iterations, '[[', 'cor')
n_iterations[[which.max(cor_values)]]
$vec
  [1] 0.89799621 0.67532970 0.46456652 0.75948043 0.30868672 0.83244284
  [7] 0.86719899 0.55623242 0.63877904 0.73996913 0.71901990 0.85240338
 [13] 0.81109473 0.52571331 0.82931652 0.60854785 0.19701975 0.26986325
 [19] 0.58667068 0.52140366 0.40796732 0.22736077 0.74445997 0.40125309
 [25] 0.89193212 0.52031017 0.92285322 0.91790830 0.91780357 0.49734797
 [31] 0.07729027 0.11796154 0.69356590 0.95770015 0.74418811 0.43377516
 [37] 0.55981826 0.93903143 0.30047034 0.84776853 0.32203732 0.25235660
 [43] 0.79072381 0.58727094 0.99006038 0.01581969 0.41357545 0.52590383
 [49] 0.27980561 0.50573325 0.92423638 0.11261002 0.89840404 0.15345334
 [55] 0.61397396 0.27077191 0.12056045 0.45862163 0.18885955 0.77785348
 [61] 0.23440845 0.05244047 0.25931398 0.57217746 0.35554970 0.34622648
 [67] 0.21050067 0.08890968 0.84003379 0.95114398 0.83349287 0.82437619
 [73] 0.46724807 0.02641979 0.71740151 0.74439384 0.14830718 0.82685046
 [79] 0.33821824 0.71627101 0.77182339 0.72147500 0.08801319 0.08626022
 [85] 0.87332019 0.34303707 0.45393971 0.47871491 0.29254936 0.08939812
 [91] 0.35698405 0.67369873 0.27087693 0.78854984 0.87959539 0.22376395
 [97] 0.02674156 0.07077250 0.57461506 0.40616274

$cor
[1] 0.3166681

... или найти наиболее близкое значение к соотношению 0,2:

n_iterations[[which.min(abs(cor_values - 0.2))]]
$vec
  [1] 0.02641979 0.49734797 0.32203732 0.95770015 0.82931652 0.52571331
  [7] 0.25931398 0.30047034 0.55981826 0.08801319 0.29254936 0.23440845
 [13] 0.12056045 0.89799621 0.57461506 0.99006038 0.27077191 0.08626022
 [19] 0.14830718 0.45393971 0.22376395 0.89840404 0.08890968 0.15345334
 [25] 0.87332019 0.92285322 0.50573325 0.40796732 0.91780357 0.57217746
 [31] 0.52590383 0.84003379 0.52031017 0.67532970 0.83244284 0.95114398
 [37] 0.81109473 0.35554970 0.92423638 0.83349287 0.34622648 0.18885955
 [43] 0.61397396 0.89193212 0.74445997 0.46724807 0.72147500 0.33821824
 [49] 0.71740151 0.75948043 0.52140366 0.69356590 0.41357545 0.21050067
 [55] 0.87959539 0.11796154 0.73996913 0.30868672 0.47871491 0.63877904
 [61] 0.22736077 0.40125309 0.02674156 0.26986325 0.43377516 0.07077250
 [67] 0.79072381 0.08939812 0.86719899 0.55623242 0.60854785 0.71627101
 [73] 0.40616274 0.35698405 0.67369873 0.82437619 0.27980561 0.77182339
 [79] 0.19701975 0.82685046 0.74418811 0.58667068 0.93903143 0.74439384
 [85] 0.46456652 0.85240338 0.34303707 0.45862163 0.91790830 0.84776853
 [91] 0.78854984 0.05244047 0.58727094 0.77785348 0.01581969 0.27087693
 [97] 0.07729027 0.71901990 0.25235660 0.11261002

$cor
[1] 0.2000199

Чтобы получить более высокую корреляцию, вам нужно увеличить количество итераций.

$Y_1$$Y_2,\dots,Y_n$$R$

Решение:

1. $CC^T=R$
2. $X_2,\dots,X_n$$Y_1$
3. $Y_1$
4. $Y=CX$$Y_i$$Y_1$

Код Python:

import numpy as np
import math
from scipy.linalg import toeplitz, cholesky
from statsmodels.stats.moment_helpers import cov2corr

# create the large correlation matrix R
p = 4
h = 2/p
v = np.linspace(1,-1+h,p)
R = cov2corr(toeplitz(v))

# create the first variable
T = 1000;
y = np.random.randn(T)

# generate p-1 correlated randoms
X = np.random.randn(T,p)
X[:,0] = y
C = cholesky(R)
Y = np.matmul(X,C)

# check that Y didn't change
print(np.max(np.abs(Y[:,0]-y)))

# check the correlation matrix
print(R)
print(np.corrcoef(np.transpose(Y)))

Тестовый вывод:

0.0
[[ 1.   0.5  0.  -0.5]
 [ 0.5  1.   0.5  0. ]
 [ 0.   0.5  1.   0.5]
 [-0.5  0.   0.5  1. ]]
[[ 1.          0.50261766  0.02553882 -0.46259665]
 [ 0.50261766  1.          0.51162821  0.05748082]
 [ 0.02553882  0.51162821  1.          0.51403266]
 [-0.46259665  0.05748082  0.51403266  1.        ]]

Генерация нормальных переменных с ковариационной матрицей SAMPLING, как указано

covsam <- function(nobs,covm, seed=1237) {; 
          library (expm);
          # nons=number of observations, covm = given covariance matrix ; 
          nvar <- ncol(covm); 
          tot <- nvar*nobs;
          dat <- matrix(rnorm(tot), ncol=nvar); 
          covmat <- cov(dat); 
          a2 <- sqrtm(solve(covmat)); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% a2 %*% m2 ; 
          rc <- cov(dat2);};
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covsam(10,cm)  ;
          res;

Генерация нормальных переменных с ковариационной матрицей численности населения, как указано

covpop <- function(nobs,covm, seed=1237) {; 
          library (expm); 
          # nons=number of observations, covm = given covariance matrix;
          nvar <- ncol(covm); 
          tot <- nvar*nobs;  
          dat <- matrix(rnorm(tot), ncol=nvar); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% m2;  
          rc <- cov(dat2); }; 
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covpop(10,cm); 
          res

Просто создайте случайный вектор и сортируйте, пока не получите желаемый r.