Gere uma variável aleatória com uma correlação definida para uma (s) variável (s) existente (s)

Para um estudo de simulação, eu tenho que gerar variáveis ​​aleatórias que mostram uma correlação pré-definida (população) com uma variável existente $Y$.

Examinei os Rpacotes copulae CDVineque podem produzir distribuições multivariadas aleatórias com uma determinada estrutura de dependência. No entanto, não é possível corrigir uma das variáveis ​​resultantes em uma variável existente.

Todas as idéias e links para as funções existentes são apreciadas!

Conclusão: Surgiram duas respostas válidas, com diferentes soluções:

1. Um R script de caracal, que calcula uma variável aleatória com uma correlação exata (de amostra) a uma variável predefinida
2. Uma R função que encontrei, que calcula uma variável aleatória com uma correlação populacional definida com uma variável predefinida

[Além do @ttnphns: tomei a liberdade de expandir o título da pergunta de uma única variável fixa para um número arbitrário de variáveis ​​fixas; ou seja, como gerar uma variável com corretação (ões) predefinida (s) com algumas variáveis ​​fixas existentes

Aqui está outro: para vetores com média 0, sua correlação é igual ao cosseno de seu ângulo. Então, uma maneira de encontrar um vetor com exatamente a correlação desejada r , correspondente a um ângulo θ :$x$$r$$\theta$

1. obter vetor fixo e um vetor aleatório x $x_1$$x_2$
2. centralize ambos os vetores (média 0), fornecendo vetores , ˙ x$\dot{x}_{1}$$\dot{x}_{2}$
3. faça ortogonal a ˙ x 1 (projeção no subespaço ortogonal), dando ˙ x ⊥ $\dot{x}_{2}$$\dot{x}_{1}$$\dot{x}_{2}^{\perp}$
4. escalar e ˙ x ⊥ 2 para o comprimento 1, dando ˉ x 1 e ˉ x ⊥ $\dot{x}_{1}$$\dot{x}_{2}^{\perp}$$\bar{x}_{1}$$\bar{x}_{2}^{\perp}$
5. é o vector cujo ângulo a ° x 1éθ, e cuja correlação com ˉ x 1 , portanto, ér. Essa também é a correlação parax1,pois as transformações lineares mantêm a correlação inalterada.$\bar{x}_{2}^{\perp} + (1/\tan(\theta)) \cdot \bar{x}_{1}$$\bar{x}_{1}$$\theta$$\bar{x}_{1}$$r$$x_1$

Aqui está o código:

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

Para a projecção ortogonal , I utilizado o Q R -decomposition para melhorar a estabilidade numérico, uma vez que, em seguida, simplesmente P = Q Q ' .$P$$QR$$P = Q Q'$

Vou descrever a solução mais geral possível. Resolver o problema nessa generalidade nos permite obter uma implementação de software notavelmente compacta: apenas duas pequenas linhas de Rcódigo são suficientes.

Escolha um vetor , do mesmo comprimento que Y , de acordo com a distribuição que desejar. Vamos Y ⊥ ser os resíduos da regressão de mínimos quadrados de X contra Y : esta extrai a Y componente de X . Por via da adição de um múltiplo apropriado de Y para Y ⊥ , que pode produzir um vector com qualquer correlação desejada ρ com Y . Até uma constante aditiva arbitrária e uma constante multiplicativa positiva - que você pode escolher de qualquer forma - a solução é$X$$Y$$Y^\perp$$X$$Y$$Y$$X$$Y$$Y^\perp$$\rho$$Y$

(" " significa qualquer cálculo proporcional a um desvio padrão.)$\operatorname{SD}$

Aqui está o Rcódigo de trabalho . Se você não fornecer , o código extrairá seus valores da distribuição normal padrão multivariada.$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)
}

Para ilustrar, gerei um aleatório com 50 componentes e produzi X Y ; ρ tendo várias correlações indicados com esta Y . Todos foram criados com o mesmo vetor inicial X = ( 1 , 2 , … , 50 ) . Aqui estão os gráficos de dispersão. Os "rugplots" na parte inferior de cada painel mostram o vetor Y comum .$Y$$50$$X_{Y;\rho}$$Y$$X=(1,2,\ldots, 50)$$Y$

