Kapan interval kepercayaan "masuk akal" tetapi interval kredibel yang sesuai tidak?

14

Sering terjadi bahwa interval kepercayaan dengan cakupan 95% sangat mirip dengan interval kredibel yang mengandung 95% dari kepadatan posterior. Ini terjadi ketika yang sebelumnya seragam atau hampir seragam dalam kasus yang terakhir. Dengan demikian interval kepercayaan sering dapat digunakan untuk memperkirakan interval yang kredibel dan sebaliknya. Yang penting, kita dapat menyimpulkan dari hal ini bahwa kesalahan interpretasi yang keliru dari interval kepercayaan sebagai interval yang kredibel memiliki sedikit atau tidak ada kepentingan praktis untuk banyak kasus penggunaan sederhana.

Ada sejumlah contoh di luar sana kasus di mana ini tidak terjadi, namun mereka semua tampaknya dicintai oleh pendukung statistik Bayesian dalam upaya untuk membuktikan ada sesuatu yang salah dengan pendekatan yang sering terjadi. Dalam contoh-contoh ini, kita melihat interval kepercayaan berisi nilai-nilai yang tidak mungkin, dll yang seharusnya menunjukkan bahwa mereka tidak masuk akal.

Saya tidak ingin membahas kembali contoh-contoh itu, atau diskusi filosofis tentang Bayesian vs Frequentist.

Saya hanya mencari contoh yang sebaliknya. Apakah ada kasus di mana interval kepercayaan dan kredibilitas berbeda secara substansial, dan interval yang disediakan oleh prosedur kepercayaan jelas lebih unggul?

Untuk memperjelas: Ini tentang situasi ketika interval yang kredibel biasanya diharapkan bertepatan dengan interval kepercayaan yang sesuai, yaitu ketika menggunakan flat, seragam, dll. Saya tidak tertarik pada kasus di mana seseorang memilih yang sebelumnya buruk.

EDIT: Menanggapi jawaban @JaeHyeok Shin di bawah, saya harus tidak setuju bahwa contohnya menggunakan kemungkinan yang benar. Saya menggunakan perhitungan bayesian perkiraan untuk memperkirakan distribusi posterior yang benar untuk theta di bawah ini di R:

### Methods ###
# Packages
require(HDInterval)

# Define the likelihood
like <- function(k = 1.2, theta = 0, n_print = 1e5){
  x    = NULL
  rule = FALSE
  while(!rule){
    x     = c(x, rnorm(1, theta, 1))
    n     = length(x)
    x_bar = mean(x)

    rule = sqrt(n)*abs(x_bar) > k

    if(n %% n_print == 0){ print(c(n, sqrt(n)*abs(x_bar))) }
  }
  return(x)
}

# Plot results
plot_res <- function(chain, i){
    par(mfrow = c(2, 1))
    plot(chain[1:i, 1], type = "l", ylab = "Theta", panel.first = grid())
    hist(chain[1:i, 1], breaks = 20, col = "Grey", main = "", xlab = "Theta")
}


### Generate target data ### 
set.seed(0123)
X = like(theta = 0)
m = mean(X)


### Get posterior estimate of theta via ABC ###
tol   = list(m = 1)
nBurn = 1e3
nStep = 1e4


# Initialize MCMC chain
chain           = as.data.frame(matrix(nrow = nStep, ncol = 2))
colnames(chain) = c("theta", "mean")
chain$theta[1]  = rnorm(1, 0, 10)

# Run ABC
for(i in 2:nStep){
  theta = rnorm(1, chain[i - 1, 1], 10)
  prop  = like(theta = theta)

  m_prop = mean(prop)


  if(abs(m_prop - m) < tol$m){
    chain[i,] = c(theta, m_prop)
  }else{
    chain[i, ] = chain[i - 1, ]
  }
  if(i %% 100 == 0){ 
    print(paste0(i, "/", nStep)) 
    plot_res(chain, i)
  }
}

# Remove burn-in
chain = chain[-(1:nBurn), ]

