Wygeneruj zmienną losową ze zdefiniowaną korelacją z istniejącą zmienną (zmiennymi)

Dla badań symulacyjnych mam do generowania zmiennych losowych, które wykazują prefined (populacji) korelację do istniejącej zmiennej .$Y$

I spojrzał w Ropakowaniach copula, a CDVinektóre mogą powodować przypadkowe wielowymiarowych rozkładów danej struktury zależności. Nie można jednak naprawić jednej z powstałych zmiennych do istniejącej zmiennej.

Wszelkie pomysły i linki do istniejących funkcji są mile widziane!

Wniosek: pojawiły się dwie ważne odpowiedzi, z różnymi rozwiązaniami:

1. R Scenariusz według caracal, co wylicza zmienną losową z dokładnym (próbki) korelacji ustalonej zmiennej
2. R Funkcja znalazłem się, co wylicza zmienną losową o określonej populacji korelacji do predefiniowanej zmiennej

[Dodanie @ttnphns: Zezwoliłem na rozszerzenie tytułu pytania z pojedynczej stałej zmiennej na dowolną liczbę stałych zmiennych; tj. jak wygenerować zmienną mającą predefiniowane korelacje z pewnymi stałymi, istniejącymi zmiennymi]

Oto kolejny: dla wektorów ze średnią 0 ich korelacja jest równa cosinus ich kąta. Zatem jednym sposobem na znalezienie wektora o dokładnie pożądanej korelacji , odpowiadającej kątowi :r $x$$r$$\theta$

1. uzyskaj ustalony wektor i losowy wektorx $x_1$$x_2$
2. wyśrodkuj oba wektory (średnia 0), dając wektory , ˙ x$\dot{x}_{1}$$\dot{x}_{2}$
3. uczyń prostopadłą do (rzut na podprzestrzeń ortogonalną), dając ˙ x 1 ˙ x ⊥ $\dot{x}_{2}$$\dot{x}_{1}$$\dot{x}_{2}^{\perp}$
4. przeskaluj i do długości 1, dając i ˙ 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}$ to wektor, którego kąt do to , a którego korelacja z jest więc . Jest to również korelacja z ponieważ transformacje liniowe pozostawiają korelację bez zmian.$\bar{x}_{1}$$\theta$$\bar{x}_{1}$$r$$x_1$

Oto kod:

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

Na rzut prostopadły , to używany Q R -decomposition poprawić stabilność liczbową, gdyż wtedy po prostu P = Q, Q ' .$P$$QR$$P = Q Q'$

Opiszę najogólniejsze możliwe rozwiązanie. Rozwiązanie problemu w tej ogólności pozwala nam osiągnąć niezwykle kompaktową implementację oprogramowania: wystarczy tylko dwie krótkie linie Rkodu.

Wybierz wektor o tej samej długości co , zgodnie z dowolnym rozkładem. Niech być pozostałości z regresji metodą najmniejszych kwadratów z z : ten wyodrębnia elementu z . Przez ponowne dodanie odpowiedniego wielokrotność do możemy wytworzenia wektora posiadającego dowolną korelacji z . Rozwiązaniem jest dowolna dowolna stała addytywna i dodatnia stała mnożąca - którą możesz dowolnie wybraćY Y ⊥ X Y Y X Y Y ⊥ ρ $X$$Y$$Y^\perp$$X$$Y$$Y$$X$$Y$$Y^\perp$$\rho$$Y$

(„ nazwa ” oznacza wszelkie obliczenia proporcjonalne do odchylenia standardowego.)$\operatorname{SD}$

Oto działający Rkod. Jeśli nie podasz , kod pobierze swoje wartości ze standardowego rozkładu normalnego na wielu odmianach.$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)
}

W celu zilustrowania, że generowane losowo z elementów, a wytwarzane o różnych określonych korelacji z tym . Wszystkie zostały utworzone przy użyciu tego samego wektora początkowego . Oto ich wykresy rozrzutu. „Wykresy rugowe” u dołu każdego panelu pokazują wspólny wektor50 X Y ; ρ Y X = ( 1 , 2 , … , 50 ) $Y$$50$$X_{Y;\rho}$$Y$$X=(1,2,\ldots, 50)$$Y$

Istnieje niezwykłe podobieństwo między fabułami, czyż nie :-).

Jeśli chcesz eksperymentować, oto kod, który wygenerował te dane i rysunek. (Nie zawracałem sobie głowy skorzystaniem ze swobody, aby przesuwać i skalować wyniki, które są łatwymi operacjami).

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

