Hasilkan variabel acak dengan korelasi yang ditentukan dengan variabel yang ada

Untuk studi simulasi saya harus membuat variabel acak yang menunjukkan korelasi (populasi) prefined ke variabel .$Y$

Saya melihat ke dalam Rpaket copuladan CDVineyang dapat menghasilkan distribusi multivarian acak dengan struktur ketergantungan yang diberikan. Namun, tidak mungkin untuk memperbaiki salah satu variabel yang dihasilkan ke variabel yang ada.

Setiap ide dan tautan ke fungsi yang ada dihargai!

Kesimpulan: Dua jawaban yang valid muncul, dengan solusi yang berbeda:

1. Sebuah R naskah oleh caracal, yang menghitung variabel acak dengan tepat (sampel) korelasi untuk variabel yang telah ditetapkan
2. Suatu R fungsi yang saya temukan sendiri, yang menghitung variabel acak dengan korelasi populasi yang ditentukan dengan variabel yang telah ditentukan

[@ttnphns 'tambahan: Saya mengambil kebebasan untuk memperluas judul pertanyaan dari kasus variabel tetap tunggal ke jumlah variabel tetap sewenang-wenang; yaitu cara membuat variabel yang memiliki corretation yang telah ditentukan sebelumnya dengan beberapa variabel tetap yang ada]

Ini satu lagi: untuk vektor dengan rata-rata 0, korelasinya sama dengan kosinus sudutnya. Jadi salah satu cara untuk menemukan vektor dengan korelasi yang diinginkan , sesuai dengan sudut :r $x$$r$$\theta$

1. dapatkan vektor tetap dan vektor acakx $x_1$$x_2$
2. pusatkan kedua vektor (rata-rata 0), memberikan vektor , ˙ x$\dot{x}_{1}$$\dot{x}_{2}$
3. buat orthogonal ke (proyeksi ke subruang ortogonal), memberikan ˙ x 1 ˙ x ⊥ $\dot{x}_{2}$$\dot{x}_{1}$$\dot{x}_{2}^{\perp}$
4. skala dan dengan panjang 1, memberikan dan ˙ 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}$ adalah vektor yang sudutnya ke adalah , dan yang hubungannya dengan dengan demikian adalah . Ini juga merupakan korelasi dengan karena transformasi linear membuat korelasinya tidak berubah.$\bar{x}_{1}$$\theta$$\bar{x}_{1}$$r$$x_1$

Ini kodenya:

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

Untuk proyeksi ortogonal , saya menggunakan dekomposisi untuk meningkatkan stabilitas numerik, sejak saat itu cukup .Q R P = Q Q $P$$QR$$P = Q Q'$

Saya akan menjelaskan solusi yang paling umum. Memecahkan masalah dalam generalitas ini memungkinkan kita untuk mencapai implementasi perangkat lunak yang sangat kompak: cukup dua baris Rkode saja.

Pilih vektor , dengan panjang yang sama dengan , sesuai dengan distribusi yang Anda suka. Mari menjadi residual kuadrat regresi setidaknya dari terhadap : ini ekstrak komponen dari . Dengan menambahkan kembali kelipatan cocok ke , kita dapat menghasilkan vektor memiliki apapun yang diinginkan korelasi dengan . Hingga konstanta multiplikatif aditif sewenang-wenang dan positif - yang bebas Anda pilih dengan cara apa pun - solusinya adalahY Y ⊥ X Y Y X Y Y ⊥ ρ $X$$Y$$Y^\perp$$X$$Y$$Y$$X$$Y$$Y^\perp$$\rho$$Y$

(" " adalah singkatan dari setiap perhitungan yang sebanding dengan standar deviasi.)$\operatorname{SD}$

Ini Rkode kerjanya . Jika Anda tidak menyediakan , kode akan mengambil nilainya dari distribusi Normal standar multivariat.$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)
}

Sebagai ilustrasi, saya membuat acak dengan komponen dan menghasilkan memiliki berbagai korelasi spesifik dengan ini . Mereka semua dibuat dengan vektor awal yang sama . Berikut adalah sebar plot mereka. "Rugplots" di bagian bawah setiap panel menunjukkan vektor umum .50 X Y ; ρ Y X = ( 1 , 2 , … , 50 ) $Y$$50$$X_{Y;\rho}$$Y$$X=(1,2,\ldots, 50)$$Y$

Ada kesamaan yang luar biasa di antara plot, tidak ada :-).

Jika Anda ingin bereksperimen, berikut adalah kode yang menghasilkan data ini dan gambar. (Saya tidak repot-repot menggunakan kebebasan untuk mengubah dan mengukur hasilnya, yang merupakan operasi yang mudah.)

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, metode ini siap digeneralisasi menjadi lebih dari satu : jika secara matematis memungkinkan, ia akan menemukan memiliki korelasi yang ditentukan dengan keseluruhan set . Cukup gunakan kuadrat terkecil biasa untuk menghilangkan efek semua dari dan membentuk kombinasi linear yang sesuai dari dan residu. (Ini membantu untuk melakukan ini dalam hal basis ganda untuk , yang diperoleh dengan menghitung pseudo-invers. Kode follownig menggunakan SVD untuk mencapai itu.)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$

Berikut ini sketsa algoritme dalam R, di mana diberikan sebagai kolom dari sebuah matriks :$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)

Berikut ini adalah implementasi yang lebih lengkap bagi mereka yang ingin bereksperimen.

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

