Bagaimana cara mensimulasikan data yang memenuhi kendala spesifik seperti memiliki rata-rata spesifik dan standar deviasi?

56

Pertanyaan ini dimotivasi oleh pertanyaan saya tentang meta-analisis . Tapi saya membayangkan bahwa itu juga akan berguna dalam konteks pengajaran di mana Anda ingin membuat dataset yang persis mencerminkan dataset yang sudah ada diterbitkan.

Saya tahu cara menghasilkan data acak dari distribusi yang diberikan. Jadi misalnya, jika saya membaca tentang hasil studi yang memiliki:

  • rata-rata 102,
  • standar deviasi 5,2, dan
  • ukuran sampel 72.

Saya bisa menghasilkan data serupa menggunakan rnormdalam R. Misalnya,

set.seed(1234)
x <- rnorm(n=72, mean=102, sd=5.2)

Tentu saja mean dan SD tidak akan persis sama dengan 102 dan 5.2 masing-masing:

round(c(n=length(x), mean=mean(x), sd=sd(x)), 2)
##     n   mean     sd 
## 72.00 100.58   5.25 

Secara umum saya tertarik pada bagaimana mensimulasikan data yang memenuhi serangkaian kendala. Dalam kasus di atas, kendala adalah ukuran sampel, rata-rata, dan standar deviasi. Dalam kasus lain, mungkin ada kendala tambahan. Sebagai contoh,

  • minimum dan maksimum dalam data atau variabel yang mendasarinya mungkin diketahui.
  • variabel mungkin diketahui hanya mengambil nilai integer atau hanya nilai-nilai non-negatif.
  • data mungkin mencakup banyak variabel dengan inter-korelasi yang diketahui.

Pertanyaan

  • Secara umum, bagaimana saya bisa mensimulasikan data yang benar-benar memenuhi serangkaian kendala?
  • Apakah ada artikel yang ditulis tentang ini? Apakah ada program di R yang melakukan ini?
  • Sebagai contoh, bagaimana saya bisa mensimulasikan variabel sehingga memiliki rata-rata dan sd tertentu?
Jeromy Anglim
sumber
1
Mengapa Anda ingin mereka persis seperti hasil yang dipublikasikan? Bukankah ini perkiraan rata-rata populasi dan standar deviasi yang diberikan sampel data mereka. Mengingat ketidakpastian dalam estimasi tersebut, siapa yang dapat mengatakan bahwa sampel yang Anda perlihatkan di atas tidak konsisten dengan pengamatan mereka?
Pasang kembali Monica - G. Simpson
4
Karena pertanyaan ini tampaknya mengumpulkan balasan yang tidak sesuai dengan standar (IMHO), saya ingin menunjukkan bahwa secara konseptual jawabannya langsung: kendala kesetaraan diperlakukan seperti distribusi marjinal dan kendala ketidaksetaraan adalah analog pemotongan multivarian. Pemotongan relatif mudah ditangani (seringkali dengan sampel penolakan); masalah yang lebih sulit adalah menemukan cara untuk mencicipi distribusi marjinal ini. Ini berarti baik sampel marjinal yang diberikan distribusi dan kendala, atau berintegrasi untuk menemukan distribusi marginal dan pengambilan sampel dari itu.
whuber
4
BTW, pertanyaan terakhir adalah sepele untuk keluarga distribusi skala lokasi. Misalnya, x<-rnorm(72);x<-5.2*(x-mean(x))/sd(x)+102lakukan triknya.
whuber
1
@whuber, sebagai kardinal menyinggung dalam komentar untuk jawaban saya (yang menyebutkan "trik" ini) dan komentar untuk jawaban lain - metode ini, secara umum, tidak akan menyimpan variabel dalam keluarga distribusi yang sama, karena Anda membagi oleh standar deviasi sampel.
Makro
5
@ Macro Ini adalah poin yang bagus, tetapi mungkin jawaban terbaik adalah, "tentu saja mereka tidak akan memiliki distribusi yang sama"! Distribusi yang Anda inginkan adalah distribusi yang tergantung pada kendala. Secara umum itu tidak akan berasal dari keluarga yang sama dengan distribusi orang tua. Misalnya, setiap elemen sampel berukuran 4 dengan rata-rata 0 dan SD 1 yang diambil dari distribusi normal akan memiliki probabilitas yang hampir seragam pada [-1,5, 1,5], karena kondisi menempatkan batas atas dan bawah pada nilai yang mungkin.
Whuber