# Results
plot_res(chain, nrow(chain))
as.numeric(hdi(chain[, 1], credMass = 0.95))

Ini adalah interval kredibel 95%:

> as.numeric(hdi(chain[, 1], credMass = 0.95))
[1] -1.400304  1.527371

masukkan deskripsi gambar di sini

EDIT # 2:

Ini adalah pembaruan setelah komentar @JaeHyeok Shin. Saya mencoba untuk membuatnya sesederhana mungkin tetapi skrip menjadi sedikit lebih rumit. Perubahan utama:

  1. Sekarang menggunakan toleransi 0,001 untuk rata-rata (itu 1)
  2. Peningkatan jumlah langkah hingga 500 ribu untuk memperhitungkan toleransi yang lebih kecil
  3. Mengurangi sd distribusi proposal menjadi 1 untuk memperhitungkan toleransi yang lebih kecil (10)
  4. Menambahkan kemungkinan rnorm sederhana dengan n = 2k untuk perbandingan
  5. Menambahkan ukuran sampel (n) sebagai statistik ringkasan, atur toleransi ke 0,5 * n_target

Ini kodenya:

### Methods ###
# Packages
require(HDInterval)

# Define the likelihood
like <- function(k = 1.3, theta = 0, n_print = 1e5, n_max = Inf){
  x    = NULL
  rule = FALSE
  while(!rule){
    x     = c(x, rnorm(1, theta, 1))
    n     = length(x)
    x_bar = mean(x)
    rule  = sqrt(n)*abs(x_bar) > k
    if(!rule){
     rule = ifelse(n > n_max, TRUE, FALSE)
    }

    if(n %% n_print == 0){ print(c(n, sqrt(n)*abs(x_bar))) }
  }
  return(x)
}


# Define the likelihood 2
like2 <- function(theta = 0, n){
  x = rnorm(n, theta, 1)
  return(x)
}



# Plot results
plot_res <- function(chain, chain2, i, main = ""){
    par(mfrow = c(2, 2))
    plot(chain[1:i, 1],  type = "l", ylab = "Theta", main = "Chain 1", panel.first = grid())
    hist(chain[1:i, 1],  breaks = 20, col = "Grey", main = main, xlab = "Theta")
    plot(chain2[1:i, 1], type = "l", ylab = "Theta", main = "Chain 2", panel.first = grid())
    hist(chain2[1:i, 1], breaks = 20, col = "Grey", main = main, xlab = "Theta")
}


### Generate target data ### 
set.seed(01234)
X    = like(theta = 0, n_print = 1e5, n_max = 1e15)
m    = mean(X)
n    = length(X)
main = c(paste0("target mean = ", round(m, 3)), paste0("target n = ", n))



### Get posterior estimate of theta via ABC ###
tol   = list(m = .001, n = .5*n)
nBurn = 1e3
nStep = 5e5

# Initialize MCMC chain
chain           = chain2 = as.data.frame(matrix(nrow = nStep, ncol = 2))
colnames(chain) = colnames(chain2) = c("theta", "mean")
chain$theta[1]  = chain2$theta[1]  = rnorm(1, 0, 1)

# Run ABC
for(i in 2:nStep){
  # Chain 1
  theta1 = rnorm(1, chain[i - 1, 1], 1)
  prop   = like(theta = theta1, n_max = n*(1 + tol$n))
  m_prop = mean(prop)
  n_prop = length(prop)
  if(abs(m_prop - m) < tol$m &&
     abs(n_prop - n) < tol$n){
    chain[i,] = c(theta1, m_prop)
  }else{
    chain[i, ] = chain[i - 1, ]
  }

  # Chain 2
  theta2  = rnorm(1, chain2[i - 1, 1], 1)
  prop2   = like2(theta = theta2, n = 2000)
  m_prop2 = mean(prop2)
  if(abs(m_prop2 - m) < tol$m){
    chain2[i,] = c(theta2, m_prop2)
  }else{
    chain2[i, ] = chain2[i - 1, ]
  }

  if(i %% 1e3 == 0){ 
    print(paste0(i, "/", nStep)) 
    plot_res(chain, chain2, i, main = main)
  }
}