Há uma notável semelhança entre as parcelas, não é :-).

Se você deseja experimentar, aqui está o código que produziu esses dados e a figura. (Não me preocupei em usar a liberdade para alterar e escalar os resultados, que são operações fáceis.)

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

Aliás, esse método generaliza prontamente para mais de um : se for matematicamente possível, ele encontrará um X Y 1 , Y 2 , ... , Y k ; ρ 1 , ρ 2 , … , ρ k tendo correlações especificadas com um conjunto inteiro de Y i . Apenas use mínimos quadrados comuns para eliminar os efeitos de todo o Y i de X e formar uma combinação linear adequada do Y $Y$$X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$$Y_i$$Y_i$$X$$Y_i$e os resíduos. (Ajuda a fazer isso em termos de uma base dupla para , que é obtida computando uma pseudo-inversa. O código a seguir usa o SVD de Y para fazer isso.)$Y$$Y$

Aqui está um esboço do algoritmo em R, onde o são dadas como colunas de uma matriz :$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)

A seguir, é apresentada uma implementação mais completa para aqueles que desejam experimentar.

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

Aqui está outra abordagem computacional (a solução é adaptada de uma postagem no fórum de Enrico Schumann). Segundo Wolfgang (ver comentários), isso é computacionalmente idêntico à solução proposta por ttnphns.

Em contraste com a solução de caracal, não produz uma amostra com a correlação exata de , mas dois vetores cuja correlação populacional é igual a ρ .$\rho$$\rho$

A função a seguir pode calcular uma distribuição de amostra bivariada obtida de uma população com um dado . Ele calcula duas variáveis ​​aleatórias ou pega uma variável existente (passada como parâmetro ) e cria uma segunda variável com a correlação desejada:$\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)
}

A função também pode usar distribuições marginais não normais ajustando o parâmetro mar.fun. Note, no entanto, que a fixação de uma variável única parece funcionar com uma variável normalmente distribuída x! (que pode estar relacionado ao comentário da macro).

Observe também que o "pequeno fator de correção" do post original foi removido, pois parece influenciar as correlações resultantes, pelo menos no caso de distribuições gaussianas e correlações de Pearson (também ver comentários).

Seja sua variável fixa e você deseja gerar a variável Y que se correlaciona com X pela quantidade r . Se X é padronizado, então (porque r é coeficiente beta em regressão simples) Y = r X + E , onde E é variável aleatória da distribuição normal com média 0 e sd = $X$$Y$$X$$r$$X$$r$$Y= rX+E$$E$$0$ . A correlação observada entre osdadosXeYserá aproximadamenter; XeYpodem ser vistos como amostras aleatórias dapopulaçãonormal bivariada(seXé normal) comρ=r.$\text{sd}=\sqrt{1-r^2}$$X$$Y$$r$$X$$Y$$X$$\rho=r$

Agora, se você deseja obter a correlação na sua amostra bivariada exatamente , é preciso prever que E tem de zero correlação com X . Esse aperto em zero pode ser alcançado com a modificação E iterativamente. Bem, com apenas duas variáveis, uma dada ( X ) e um para gerar ( Y ), o número suficiente de iterações é, na verdade, um, mas com múltiplas variáveis de dados ( X 1 , X 2 , X 3 , . . . ) Iterações ser necessário.$r$$E$$X$$E$$X$$Y$$X_1, X_2, X_3,...$

Deve-se notar que, se é normal, no primeiro procedimento (" r aproximado ") Y também será normal; no entanto, no ajuste iterativo de Y para o " r exato " Y provavelmente perderá a normalidade porque o ajuste explora os valores de caso seletivamente.$X$$r$$Y$$Y$$r$$Y$

Atualize 11 de novembro de 2017. Encontrei esse tópico antigo hoje e decidi expandir minha resposta, mostrando o algoritmo do ajuste iterativo sobre o qual eu estava falando inicialmente.

Aqui está uma solução iterativa de como treinar uma variável simulada ou preexistente aleatoriamente para correlacionar ou covaria exatamente como desejamos (ou muito próximo disso - dependendo do número de iterações) com um conjunto de variáveis X s (elas não podem ser modificadas).$Y$ $X$