Jawaban:

26

Secara umum, untuk membuat mean dan varians sampel Anda persis sama dengan nilai yang ditentukan sebelumnya, Anda dapat dengan tepat menggeser dan menskalakan variabel. Secara khusus, jika adalah sampel, maka variabel baruX1,X2,...,Xn

Zsaya=c1(Xsaya-X¯sX)+c2

di mana adalah mean sampel dans 2 X =1X¯=1nsaya=1nXsayaadalah varians sampel adalah seperti bahwa mean sampel dariZi's persisc2dan varians sampel mereka adalah persisc1. Contoh yang dibuat serupa dapat membatasi rentang -sX2=1n-1saya=1n(Xsaya-X¯)2Zsayac2c1

Bsaya=Sebuah+(b-Sebuah)(Xsaya-min({X1,...,Xn})maks({X1,...,Xn})-min({X1,...,Xn}))

akan menghasilkan kumpulan data yang terbatas pada interval ( a , b ) . B1,...,Bn(Sebuah,b)

Catatan: Jenis pemindahan / penskalaan ini akan, secara umum, mengubah keluarga distribusi data, bahkan jika data asli berasal dari keluarga skala lokasi.

Dalam konteks distribusi normal , mvrnormfungsi in R memungkinkan Anda untuk mensimulasikan data normal (atau multivarian normal) dengan sampel rata-rata / kovarians sampel dengan menetapkan empirical=TRUE. Secara khusus, fungsi ini mensimulasikan data dari distribusi kondisional dari variabel yang terdistribusi normal, mengingat rerata sampel dan (co) varians sama dengan nilai yang ditentukan sebelumnya . Perhatikan bahwa distribusi marjinal yang dihasilkan tidak normal, seperti yang ditunjukkan oleh @whuber dalam komentar untuk pertanyaan utama.

Berikut adalah contoh univariat sederhana di mana mean sampel (dari sampel ) dibatasi menjadi 0 dan standar deviasi sampel adalah 1. Kita dapat melihat bahwa elemen pertama jauh lebih mirip dengan distribusi seragam daripada normal. distribusi:n=4

library(MASS)
 z = rep(0,10000)
for(i in 1:10000)
{
    x = mvrnorm(n = 4, rep(0,1), 1, tol = 1e-6, empirical = TRUE)
    z[i] = x[1]
}
hist(z, col="blue")

                  masukkan deskripsi gambar di sini

Makro
sumber
1
The akan tidak terdistribusi secara normal, meskipun mereka mungkin kurang jadi jika ukuran sampel besar. Komentar pertama untuk jawaban @ Sean menyinggung ini. Zsaya
kardinal
1
Ya, itu hal yang wajar untuk dilakukan ... dan seringkali tidak menyebabkan terlalu banyak masalah.
kardinal
1
+1. Dalam contoh, seragam adalah jawaban yang tepat . (Penurunan yang jelas di ujung plot adalah artefak tentang bagaimana R menggambar histogram.)
whuber
1
@whuber, terima kasih telah memotivasi contoh ini. Mengingat fakta bahwa distribusi marjinal berubah begitu Anda mengkondisikan pada mean / varians sampel, tampaknya "jawaban" terbaik dalam semangat pertanyaan OP adalah hanya mensimulasikan data dengan mean / varians populasi yang sama dengan yang dilaporkan sebagai sampel. kuantitas (seperti yang disarankan oleh OP sendiri), bukan? Dengan begitu, Anda mendapatkan jumlah sampel "mirip" dengan yang diinginkan, dan distribusi marginal adalah apa yang Anda inginkan.
Makro
1
@whuber, Jika sampel Anda normal, maka memiliki distribusi- t , ya? Variabel "baru" yang dipermasalahkan hanya akan menjadi kombinasi linear dari T i . Tsaya=(Xsaya-X¯)/stTsaya
Makro
22

Mengenai permintaan Anda untuk makalah, ada:

Ini bukan apa yang Anda cari, tetapi mungkin berfungsi sebagai gandum bagi pabrik.


Ada strategi lain yang tampaknya tidak ada yang disebutkan. Dimungkinkan untuk menghasilkan data acak (pseudo) dari satu set ukuran N sehingga seluruh set memenuhi kendala k selama data k yang tersisa ditetapkan pada nilai yang sesuai. Nilai yang diperlukan harus dipecahkan dengan sistem persamaan k , aljabar, dan beberapa minyak siku. N-kNkkk