# Remove burn-in
nBurn  = max(which(is.na(chain$mean) | is.na(chain2$mean)))
chain  = chain[ -(1:nBurn), ]
chain2 = chain2[-(1:nBurn), ]


# Results
plot_res(chain, chain2, nrow(chain), main = main)
hdi1 = as.numeric(hdi(chain[, 1],  credMass = 0.95))
hdi2 = as.numeric(hdi(chain2[, 1], credMass = 0.95))


2*1.96/sqrt(2e3)
diff(hdi1)
diff(hdi2)

Hasilnya, di mana hdi1 adalah "kemungkinan" saya dan hdi2 adalah rnorm sederhana (n, theta, 1):

> 2*1.96/sqrt(2e3)
[1] 0.08765386
> diff(hdi1)
[1] 1.087125
> diff(hdi2)
[1] 0.07499163

Jadi setelah menurunkan toleransi secukupnya, dan dengan mengorbankan banyak lagi langkah MCMC, kita dapat melihat lebar CRI yang diharapkan untuk model rnorm.

masukkan deskripsi gambar di sini

Marah
sumber
Tidak duplikat, tetapi memiliki hubungan dekat dengan stats.stackexchange.com/questions/419916/…
user158565
6
Secara umum, ketika Anda memiliki informasi sebelumnya yang sangat salah, dalam arti informal, misalnya, Normal (0,1) ketika nilai aktualnya adalah -3,6, interval kredibel Anda tanpa adanya banyak data akan sangat buruk ketika memandang dari sudut pandang sering.
jbowman
@jbowman Ini khusus tentang kasus ketika menggunakan seragam sebelumnya atau sesuatu seperti N (0, 1e6).
Livid
Beberapa dekade yang lalu, Bayesian yang sebenarnya bernama ahli statistik yang menggunakan non-informatif sebelumnya sebagai pseudo- (atau palsu-) Bayesian.
user158565
@ user158565 Ini offtopic tapi seragam sebelumnya hanya perkiraan. Jika p (H_0) = p (H_1) = p (H_2) = ... = p (H_n) maka semua prior dapat keluar dari aturan Bayes yang membuat perhitungan menjadi lebih mudah. Ini tidak lebih salah daripada menjatuhkan istilah-istilah kecil dari penyebut ketika itu masuk akal.
Livid

Jawaban:

6

Dalam desain eksperimental berurutan, interval kredibel dapat menyesatkan.

(Penafian: Saya tidak berpendapat bahwa ini tidak masuk akal - sangat masuk akal dalam pemikiran Bayesian dan tidak menyesatkan dalam perspektif sudut pandang Bayesian.)

Sebagai contoh sederhana, katakanlah kita memiliki mesin yang memberi kita sampel acak dari dengan tidak diketahui . Alih-alih menggambar iid sampel, kami mengambil sampel sampai untuk tetap . Artinya, jumlah sampel adalah waktu berhenti ditentukan oleh XN(θ,1)θnnX¯n>kkN

N=inf{n1:nX¯n>k}.

Dari hukum iterasi logaritma, kita tahu untuk setiap . Jenis penghentian ini biasanya digunakan dalam pengujian / estimasi berurutan untuk mengurangi jumlah sampel yang harus dilakukan inferensi.Pθ(N<)=1θR

Prinsip kemungkinan menunjukkan bahwa posterior tidak dipengaruhi oleh aturan pemberhentian dan dengan demikian untuk setiap kelancaran wajar sebelumnya (misalnya, , jika kita menetapkan cukup besar , posterior kira-kira dan dengan demikian interval kredibel kira-kira diberikan sebagai Namun, dari definisi , kita tahu bahwa interval yang kredibel ini tidak mengandung jika adalah besar θπ(θ)θN(0,10))kθN(X¯N,1/N)