BTW, ta metoda z łatwością uogólnia na więcej niż jedno : jeśli jest to matematycznie możliwe, znajdzie po określeniu korelacji z całością zestaw . Wystarczy użyć zwykłych najmniejszych kwadratów, aby wyjąć efekty wszystkich z i utworzyć odpowiednią liniową kombinację i reszt. (Pomaga to zrobić w kategoriach podwójnej podstawy dla , która jest uzyskiwana przez obliczenie pseudo-odwrotności. Poniższy kod używa SVD dla osiągnięcia tego.)X Y 1 , Y 2 , … , Y k ; ρ 1 , ρ 2 , … , ρ k Y i Y i X Y i Y $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$

Oto szkic algorytmu, w Rktórym są podane jako kolumny macierzy :$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)

Poniżej znajduje się pełniejsza implementacja dla tych, którzy chcieliby eksperymentować.

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

Oto inne podejście obliczeniowe (rozwiązanie zostało zaadaptowane z postu na forum Enrico Schumanna). Według Wolfganga (patrz komentarze) jest to obliczeniowo identyczne z rozwiązaniem zaproponowanym przez ttnphns.

W przeciwieństwie do rozwiązania karakala nie wytwarza próbki o dokładnej korelacji , ale dwa wektory, których korelacja populacji jest równa .$\rho$$\rho$

Poniższa funkcja może obliczyć dwuwymiarowy rozkład próbek pobranych z populacji o danym . Oblicza dwie zmienne losowe lub pobiera jedną istniejącą zmienną (przekazaną jako parametr ) i tworzy drugą zmienną o pożądanej korelacji:$\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)
}

Funkcja może również wykorzystywać niestandardowe rozkłady brzeżne poprzez dostosowanie parametru mar.fun. Należy jednak pamiętać, że ustalenie jednej zmiennej tylko wydaje się działać z zmiennej o rozkładzie normalnym x! (co może odnosić się do komentarza Makra).

Należy również zauważyć, że „mały współczynnik korygujący” z pierwotnego postu został usunięty, ponieważ wydaje się, że przesądza powstałe korelacje, przynajmniej w przypadku rozkładów Gaussa i korelacji Pearsona (patrz także komentarze).

$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$

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

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

Zaktualizuj 11 listopada 2017 r. Dzisiaj spotkałem ten stary wątek i postanowiłem rozszerzyć moją odpowiedź, pokazując algorytm iteracyjnego dopasowania, o którym mówiłem na początku.

$Y$ $X$

Disclamer: To iteracyjne rozwiązanie, które znalazłem gorsze od doskonałego, oparte na znalezieniu podwójnej podstawy i zaproponowane przez @whuber w tym wątku dzisiaj. @ rozwiązanie Whubera nie jest iteracyjne i, co ważniejsze, wydaje mi się, że wpływa na wartości wejściowej zmiennej „świnia” nieco mniej niż algorytm „mój” (byłoby to atutem, gdyby zadaniem było „poprawić” istniejąca zmienna i nie generować losowych zmiennych od zera). Nadal publikuję moje z ciekawości i dlatego, że to działa (patrz także przypis).

$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}$$

   (mianownik nie zmienia się w iteracjach, oblicz go wcześniej)

   $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}$

   $m$$r$$SS_S$$n$

10. $C$$E$$r$$Y$$Y[\text{corrected}]=\hat{Y}+E$

11. $Y$

12. $Y$$r$

$Y$$r$$Y$

$^1$$Y$$X$

Miałem ochotę trochę programować, więc wziąłem usuniętą odpowiedź @ Adama i postanowiłem napisać fajną implementację w języku R. Skupiam się na używaniu stylu zorientowanego funkcjonalnie (tj. Pętli stylu lapply). Ogólna idea polega na pobraniu dwóch wektorów, losowym permutacji jednego z wektorów, dopóki nie zostanie osiągnięta pewna korelacja między nimi. To podejście jest bardzo brutalne, ale łatwe do wdrożenia.

Najpierw tworzymy funkcję, która losowo permutuje wektor wejściowy:

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

... i utwórz przykładowe dane

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

... napisz funkcję, która permutuje wektor wejściowy i koreluje go z wektorem referencyjnym:

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

... i iteruj tysiąc razy:

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

Zauważ, że reguły określania zakresu R zapewniają vec1i vec2znajdują się w środowisku globalnym, poza anonimową funkcją używaną powyżej. Zatem permutacje są względne w stosunku do oryginalnych zestawów danych testowych, które wygenerowaliśmy.

Następnie znajdujemy maksymalną korelację:

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

... lub znajdź wartość najbliższą korelacji 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

Aby uzyskać wyższą korelację, musisz zwiększyć liczbę iteracji.

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

Rozwiązanie:

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

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

Wyjście testowe:

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

Wygeneruj zmienne normalne z podaną macierzą kowariancji 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;

Wygeneruj zmienne normalne z podaną macierzą kowariancji LUDNOŚCI

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

Po prostu utwórz losowy wektor i sortuj, aż uzyskasz pożądane r.