Menggunakan Distribusi yang Seragam untuk Menghasilkan Sampel Acak yang Berhubungan di R

8

[Pada pertanyaan baru-baru ini saya sedang mencari untuk menghasilkan vektor acak dalam R , dan saya ingin membagikan "penelitian" itu sebagai tanya jawab independen pada titik tertentu.]

Menghasilkan data acak dengan korelasi dapat dilakukan dengan menggunakan dekomposisi Cholesky dari matriks korelasi sini , sebagaimana tercermin pada posting sebelumnya di sini dan di sini .C=LLT

Pertanyaan yang saya ingin alamat adalah bagaimana menggunakan distribusi Uniform untuk menghasilkan angka acak berkorelasi dari distribusi marjinal berbeda di R .

Antoni Parellada
sumber
2
Anda tampaknya telah menemukan kembali Gaussian copula misalnya, lihat pertanyaan terkait di sini . Ada banyak kopula lain yang populer digunakan, tetapi Gaussian cukup nyaman dan bisa sangat cocok untuk beberapa situasi.
Glen_b -Reinstate Monica

Jawaban:

8

Karena pertanyaannya adalah

"cara menggunakan distribusi Seragam untuk menghasilkan angka acak berkorelasi dari distribusi marginal yang berbeda di "R

dan tidak hanya varian acak normal, jawaban di atas tidak menghasilkan simulasi dengan korelasi yang dimaksudkan untuk sepasang distribusi marginal sewenang-wenang di .R

Alasannya adalah bahwa, untuk sebagian besar cdf dan , when mana menunjukkan standar cdf normal.GXGY