CIbayes:=[X¯N1.96N,X¯N+1.96N].
N0k
0<X¯N-kNX¯N-1.96N
untuk . Oleh karena itu, cakupan frekuensi adalah nol sejak dan diperoleh ketika adalah . Sebaliknya, cakupan Bayesian kira-kira selalu sama dengan karena k0CsayabSebuahyes
infθPθ(θCsayabSebuahyes)=0,
0θ00,95
P(θCsayabSebuahyes|X1,...,XN)0,95.

Pesan rumah: Jika Anda tertarik untuk memiliki jaminan yang sering, Anda harus berhati-hati menggunakan alat inferensi Bayesian yang selalu berlaku untuk jaminan Bayesian tetapi tidak selalu untuk yang sering.

(Saya belajar contoh ini dari kuliah Larry yang luar biasa. Catatan ini berisi banyak diskusi menarik tentang perbedaan halus antara kerangka kerja frequentist dan Bayesian. Http://www.stat.cmu.edu/~larry/=stat705/Lecture14.pdf )

EDIT Dalam Livid's ABC, nilai toleransi terlalu besar, sehingga bahkan untuk pengaturan standar di mana kami mencicipi sejumlah pengamatan tetap, itu tidak memberikan CR yang benar. Saya tidak terbiasa dengan ABC tetapi jika saya hanya mengubah nilai tol menjadi 0,05, kita dapat memiliki CR yang sangat miring sebagai berikut

> X = like(theta = 0)
> m = mean(X)
> print(m)
[1] 0.02779672

masukkan deskripsi gambar di sini

> as.numeric(hdi(chain[, 1], credMass = 0.95)) [1] -0.01711265 0.14253673

Tentu saja, rantai tidak distabilkan dengan baik tetapi bahkan jika kita menambah panjang rantai, kita bisa mendapatkan CR serupa - condong ke bagian positif.

Bahkan, saya pikir aturan penolakan berdasarkan perbedaan rata-rata tidak cocok dalam pengaturan ini karena dengan probabilitas tinggi dekat dengan jika dan dekat dengan if .NX¯Nk0<θk-k-kθ<0

JaeHyeok Shin
sumber
"Jika kita menetapkan k yang cukup besar, posterior θ kira-kira N (X_N, 1 / N)" . Sepertinya saya yang jelas Pr (X | theta)! = Normal (theta, 1). Yaitu, itu adalah kemungkinan yang salah untuk proses yang menghasilkan urutan Anda. Juga, ada salah ketik. Dalam contoh asli Anda berhenti pengambilan sampel ketika sqrt (n) * abs (rata-rata (x))> k.
Livid
Terima kasih atas komentarnya. Bahkan di bawah aturan stopping, kemungkinan diberikan sebagai . Jadi itu sama dengan produk dari pengamatan N normal meskipun N sekarang acak. Contoh ini masih berlaku untuk aturan penghentian saat ini meskipun yang Anda tunjukkan adalah contoh asli dan historis - Ada sejarah yang menarik tentang seberapa sering dan Bayesian berdebat dengan aturan penghentian ini. Anda dapat melihat, misalnya, errorstatistics.com/2013/04/06/...saya=1Nϕ(Xsaya-θ)
JaeHyeok Shin
Silakan lihat edit saya di pertanyaan. Saya masih berpikir interval kredibel Anda tidak masuk akal karena menggunakan kemungkinan yang salah. Saat menggunakan kemungkinan yang benar seperti dalam kode saya, kami mendapatkan interval yang masuk akal.
Livid
Terima kasih atas percobaan terperinci. Dalam pengaturan Anda, terlalu kecil untuk memenuhi ketidaksetaraan . Perhatikan bahwa harus lebih besar dari 1,96, dan untuk menambah kesalahan perkiraan, saya pikir akan menjadi pilihan yang aman. k0<X¯N-k/NX¯N-1.96/Nkk>10
JaeHyeok Shin
Juga, dapatkah Anda memeriksa kembali apakah komputasi ABC ini berfungsi untuk case standar? Saya mengganti fungsi \ code {like} untuk membuatnya menarik sejumlah pengamatan tetap (n = 2000) Kemudian, panjang teoritis CR harus sekitar tapi saya selalu menjadi sangat CR lebih luas (panjangnya sekitar 2). Saya pikir tingkat toleransi terlalu besar. 2×1.96/2000=0,0876
JaeHyeok Shin
4

Karena interval kredibel dibentuk dari distribusi posterior, berdasarkan pada distribusi sebelumnya yang ditentukan, Anda dapat dengan mudah membangun interval kredibel yang sangat buruk dengan menggunakan distribusi sebelumnya yang sangat terkonsentrasi pada nilai parameter yang sangat tidak masuk akal. Anda dapat membuat interval kredibel yang tidak "masuk akal" dengan menggunakan distribusi sebelumnya yang sepenuhnya terkonsentrasi pada nilai parameter yang tidak mungkin .

Pasang kembali Monica
sumber
1
Atau lebih baik lagi, kredibel yang dibangun oleh prior yang tidak setuju dengan prior Anda (meskipun itu milik orang lain) memiliki peluang bagus untuk menjadi tidak masuk akal bagi Anda. Ini tidak biasa dalam sains; Saya sudah punya peneliti mengatakan mereka tidak ingin memasukkan pendapat ahli, karena dalam pengamatan mereka, para ahli selalu sangat percaya diri.
Cliff AB
1
Ini secara khusus tentang seragam, atau "datar", prior.
Livid
1
@Livid: Anda harus memasukkan bahwa Anda sedang berbicara tentang prior flat dalam pertanyaan Anda. Itu benar-benar mengubah segalanya.
Cliff AB
1
@CliffAB Dalam dua kalimat pertama, tapi saya akan menjelaskan, terima kasih.
Livid
1

Jika kami menggunakan flat sebelumnya, ini hanya permainan di mana kami mencoba membuat flat sebelumnya dalam reparameterisasi yang tidak masuk akal.

Sebagai contoh, misalkan kita ingin membuat kesimpulan tentang suatu probabilitas. Jika kita meletakkan flat sebelum peluang log probabilitas, interval kredibel 95% untuk probabilitas sebenarnya adalah dua poin bahkan sebelum kita melihat data! Jika kita mendapatkan satu titik data positif dan membangun interval yang kredibel 95%, sekarang menjadi titik tunggal .{0,1} {1}

Inilah sebabnya mengapa banyak orang Bayesian keberatan dengan tanah datar.

Cliff AB
sumber
Saya menjelaskan motivasi saya dengan cukup jelas. Saya menginginkan sesuatu seperti contoh di mana interval kepercayaan mencakup nilai yang tidak mungkin, tetapi di mana interval yang kredibel berperilaku baik. Jika contoh Anda bergantung pada melakukan sesuatu yang tidak masuk akal, seperti misalnya memilih kemungkinan yang salah, lalu mengapa itu menarik bagi siapa pun?
Livid
1
@Livid: fungsi kemungkinan sangat masuk akal. The datar sebelum di log-kemungkinan tidak. Dan ini adalah keseluruhan argumen yang digunakan Bayesian untuk mengatakan Anda tidak harus menggunakan prior flat; mereka sebenarnya bisa sangat informatif dan seringkali tidak sesuai keinginan pengguna!
Cliff AB
1
Inilah Andrew Gelman yang membahas beberapa masalah soal harga flat.
Cliff AB
"Flat sebelum log-odds tidak." Maksud saya menempatkan flat di atas peluang log-transformed tampaknya tidak masuk akal bagi Anda, seperti menggunakan kemungkinan yang salah. Maaf, tapi saya tidak terbiasa dengan contoh ini. Apa yang seharusnya dilakukan model ini?
Livid
@Livid: ini mungkin tampak tidak biasa, tetapi sebenarnya tidak! Sebagai contoh, regresi logistik biasanya mempertimbangkan semua parameter pada skala log-odds. Jika Anda memiliki variabel dummy untuk semua grup Anda dan menggunakan prior flat pada parameter regresi Anda, Anda akan mengalami masalah ini persis .
Cliff AB