Bagaimana cara menemukan interval kredibel 95%?

13

Saya mencoba menghitung interval kredibel 95% dari distribusi posterior berikut. Saya tidak dapat menemukan fungsi di R untuk itu tetapi apakah pendekatan di bawah ini benar?

x <- seq(0.4,12,0.4)
px <-  c(0,0, 0, 0, 0, 0, 0.0002, 0.0037, 0.018, 0.06, 0.22 ,0.43, 0.64,0.7579, 0.7870, 0.72, 0.555, 0.37, 0.24, 0.11, 0.07, 0.02, 0.009, 0.005, 0.0001, 0,0.0002, 0, 0, 0)
plot(x,px, type="l")
mm <- sum(x*px)/sum(px)
var <- (sum((x)^2*px)/sum(px)) - (mm^2)
cat("95% credible interval: ", round(mm -1.96*sqrt(var),3), "-", round(mm + 1.96*sqrt(var),3),"\n")
pengguna19758
sumber
1
Tidak juga - Anda telah mengasumsikan distribusi normal dan interval yang sama tentang rata-rata, yang keduanya tidak dapat dibenarkan dalam konteks ini. Bahkan Anda telah menangkap sekitar dari probabilitas, dengan asumsi ini adalah distribusi diskrit, dan perlu sedikit memperluas interval Anda untuk mendapatkan . Mungkin lebih baik mengambil daerah dengan kepadatan tertinggi jika ini merupakan distribusi diskrit. Sebagai alternatif ambil interval sehingga probabilitas berada di bawahnya adalah atau kurang, dan probabilitas berada di atas itu adalah atau kurang, juga sini. 95 % [ 4.4 , 8.0 ] 2.5 % 2.5 % [ 4.4 , 8.0 ]94%95%[4.4,8.0]2.5%2.5%[4.4,8.0]
Henry

Jawaban:

24

Seperti dicatat oleh Henry , Anda mengasumsikan distribusi normal dan itu sangat oke jika data Anda mengikuti distribusi normal, tetapi akan salah jika Anda tidak dapat mengasumsikan distribusi normal untuk itu. Di bawah ini saya jelaskan dua pendekatan berbeda yang dapat Anda gunakan untuk distribusi yang tidak diketahui dengan hanya memberikan titik data xdan perkiraan kepadatan yang menyertainya px.

Hal pertama yang perlu dipertimbangkan adalah apa yang ingin Anda rangkum menggunakan interval Anda. Misalnya, Anda mungkin tertarik pada interval yang diperoleh dengan menggunakan kuantil, tetapi Anda juga bisa tertarik pada wilayah kepadatan tertinggi (lihat di sini , atau di sini ) dari distribusi Anda. Walaupun ini seharusnya tidak membuat banyak (jika ada) perbedaan dalam kasus-kasus sederhana seperti distribusi simetris, unimodal, ini akan membuat perbedaan untuk distribusi yang lebih "rumit". Secara umum, kuantil akan memberikan Anda interval berisi massa probabilitas yang terkonsentrasi di sekitar median ( dari distribusi Anda), sedangkan wilayah kepadatan tertinggi adalah wilayah di sekitar mode100α%dari distribusi. Ini akan lebih jelas jika Anda membandingkan dua plot pada gambar di bawah ini - kuantil "memotong" distribusi secara vertikal, sedangkan wilayah dengan kepadatan tertinggi "memotong" secara horizontal.

Interval kuantitas vs HDR

Hal berikutnya yang perlu dipertimbangkan adalah bagaimana menangani fakta bahwa Anda memiliki informasi yang tidak lengkap tentang distribusi (dengan asumsi bahwa kita berbicara tentang distribusi berkelanjutan, Anda hanya memiliki banyak poin daripada fungsi). Apa yang dapat Anda lakukan adalah mengambil nilai "apa adanya", atau menggunakan semacam interpolasi, atau menghaluskan, untuk mendapatkan nilai "di antara".