cor(X,Y)cor(GX1(Φ(X),GY1(Φ(Y)),
(X,Y)N2(0,Σ),
Φ

Intinya, berikut adalah contoh tandingan dengan Exp (1) dan Gamma (.2,1) sebagai pasangan distribusi marginal saya di .R

library(mvtnorm)
#correlated normals with correlation 0.7
x=rmvnorm(1e4,mean=c(0,0),sigma=matrix(c(1,.7,.7,1),ncol=2),meth="chol")
cor(x[,1],x[,2])
  [1] 0.704503
y=pnorm(x) #correlated uniforms
cor(y[,1],y[,2])
  [1] 0.6860069
#correlated Exp(1) and Ga(.2,1)
cor(-log(1-y[,1]),qgamma(y[,2],shape=.2))
  [1] 0.5840085

Contoh tandingan lain yang jelas adalah ketika adalah cdf, dalam hal ini korelasinya tidak didefinisikan.GX

Untuk memberikan gambaran yang lebih luas, berikut ini adalah kode R di mana dan keduanya berubah-ubah:GXGY

etacor=function(rho=0,nsim=1e4,fx=qnorm,fy=qnorm){
  #generate a bivariate correlated normal sample
  x1=rnorm(nsim);x2=rnorm(nsim)
  if (length(rho)==1){
    y=pnorm(cbind(x1,rho*x1+sqrt((1-rho^2))*x2))
    return(cor(fx(y[,1]),fy(y[,2])))
    }
  coeur=rho
  rho2=sqrt(1-rho^2)
  for (t in 1:length(rho)){
     y=pnorm(cbind(x1,rho[t]*x1+rho2[t]*x2))
     coeur[t]=cor(fx(y[,1]),fy(y[,2]))}
  return(coeur)
  }

masukkan deskripsi gambar di sini

Bermain-main dengan berbagai cdf membuat saya memilih kasus khusus untuk dan distribusi log-Normal untuk :χ32GXGY

rhos=seq(-1,1,by=.01)
trancor=etacor(rho=rhos,fx=function(x){qchisq(x,df=3)},fy=qlnorm)
plot(rhos,trancor,ty="l",ylim=c(-1,1))
abline(a=0,b=1,lty=2)

yang menunjukkan seberapa jauh dari diagonal korelasinya.

Peringatan terakhir Mengingat dua distribusi sewenang-wenang dan , rentang nilai yang mungkin dari tidak harus . Masalahnya mungkin tidak punya solusi.GXGYcor(X,Y)(1,1)

Xi'an
sumber
Fantastis! Ty! Apakah ada cara agar kita dapat menemukan segmen perkiraan di mana keberangkatan tidak ditandai, seperti tampaknya normal, masih masuk akal untuk aplikasi praktis?
Antoni Parellada
5

Saya menulis correlatepaket. Orang-orang mengatakan itu menjanjikan (layak diterbitkan dalam Journal of Statistical Software), tetapi saya tidak pernah menulis makalah untuk itu karena saya memilih untuk tidak mengejar karir akademik.

Saya percaya correlatepaket yang tidak dikelola masih pada CRAN.

Ketika Anda menginstalnya, Anda dapat melakukan hal berikut:

require('correlate')
a <- rnorm(100)
b <- runif(100)
newdata <- correlate(cbind(a,b),0.5)

Hasilnya adalah bahwa data baru akan memiliki korelasi 0,5, tanpa mengubah distribusi univariat adan b(nilai-nilai yang sama ada, mereka hanya bergerak di sekitar sampai korelasi multivariat 0,5 telah tercapai.

Saya akan menjawab pertanyaan di sini, maaf karena kurangnya dokumentasi.

PascalVKooten
sumber
Bravo, ini jawaban yang sempurna! Apakah Anda memiliki cara untuk mendeteksi nilai-nilai korelasi yang tidak mungkin dijangkau?
Xi'an
@ Xi'an Ada beberapa kemustahilan, seperti beberapa titik data dan korelasi spesifik yang dicari yang tidak dapat dicapai. misalnya hanya memiliki 3 nilai berpasangan.
PascalVKooten
Perhatikan juga kemungkinan untuk lebih dari 2 variabel, misalnya untuk 3 variabel Anda dapat mendefinisikan matriks korelasi 3x3, 4 variabel a 4x4.
PascalVKooten
Secara umum itu akan bekerja selama Anda tidak ingin yang mustahil, tetapi sebelum Anda melakukan pekerjaan yang serius dengan itu disarankan untuk melakukan beberapa tes berjalan.
PascalVKooten
Orang-orang yang tertarik menggunakan data pendapatan; banyak nol dan distribusi gaussian-ish untuk pendapatan yang tidak nol.
PascalVKooten
1
  1. Hasilkan dua sampel data berkorelasi dari distribusi acak normal standar setelah korelasi yang telah ditentukan .

    Sebagai contoh, mari kita pilih korelasi r = 0,7 , dan kode matriks korelasi seperti:

    (C <- matrix(c(1,0.7,0.7,1), nrow = 2)) [,1] [,2] [1,] 1.0 0.7 [2,] 0.7 1.0

    Kita dapat menggunakan mvtnormuntuk menghasilkan sekarang dua sampel ini sebagai vektor acak bivariat:

    set.seed(0)

    SN <- rmvnorm(mean = c(0,0), sig = C, n = 1e5)menghasilkan dua komponen vektor yang didistribusikan sebagai ~ dan dengan a . Kedua komponen dapat diekstrak sebagai berikut:N(0,1)cor(SN[,1],SN[,2])= 0.6996197 ~ 0.7

    X1 <- SN[,1]; X2 <- SN[,2]

    Berikut plot dengan garis regresi yang tumpang tindih:

  2. Gunakan Probability Integral Transform di sini untuk mendapatkan vektor acak bivariat dengan distribusi marginal ~U(0,1) dan korelasi yang sama :

    U <- pnorm(SN)- jadi kami memasukkan ke pnormdalam SNvektor untuk menemukanerf(SN) (atau Φ(SN)). Dalam prosesnya, kami melestarikan cor(U[,1], U[,2]) = 0.6816123 ~ 0.7.

    Sekali lagi kita dapat menguraikan vektor U1 <- U[,1]; U2 <- U[,2]dan menghasilkan sebar dengan distribusi marjinal di tepi, jelas menunjukkan sifat seragam mereka:

  3. Terapkan metode sampling transformasi terbalik di sini untuk akhirnya mendapatkan bivektor dari titik berkorelasi sama yang dimiliki keluarga distribusi mana pun yang kami tuju untuk mereproduksi.

    Dari sini kita hanya dapat menghasilkan dua vektor yang terdistribusi secara normal dan dengan varian yang sama atau berbeda . Misalnya: Y1 <- qnorm(U1, mean = 8,sd = 10)dan Y2 <- qnorm(U2, mean = -5, sd = 4), yang akan mempertahankan korelasi yang diinginkan cor(Y1,Y2) = 0.6996197 ~ 0.7,.

    Atau pilih distribusi yang berbeda. Jika distribusi yang dipilih sangat berbeda, korelasinya mungkin tidak tepat. Sebagai contoh, mari kita U1ikuti atdistribusi dengan 3 df, dan U2eksponensial dengan aλ= 1: Z1 <- qt(U1, df = 3)dan Z2 <- qexp(U2, rate = 1)The cor(Z1,Z2) [1] 0.5941299 < 0.7. Berikut adalah histogram masing-masing:

Berikut adalah contoh kode untuk seluruh proses dan marginal normal:

Cor_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
require(mvtnorm)
SN <- rmvnorm(mean = c(0,0), sig = C, n = n)
U <- pnorm(SN)
U1 <- U[,1]
U2 <- U[,2]

 Y1 <<- qnorm(U1, mean = mean1,sd = sd1) 
 Y2 <<- qnorm(U2, mean = mean2,sd = sd2) 

sample_measures <<- as.data.frame(c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1,Y2)), names<-c("mean Y1", "mean Y2", "SD Y1", "SD Y2", "Cor(Y1,Y2)"))
sample_measures
}