Disclamer: Esta solução iterativa que eu achei inferior à excelente, com base em encontrar a base dupla e proposta por @whuber neste tópico hoje. A solução @ whuber não é iterativa e, mais importante para mim, parece estar afetando os valores da variável "pig" de entrada um pouco menos do que o algoritmo "my" (seria um recurso, se a tarefa fosse "corrigir" variável existente e não gerar variável aleatória a partir do zero). Ainda assim, estou publicando o meu por curiosidade e porque funciona (veja também Nota de rodapé).

Então, temos dado (fixos) variáveis , e varible Y que é ou apenas gerado aleatoriamente "porco" de valores ou é uma variável de dados existente quais os valores que precisam de "correcta" - para trazer Y exactamente a correlações (ou pode ser covariâncias) r 1 , r 2 , . . . , r m com os X s. Todos os dados devem ser contínuos; em outras palavras, deve haver uma grande quantidade de valores únicos.$X_1, X_2,...,X_m$$Y$$Y$$r_1, r_2,...,r_m$$X$

A ideia: realizar ajuste iterativo de resíduos. Conhecendo as correlações / covariâncias desejadas (alvo), podemos calcular os valores previstos para o usando os Xs como múltiplos preditores lineares. Depois de obter os resíduos iniciais (do Y atual e da previsão ideal), treine-os iterativamente para não se correlacionar com os preditores. No final, recupere Y com os resíduos. (O procedimento foi minha própria invenção experimental da roda, muitos anos atrás, quando eu não conhecia a teoria; eu a codifiquei no SPSS.)$Y$$X$$Y$$Y$

1. Converta o destino s em somas de produtos cruzados multiplicando-os por df = n - 1 : S j = r j df . ( j é um índice de variável X ).$r$$\text{df}=n-1$$S_j=r_j \text{df}$$j$$X$

2. Padronize Z todas as variáveis ​​(centralize cada uma delas e divida pelo desvio padrão calculado sobre aquela acima ). Y e X s são, portanto, padrão. As somas de quadrados observadas são agora = df .$\text{df}$$Y$$X$$\text{df}$

3. Calcular os coeficientes de predição regressivo por X s de acordo com o alvo r s: b = ( X ' x ) - 1 S .$Y$$X$$r$$\bf b=(X'X)^{-1} S$

4. Calcular os valores previstos para : Y = X b .$Y$$\hat{Y}=\bf Xb$

5. Calcular resíduos .$E=Y-\hat{Y}$

6. Calcular a soma necessária (alvo) de quadrados para resíduos: .$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}$$

   (o denominador não muda nas iterações, calcule-o com antecedência)

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

Eu estava com vontade de fazer alguma programação, então peguei a resposta excluída de @ Adam e decidi escrever uma boa implementação em R. Eu me concentro em usar um estilo orientado a funções (por exemplo, looping no estilo lapply). A idéia geral é pegar dois vetores, permutar aleatoriamente um dos vetores até que uma certa correlação seja alcançada entre eles. Essa abordagem é muito bruta, mas é simples de implementar.

Primeiro, criamos uma função que permite aleatoriamente o vetor de entrada:

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

... e crie alguns dados de exemplo

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

... escreva uma função que permita o vetor de entrada e a correlacione com um vetor de referência:

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

... e itere mil vezes:

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

Observe que as regras de escopo de R garantem isso vec1e vec2são encontradas no ambiente global, fora da função anônima usada acima. Portanto, todas as permutações são relativas aos conjuntos de dados de teste originais que geramos.

Em seguida, encontramos a correlação máxima:

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

... ou encontre o valor mais próximo de uma correlação de 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

Para obter uma correlação mais alta, você precisa aumentar o número de iterações.

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

Solução:

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

Código 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)))

Saída de teste:

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

Gere variáveis ​​normais com a matriz de covariância SAMPLING, conforme indicado

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;

Gere variáveis ​​normais com a matriz de covariância POPULATION, conforme indicado

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

Basta criar um vetor aleatório e classificar até obter o r desejado.