Salah satu pendekatan akan menggunakan interpolasi linier (lihat ?approxfundi R), atau alternatifnya sesuatu yang lebih halus seperti splines (lihat ?splinefundi R). Jika Anda memilih pendekatan seperti itu, Anda harus ingat bahwa algoritma interpolasi tidak memiliki pengetahuan domain tentang data Anda dan dapat mengembalikan hasil yang tidak valid seperti nilai di bawah nol dll.

# grid of points
xx <- seq(min(x), max(x), by = 0.001)

# interpolate function from the sample
fx <- splinefun(x, px) # interpolating function
pxx <- pmax(0, fx(xx)) # normalize so prob >0

Pendekatan kedua yang dapat Anda pertimbangkan adalah menggunakan distribusi kepadatan / campuran kernel untuk memperkirakan distribusi Anda menggunakan data yang Anda miliki. Bagian rumit di sini adalah untuk memutuskan bandwidth yang optimal.

# density of kernel density/mixture distribution
dmix <- function(x, m, s, w) {
  k <- length(m)
  rowSums(vapply(1:k, function(j) w[j]*dnorm(x, m[j], s[j]), numeric(length(x))))
}

# approximate function using kernel density/mixture distribution
pxx <- dmix(xx, x, rep(0.4, length.out = length(x)), px) # bandwidth 0.4 chosen arbitrary

Selanjutnya, Anda akan menemukan interval minat. Anda dapat melanjutkan secara numerik, atau dengan simulasi.

1a) Pengambilan sampel untuk mendapatkan interval kuantil

# sample from the "empirical" distribution
samp <- sample(xx, 1e5, replace = TRUE, prob = pxx)

# or sample from kernel density
idx <- sample.int(length(x), 1e5, replace = TRUE, prob = px)
samp <- rnorm(1e5, x[idx], 0.4) # this is arbitrary sd

# and take sample quantiles
quantile(samp, c(0.05, 0.975)) 

1b) Pengambilan sampel untuk mendapatkan wilayah dengan kepadatan tertinggi

samp <- sample(pxx, 1e5, replace = TRUE, prob = pxx) # sample probabilities
crit <- quantile(samp, 0.05) # boundary for the lower 5% of probability mass

# values from the 95% highest density region
xx[pxx >= crit]

2a) Temukan kuantil secara numerik

cpxx <- cumsum(pxx) / sum(pxx)
xx[which(cpxx >= 0.025)[1]]   # lower boundary
xx[which(cpxx >= 0.975)[1]-1] # upper boundary

2b) Temukan daerah kepadatan tertinggi secara numerik

const <- sum(pxx)
spxx <- sort(pxx, decreasing = TRUE) / const
crit <- spxx[which(cumsum(spxx) >= 0.95)[1]] * const

Seperti yang dapat Anda lihat pada plot di bawah ini, jika unimodal, distribusi simetris kedua metode mengembalikan interval yang sama.

Dua macam interval

Tentu saja, Anda juga dapat mencoba menemukan interval sekitar beberapa nilai sentral sehingga dan menggunakan beberapa jenis optimasi untuk menemukan sesuai , tetapi dua pendekatan yang dijelaskan di atas tampaknya digunakan lebih umum dan lebih intuitif.100α%ζPr(Xμ±ζ)αζ

Tim
sumber
Mengapa Anda mengambil sampel ketika Anda bisa menghitung kuantil langsung dari informasi yang diberikan (menggunakan salah satu metode)?
whuber
1
@whuber karena murah dan mudah, tapi saya akan mengedit untuk menjelaskan perhitungan non-simulasi besok.
Tim
Hai Tim, ini sangat membantu. Bukankah tidak benar juga mengambil kuantil dari distibrution. (lebih rendah <- x [yang (as.logical (diff (cumsum (px) / sum (px)> 0,025)))]) (lebih tinggi <- x [yang (as.logical (diff (cumsum (px) / sum)) (px) <0,975)))])
user19758
@ user19758 silakan periksa edit saya.
Tim
+1 Penjelasan tambahan, ilustrasi, dan kode menetapkan standar tinggi untuk jawaban di situs ini. Terima kasih!
whuber