Misalnya, untuk menghasilkan satu set data dari distribusi normal yang akan memiliki rata-rata sampel yang diberikan, ˉ xNx¯ , dan varians, , Anda harus memperbaiki nilai dua titik: y dan z . Karena rata-rata sampel adalah: ˉ x = Σ N - 2 i = 1 x is2yz
yharus: y=N ˉ x

x¯=saya=1N-2xsaya+y+zN
y Varians sampel: s 2 = Σ N - 2 i = 1 ( x i - ˉ x ) 2
y=Nx¯-(saya=1N-2xsaya+z)
demikian (setelah menggantiydi atas, menggagalkan / mendistribusikan, & menata ulang ...) kita mendapatkan: 2(N ˉ x
s2=saya=1N-2(xsaya-x¯)2+(y-x¯)2+(z-x¯)2N-1
y Jika kita mengambil a = - 2 , b = 2 ( N ˉ x - ¢
2(Nx¯-saya=1N-2xsaya)z-2z2=Nx¯2(N-1)+saya=1N-2xsaya2+[saya=1N-2xsaya]2-2Nx¯saya=1N-2xsaya-(N-1)s2
Sebuah=-2, dancsebagai negasi dari RHS, kita dapat menyelesaikan untukzmenggunakanrumus kuadratik. Misalnya, dalam, kode berikut dapat digunakan: b=2(Nx¯-saya=1N-2xsaya)czR
find.yz = function(x, xbar, s2){
  N    = length(x) + 2
  sumx = sum(x)
  sx2  = as.numeric(x%*%x)          # this is the sum of x^2
  a    = -2
  b    = 2*(N*xbar - sumx)
  c    = -N*xbar^2*(N-1) - sx2 - sumx^2 + 2*N*xbar*sumx + (N-1)*s2
  rt   = sqrt(b^2 - 4*a*c)

  z    = (-b + rt)/(2*a)
  y    = N*xbar - (sumx + z)
  newx = c(x, y, z)
  return(newx)
}

set.seed(62)
x    = rnorm(2)
newx = find.yz(x, xbar=0, s2=1)
newx                                # [1] 0.8012701  0.2844567  0.3757358 -1.4614627
mean(newx)                          # [1] 0
var(newx)                           # [1] 1

N-2yzs2

set.seed(22)    
x    = rnorm(2)
newx = find.yz(x, xbar=0, s2=1)
Warning message:
In sqrt(b^2 - 4 * a * c) : NaNs produced
newx                                # [1] -0.5121391  2.4851837        NaN        NaN
var(c(x, mean(x), mean(x)))         # [1] 1.497324

Kedua, sedangkan standardisasi membuat distribusi marginal dari semua varian Anda lebih seragam, pendekatan ini hanya memengaruhi dua nilai terakhir, tetapi membuat distribusi marginalnya miring:

set.seed(82)
xScaled = matrix(NA, ncol=4, nrow=10000)
for(i in 1:10000){
  x           = rnorm(4)
  xScaled[i,] = scale(x)
}

(sisipkan plot)

set.seed(82)
xDf = matrix(NA, ncol=4, nrow=10000)
i   = 1
while(i<10001){
  x       = rnorm(2)
  xDf[i,] = try(find.yz(x, xbar=0, s2=2), silent=TRUE)  # keeps the code from crashing
  if(!is.nan(xDf[i,4])){ i = i+1 }                      # increments if worked
}

(sisipkan plot)

Ketiga, sampel yang dihasilkan mungkin tidak terlihat sangat normal; mungkin terlihat seperti memiliki 'outlier' (yaitu, poin yang berasal dari proses pembuatan data yang berbeda dari yang lain), karena pada dasarnya itulah masalahnya. Ini cenderung menjadi masalah dengan ukuran sampel yang lebih besar, karena statistik sampel dari data yang dihasilkan harus menyatu dengan nilai-nilai yang diperlukan dan dengan demikian membutuhkan penyesuaian yang lebih sedikit. Dengan sampel yang lebih kecil, Anda selalu bisa menggabungkan pendekatan ini dengan algoritma terima / tolak yang mencoba lagi jika sampel yang dihasilkan memiliki statistik bentuk (misalnya, skewness dan kurtosis) yang berada di luar batas yang dapat diterima (lih., Komentar @ cardinal ), atau memperluas pendekatan ini untuk menghasilkan sampel dengan mean tetap, varians, skewness, dankurtosis (saya akan membiarkan aljabar terserah Anda). Atau, Anda dapat menghasilkan sejumlah kecil sampel dan menggunakan sampel dengan statistik Kolmogorov-Smirnov terkecil.

library(moments)
set.seed(7900)  
x = rnorm(18)
newx.ss7900 = find.yz(x, xbar=0, s2=1)
skewness(newx.ss7900)                       # [1] 1.832733
kurtosis(newx.ss7900) - 3                   # [1] 4.334414
ks.test(newx.ss7900, "pnorm")$statistic     # 0.1934226

set.seed(200)  
x = rnorm(18)
newx.ss200 = find.yz(x, xbar=0, s2=1)
skewness(newx.ss200)                        # [1] 0.137446
kurtosis(newx.ss200) - 3                    # [1] 0.1148834
ks.test(newx.ss200, "pnorm")$statistic      # 0.1326304 

set.seed(4700)  
x = rnorm(18)
newx.ss4700 = find.yz(x, xbar=0, s2=1)
skewness(newx.ss4700)                       # [1]  0.3258491
kurtosis(newx.ss4700) - 3                   # [1] -0.02997377
ks.test(newx.ss4700, "pnorm")$statistic     # 0.07707929S

(tambahkan plot)

gung - Reinstate Monica
sumber
10

Teknik umum adalah 'Metode Penolakan', di mana Anda hanya menolak hasil yang tidak memenuhi kendala Anda. Kecuali jika Anda memiliki semacam panduan (seperti MCMC), maka Anda bisa menghasilkan banyak kasus (tergantung pada skenario Anda) yang ditolak!

Di mana Anda mencari sesuatu seperti mean dan deviasi standar dan Anda dapat membuat metrik jarak semacam untuk mengatakan seberapa jauh Anda jauh dari tujuan Anda, Anda dapat menggunakan optimasi untuk mencari variabel input yang memberikan Anda output yang diinginkan nilai-nilai.

Sebagai contoh jelek di mana kita akan mencari vektor seragam acak dengan panjang 100 yang memiliki mean = 0 dan standar deviasi = 1.

# simplistic optimisation example
# I am looking for a mean of zero and a standard deviation of one
# but starting from a plain uniform(0,1) distribution :-)
# create a function to optimise
fun <- function(xvec, N=100) {
  xmin <- xvec[1]
  xmax <- xvec[2]
  x <- runif(N, xmin, xmax)
  xdist <- (mean(x) - 0)^2 + (sd(x) - 1)^2
  xdist
}
xr <- optim(c(0,1), fun)

# now lets test those results
X <- runif(100, xr$par[1], xr$par[2])
mean(X) # approx 0
sd(X)   # approx 1
Sean
sumber
7
Kendala yang terjadi dengan probabilitas nol sulit untuk dipenuhi. ;-) Untuk contoh spesifik yang ada, pergeseran dan pelebaran yang tepat dengan mudah mencapai tujuan yang dinyatakan , meskipun orang mungkin ingin menganalisis sedikit lebih dalam untuk melihat bagaimana distribusi data terganggu oleh operasi seperti itu.
kardinal
Terima kasih. Tentu akan mudah untuk menolak pengamatan yang kurang dari minimum dan lebih besar dari maks. Dan saya dapat melihat bagaimana Anda dapat mendefinisikannya sebagai masalah optimisasi. Akan menyenangkan untuk melihat beberapa contoh atau mungkin memiliki beberapa saran tentang apa yang harus dibaca selanjutnya.
Jeromy Anglim
1
@ cardinal - setuju. Kita harus melihat distribusi (yaitu histogram) dari bilangan input dan output simulasi karena kadang-kadang ini bisa terlihat sangat aneh!
Sean
9

Apakah ada program di R yang melakukan ini?

Itu Runuran R berisi banyak metode untuk menghasilkan variasi acak. Ia menggunakan C library dari proyek UNU.RAN (Universal Non-Uniform RAndom Number generator) proyek. Pengetahuan saya sendiri tentang bidang pembuatan varian acak terbatas, tetapi sketsa Runuran memberikan gambaran yang bagus. Di bawah ini adalah metode yang tersedia dalam paket Runuran, diambil dari sketsa:

Distribusi berkelanjutan:

  • Sampling Penolakan Adaptif
  • Penolakan Kepadatan Transformasi Terbalik
  • Interpolasi Polinomial dari Inversi CDF
  • Metode Rasio-Seragam yang Sederhana
  • Penolakan Kepadatan Berubah

Distribusi terpisah:

  • Inversi Penolakan Otomatis Terpisah
  • Metode Alias-Urn
  • Panduan-Tabel Metode untuk Inversi Diskrit

Distribusi multivarian:

  • Algoritma Hit-and-Run dengan Metode Ratio-of-Uniforms
  • Metode Rasio-of-Seragam Naif Multivariat

Contoh:

Sebagai contoh cepat, misalkan Anda ingin menghasilkan distribusi normal yang dibatasi antara 0 dan 100:

require("Runuran")

## Normal distribution bounded between 0 and 100
d1 <- urnorm(n = 1000, mean = 50, sd = 25, lb = 0, ub = 100)

summary(d1)
sd(d1)
hist(d1)

The urnorm()fungsi fungsi pembungkus nyaman. Saya percaya bahwa di balik layar menggunakan Interpolasi Polinomial metode Invers CDF tetapi saya tidak yakin. Untuk sesuatu yang lebih kompleks, katakanlah, distribusi Normal diskrit dibatasi antara 0 dan 100:

require("Runuran")

## Discrete normal distribution bounded between 0 and 100
# Create UNU.RAN discrete distribution object
discrete <- unuran.discr.new(pv = dnorm(0:100, mean = 50, sd = 25), lb = 0, ub = 100)

# Create UNU.RAN object using the Guide-Table Method for Discrete Inversion
unr <- unuran.new(distr = discrete, method = "dgt")

# Generate random variates from the UNU.RAN object
d2 <- ur(unr = unr, n = 1000)

summary(d2)
sd(d2)
head(d2)
hist(d2)
jthetzel
sumber
3

Tampaknya ada paket R yang memenuhi kebutuhan Anda yang baru diterbitkan kemarin! simstudy Oleh Keith Goldfeld

Mensimulasikan set data untuk mengeksplorasi teknik pemodelan atau lebih memahami proses pembuatan data. Pengguna menentukan serangkaian hubungan antara kovariat, dan menghasilkan data berdasarkan spesifikasi ini. Set data akhir dapat mewakili data dari uji kontrol acak, desain ukuran berulang (longitudinal), dan uji coba acak kelompok. Hilangnya dapat dihasilkan menggunakan berbagai mekanisme (MCAR, MAR, NMAR).

Tyelcie
sumber
1
Baik dalam sketsa maupun di halaman beranda program pertemuan tepat kendala disebutkan. Menurut Anda mengapa paket ini memenuhi persyaratan menggambar dari distribusi bersyarat?
gg
2

saya=1nf(xsaya)
saya=1nxsaya=μ0saya=1nxsaya2=σ02
Xi'an
sumber
2

Jawaban ini mempertimbangkan pendekatan lain untuk kasus di mana Anda ingin memaksa variates untuk berada dalam rentang yang ditentukan dan tambahan menentukan rata-rata dan / atau varians.

[0,1]wk[0,1]k=1Nwk=1wk=1/Nμ(0,1)0<σ2<μ(1-μ)σ2x1,...,xN[0,1]

y1,...,yNN(0,1)[0,1]

xk=11+e-(ykv-h)

ykhvhvx1,...,xN

μ=k=1Nwk1+e-(ykv-h)σ2=k=1Nwk(1+e-(ykv-h))2-(k=1Nwk1+e-(ykv-h))2

vhvh

μ=0.8v=1wk=1/NN=200000N(0,1)N(0,0,1)Unif(0,1) 0.8N

Contoh 1

μ=0,2wk=1/NN=2000σ=0,1,0,05,0,01N(0,1)

Contoh 2

Perhatikan bahwa ini mungkin terlihat sedikit beta, tetapi sebenarnya tidak.

Ian Hincks
sumber
1

Dalam jawaban saya di sini , saya mencantumkan tiga paket R untuk melakukan ini:

abalter
sumber
Harus ada beberapa format untuk tautan ke referensi. Haruskah itu menjadi komentar?
Abalter