Berikut pendekatan komputasi lain (solusinya diadaptasi dari posting forum oleh Enrico Schumann). Menurut Wolfgang (lihat komentar), ini identik secara komputasi dengan solusi yang diajukan oleh ttnphns.

Berbeda dengan solusi caracal, ia tidak menghasilkan sampel dengan korelasi tepat , tetapi dua vektor yang populasinya korelasinya sama dengan .$\rho$$\rho$

Fungsi berikut dapat menghitung distribusi sampel bivariat yang diambil dari suatu populasi dengan diberikan . Entah menghitung dua variabel acak, atau mengambil satu variabel yang ada (dilewatkan sebagai parameter ) dan membuat variabel kedua dengan korelasi yang diinginkan:$\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)
}

Fungsi ini juga dapat menggunakan distribusi marginal non-normal dengan menyesuaikan parameter mar.fun. Namun, perlu diketahui bahwa memperbaiki satu variabel hanya berfungsi dengan variabel yang terdistribusi normal x! (yang mungkin terkait dengan komentar Makro).

Juga perhatikan bahwa "faktor koreksi kecil" dari pos asli telah dihapus karena tampaknya bias korelasi yang dihasilkan, setidaknya dalam kasus distribusi Gaussian dan korelasi Pearson (juga lihat komentar).

Biarkan menjadi variabel tetap Anda dan Anda ingin menghasilkan variabel yang berkorelasi dengan dengan jumlah . Jika distandarisasi maka (karena adalah koefisien beta dalam regresi sederhana) , di mana adalah variabel acak dari distribusi normal yang memiliki mean dan . Korelasi yang diamati antara data dan akan kira-kira ; dan dapat dilihat sebagai sampel acak dari populasi normal bivariat (jikaY X r X r Y = r X + E E 0 sd = $X$$Y$$X$$r$$X$$r$$Y= rX+E$$E$$0$ XYrXYXρ=$\text{sd}=\sqrt{1-r^2}$$X$$Y$$r$$X$$Y$$X$ dari normal) dengan .$\rho=r$

Sekarang, jika Anda ingin mencapai korelasi dalam sampel bivariat Anda persis , Anda perlu memberikan yang memiliki nol korelasi dengan . Pengetatan ini menjadi nol dapat dicapai dengan memodifikasi secara berulang. Nah, dengan hanya dua variabel, satu diberikan ( ) dan satu untuk menghasilkan ( ), jumlah iterasi yang cukup sebenarnya 1, tetapi dengan beberapa variabel yang diberikan ( ) iterasi akan diperlukan.E X E X Y X 1 , X 2 , X 3 , . . $r$$E$$X$$E$$X$$Y$$X_1, X_2, X_3,...$

Perlu dicatat bahwa jika normal maka pada prosedur pertama ("perkiraan ") juga akan normal; Namun, dalam pemasangan berulang ke "tepat " cenderung kehilangan normalitas karena pemasangan mengeksploitasi nilai kasus secara selektif.r Y Y r $X$$r$$Y$$Y$$r$$Y$

Pembaruan 11 Nov 2017. Saya telah menemukan utas lama ini hari ini dan memutuskan untuk memperluas jawaban saya dengan menunjukkan algoritme pengulangan yang pas tentang yang saya bicarakan pada awalnya.

$Y$ $X$

Disclamer: Ini solusi berulang yang saya temukan lebih rendah daripada yang terbaik berdasarkan menemukan basis ganda dan diusulkan oleh @whuber di utas ini hari ini. solusi @ whuber tidak iteratif dan, yang lebih penting bagi saya, tampaknya akan mempengaruhi nilai-nilai dari variabel input "babi" agak kurang dari algoritma "saya" (itu akan menjadi aset kemudian jika tugasnya adalah untuk "memperbaiki" variabel yang ada dan tidak menghasilkan variate acak dari awal). Tetap saja, saya menerbitkan buku saya untuk rasa ingin tahu dan karena itu berhasil (lihat juga Catatan Kaki).

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

   (penyebut tidak berubah pada iterasi, hitung terlebih dahulu)

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

Saya merasa ingin melakukan beberapa pemrograman, jadi saya mengambil jawaban yang dihapus @ Adam dan memutuskan untuk menulis implementasi yang bagus di R. Saya fokus menggunakan gaya berorientasi fungsional (yaitu lapply style looping). Gagasan umum adalah untuk mengambil dua vektor, secara acak mengubah salah satu vektor sampai korelasi tertentu telah tercapai di antara mereka. Pendekatan ini sangat kasar, tetapi mudah diimplementasikan.

Pertama kita membuat fungsi yang secara acak mengizinkan vektor input:

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

... dan buat beberapa contoh data

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

... tulis fungsi yang memungkinkan vektor input, dan menghubungkannya ke vektor referensi:

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

... dan ulangi ribuan kali:

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

Perhatikan bahwa aturan pelingkupan R memastikan vec1dan vec2ditemukan di lingkungan global, di luar fungsi anonim yang digunakan di atas. Jadi, permutasi semua relatif terhadap dataset uji asli yang kami hasilkan.

Berikutnya, kami menemukan korelasi maksimum:

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

... atau temukan nilai terdekat dengan korelasi 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

Untuk mendapatkan korelasi yang lebih tinggi, Anda perlu menambah jumlah iterasi.

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

Larutan:

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

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

Hasil tes:

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

Hasilkan variabel normal dengan matriks kovarians SAMPLING seperti yang diberikan

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;

Hasilkan variabel normal dengan matriks kovarians POPULASI seperti yang diberikan

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

Cukup buat vektor acak dan urutkan sampai Anda mendapatkan r yang diinginkan.