Sebagai perbandingan, saya telah mengumpulkan fungsi berdasarkan dekomposisi Cholesky:

Cholesky_samples <- function(r, n, mean1, mean2, sd1, sd2){
C <- matrix(c(1,r,r,1), nrow = 2)
L <- chol(C)
X1 <- rnorm(n)
X2 <- rnorm(n)
X <- rbind(X1,X2)

Y <- t(L)%*%X
Y1 <- Y[1,]
Y2 <- Y[2,]

N_1 <<- Y[1,] * sd1 + mean1
N_2 <<- Y[2,] * sd2 + mean2

sample_measures <<- as.data.frame(c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2)), 
                  names<-c("mean N_1", "mean N_2", "SD N_1", "SD N_2","cor(N_1,N_2)"))
sample_measures
}

Mencoba kedua metode untuk menghasilkan yang berkorelasi (katakanlah, r=0.7) sampel didistribusikan ~ N(97,23) dan N(32,8)kita dapatkan, pengaturan set.seed(99):

Menggunakan Seragam:

cor_samples(0.7, 1000, 97, 32, 23, 8)
           c(mean(Y1), mean(Y2), sd(Y1), sd(Y2), cor(Y1, Y2))
mean Y1                                            96.5298821
mean Y2                                            32.1548306
SD Y1                                              22.8669448
SD Y2                                               8.1150780
cor(Y1,Y2)                                          0.7061308

dan Menggunakan Cholesky:

Cholesky_samples(0.7, 1000, 97, 32, 23, 8)
             c(mean(N_1), mean(N_2), sd(N_1), sd(N_2), cor(N_1, N_2))
mean N_1                                                   96.4457504
mean N_2                                                   31.9979675
SD N_1                                                     23.5255419
SD N_2                                                      8.1459100
cor(N_1,N_2)                                                0.7282176
Antoni Parellada
sumber
Secara empiris, sepertinya Anda beralih dari N (0,1) ->
F1(X)
~ Unif. ->
f(F1(X))
~ terdistribusi sesuai dengan distribusi yang dipilih, korelasinya tidak berubah kecuali distribusi terakhir secara substansial berbeda dari N awal (0,1). Saya menyertakan nilai-nilai ... Dalam kasus apa pun, apakah Anda melihat masalah khusus dengan metode itu sendiri untuk aplikasi praktis?
Antoni Parellada
Saya mengubah fungsi di akhir jawaban untuk memasukkan korelasi sampel yang dihitung, sehingga dibandingkan dengan angka yang terpasang, dan mereka tampaknya cocok.
Antoni Parellada
2
Apakah ada masalah dengan aplikasi praktis tergantung pada aplikasi praktis; untuk beberapa hal ini tidak masalah. Perhatikan bahwa karena transformasi bersifat monoton, korelasi nonparametrik seperti Spearman rho dan Kendall's tau tidak akan berubah.
Glen_b -Reinstate Monica