Seperti apa model Bayesian yang kuat untuk memperkirakan skala distribusi normal?

32

Ada sejumlah penaksir skala yang kuat . Contoh penting adalah deviasi absolut median yang berhubungan dengan deviasi standar sebagai σ=MAD1.4826 . Dalam kerangka kerja Bayesian ada sejumlah cara untuk memperkirakan dengan kuat lokasi distribusi yang kira-kira normal (misalnya, Normal yang terkontaminasi oleh pencilan), misalnya, dapat diasumsikan bahwa data didistribusikan sebagai distribusi atau distribusi Laplace. Sekarang pertanyaan saya:

Akan seperti apa model Bayesian untuk mengukur skala distribusi normal dengan cara yang kuat, kuat dalam arti yang sama dengan MAD atau penduga kuat yang serupa?

Seperti halnya dengan MAD, akan lebih rapi jika model Bayesian dapat mendekati SD dari distribusi normal dalam kasus ketika distribusi data sebenarnya terdistribusi normal.

edit 1:

Sebuah contoh khas dari model yang kuat terhadap kontaminasi / outlier ketika asumsi data yi adalah kira-kira normal menggunakan di distribusi seperti:

yit(m,s,ν)

Di mana m adalah mean, s adalah skalanya, dan ν adalah derajat kebebasan. Dengan prior yang sesuai pada m,s dan ν , m akan menjadi perkiraan rata-rata yi yang akan kuat terhadap outlier. Namun, s tidak akan menjadi estimasi yang konsisten dari SD yi karena s tergantung pada ν . Misalnya, jika ν akan diperbaiki ke 4.0 dan model di atas akan dipasang ke sejumlah besar sampel dari Norm(μ=0,σ=1) distribusi makas akan menjadi sekitar 0,82. Yang saya cari adalah model yang kuat, seperti model t, tetapi untuk SD alih-alih (atau sebagai tambahan) artinya.

edit 2:

Berikut mengikuti contoh kode dalam R dan JAGS tentang bagaimana model-t yang disebutkan di atas lebih kuat sehubungan dengan mean.

# generating some contaminated data
y <- c( rnorm(100, mean=10, sd=10), 
        rnorm(10, mean=100, sd= 100))

#### A "standard" normal model ####
model_string <- "model{
  for(i in 1:length(y)) {
    y[i] ~ dnorm(mu, inv_sigma2)
  }

  mu ~ dnorm(0, 0.00001)
  inv_sigma2 ~ dgamma(0.0001, 0.0001)
  sigma <- 1 / sqrt(inv_sigma2)
}"

model <- jags.model(textConnection(model_string), list(y = y))
mcmc_samples <- coda.samples(model, "mu", n.iter=10000)
summary(mcmc_samples)

### The quantiles of the posterior of mu
##  2.5%   25%   50%   75% 97.5% 
##   9.8  14.3  16.8  19.2  24.1 

#### A (more) robust t-model ####
library(rjags)
model_string <- "model{
  for(i in 1:length(y)) {
    y[i] ~ dt(mu, inv_s2, nu)
  }

  mu ~ dnorm(0, 0.00001)
  inv_s2 ~ dgamma(0.0001,0.0001)
  s <- 1 / sqrt(inv_s2)
  nu ~ dexp(1/30) 
}"

model <- jags.model(textConnection(model_string), list(y = y))
mcmc_samples <- coda.samples(model, "mu", n.iter=1000)
summary(mcmc_samples)

### The quantiles of the posterior of mu
## 2.5%   25%   50%   75% 97.5% 
##8.03  9.35  9.99 10.71 12.14 
Rasmus Bååth
sumber
Mungkin itu tidak cukup kuat, tetapi distribusi chi-squared adalah konjugat yang biasanya dipilih sebelum untuk kebalikan dari varians.
Mike Dunlavey
Anda mungkin ingin melihat apakah jawaban pertama untuk pertanyaan ini stats.stackexchange.com/questions/6493/… sudah cukup untuk Anda; mungkin tidak, tapi mungkin memang begitu.
jbowman
Apa yang Anda lakukan sebelumnya untuk tingkat kontaminasi? Apakah kontaminasi akan sistematis? Acak? Apakah akan dihasilkan oleh satu distribusi, atau beberapa distribusi? Apakah kita tahu sesuatu tentang distribusi kebisingan? Jika setidaknya beberapa hal di atas diketahui, maka kita dapat memuat semacam model campuran. Kalau tidak, saya tidak yakin apa keyakinan Anda tentang masalah ini sebenarnya, dan jika Anda tidak memiliki apa pun dari ini sepertinya pengaturan yang sangat kabur. Anda perlu memperbaiki sesuatu, jika tidak, Anda dapat secara acak memilih suatu titik dan menyatakannya sebagai satu-satunya titik yang dihasilkan dengan Gaussian.
Berarti-makna
Tetapi secara umum, Anda bisa memasukkan distribusi-t yang lebih tahan terhadap pencilan, atau campuran distribusi-t. Saya yakin ada banyak makalah, di sini ada satu oleh Bishop research.microsoft.com/en-us/um/people/cmbishop/downloads/… dan di sini ada paket-R yang sesuai dengan campuran: maths.uq.edu. au / ~ gjm / mix_soft / EMMIX_R / EMMIX-manual.pdf
means-to-meaning
1
σ=MAD1.4826

Jawaban:

10

Bayesian inference dalam model T noise dengan prior yang tepat akan memberikan estimasi lokasi dan skala yang kuat. Kondisi tepat yang kemungkinan dan kebutuhan sebelumnya untuk dipenuhi diberikan dalam makalah Bayesian modelling parameter lokasi dan skala oleh Andrade dan O'Hagan (2011). Estimasi ini kuat dalam arti bahwa pengamatan tunggal tidak dapat membuat estimasi besar secara sewenang-wenang, seperti yang ditunjukkan pada gambar 2 dari makalah ini.

νσsss=σf(ν)f

library(stats)
library(stats4)
y = rnorm(100000, mean=0,sd=1)
nu = 4
nLL = function(s) -sum(stats::dt(y/s,nu,log=TRUE)-log(s))
fit = mle(nLL, start=list(s=1), method="Brent", lower=0.5, upper=2)
# the variance of a standard T is nu/(nu-2)
print(coef(fit)*sqrt(nu/(nu-2)))

ν=4f(ν)=1.18σ^=s/f(ν)

Tom Minka
sumber
1
Jawaban yang bagus (+1). 'dalam arti bahwa pengamatan tunggal tidak dapat membuat perkiraan besar secara sewenang-wenang,' jadi titik rinciannya adalah 2 / n (saya bertanya-tanya tentang ini) .... Sebagai titik perbandingan, untuk prosedur yang diilustrasikan dalam jawaban saya itu adalah n / 2.
user603
Wow terima kasih! Pertanyaan tindak lanjut yang kabur. Apakah kemudian benar - benar masuk akal untuk "memperbaiki" skala sehingga konsisten dengan SD dalam kasus Normal? Kasus penggunaan yang saya pikirkan adalah ketika melaporkan ukuran spread. Saya tidak akan memiliki masalah dengan skala pelaporan, tetapi akan menyenangkan untuk melaporkan sesuatu yang akan konsisten dengan SD karena itu adalah ukuran penyebaran yang paling umum (setidaknya dalam psikologi). Apakah Anda melihat situasi di mana koreksi ini akan menyebabkan perkiraan yang aneh dan tidak konsisten?
Rasmus Bååth
6

Saat Anda mengajukan pertanyaan tentang masalah yang sangat tepat (estimasi kuat), saya akan menawarkan jawaban yang sama persis. Namun, pertama-tama, saya akan mulai mencoba menghilangkan asumsi yang tidak beralasan. Tidak benar bahwa ada perkiraan bayesian lokasi yang kuat (ada perkiraan bayesian lokasi tetapi seperti yang saya ilustrasikan di bawah mereka tidak kuat dan, tampaknya , bahkan estimator kuat yang paling sederhana dari lokasi bukanlah bayesian). Menurut pendapat saya, alasan tidak adanya tumpang tindih antara paradigma 'bayesian' dan 'kuat' dalam kasus lokasi sangat membantu dalam menjelaskan mengapa tidak ada penduga penyebaran yang baik dan bayesian.

m,sνmyi

Sebenarnya tidak. Perkiraan yang dihasilkan hanya akan kuat dalam arti kata yang sangat lemah. Namun, ketika kami mengatakan bahwa median kuat untuk outlier kami berarti kata kuat dalam arti yang jauh lebih kuat. Yaitu, dalam statistik yang kuat, kekokohan median mengacu pada properti yang jika Anda menghitung median pada kumpulan data pengamatan yang diambil dari uni-modal, model kontinu dan kemudian mengganti kurang dari setengah pengamatan ini dengan nilai acak. , nilai median yang dihitung pada data yang terkontaminasi dekat dengan nilai yang Anda miliki seandainya Anda menghitungnya pada set data asli (tidak terkontaminasi). Maka, mudah untuk menunjukkan bahwa strategi estimasi yang Anda usulkan dalam paragraf yang saya kutip di atas jelas tidak kuat dalam arti bagaimana kata tersebut biasanya dipahami untuk median.

F

  1. zi=|ximed(x)|mad(x)
  2. zi>qα(z|xF)αzxFF
  3. Jalankan analisis Bayesian (biasanya, tidak kuat) pada pengamatan yang tidak ditolak.

EDIT:

Terima kasih kepada OP untuk menyediakan kode R yang lengkap untuk melakukan analisis bayesian bonna fide masalah.

n/22

N(1000,1)

n<-100
set.seed(123)
y<-rnorm(n,1000,1)

Tambahkan sejumlah kontaminan:

y[1:30]<-y[1:30]/100-1000 
w<-rep(0,n)
w[1:30]<-1

indeks w mengambil nilai 1 untuk outlier. Saya mulai dengan pendekatan yang disarankan oleh OP:

library("rjags")
model_string<-"model{
  for(i in 1:length(y)){
    y[i]~dt(mu,inv_s2,nu)
  }
  mu~dnorm(0,0.00001)
  inv_s2~dgamma(0.0001,0.0001)
  s<-1/sqrt(inv_s2)
  nu~dexp(1/30) 
}"

model<-jags.model(textConnection(model_string),list(y=y))
mcmc_samples<-coda.samples(model,"mu",n.iter=1000)
print(summary(mcmc_samples)$statistics[1:2])
summary(mcmc_samples)

Saya mendapat:

     Mean        SD 
384.2283  97.0445 

dan:

2. Quantiles for each variable:

 2.5%   25%   50%   75% 97.5% 
184.6 324.3 384.7 448.4 577.7 

(Diam jauh dari nilai target)

Untuk metode yang kuat,

z<-abs(y-median(y))/mad(y)
th<-max(abs(rnorm(length(y))))
print(c(mean(y[which(z<=th)]),sd(y[which(z<=th)])))

satu mendapat:

 1000.149 0.8827613

(sangat dekat dengan nilai target)

zthF
t

  • [1] Ricardo A. Maronna, Douglas R. Martin, Victor J. Yohai (2006). Statistik Kuat: Teori dan Metode (Seri Wiley dalam Probabilitas dan Statistik).
  • Huber, PJ (1981). Statistik yang kuat. New York: John Wiley and Sons.
pengguna603
sumber
1
Nah, t sering diusulkan sebagai alternatif yang kuat untuk distribusi normal. Saya tidak tahu apakah ini dalam arti lemah atau tidak. Lihat misalnya: Lange, KL, Little, RJ, & Taylor, JM (1989). Pemodelan statistik yang kuat menggunakan distribusi t. Jurnal Asosiasi Statistik Amerika , 84 (408), 881-896. pdf
Rasmus Bååth
1
Ini adalah indra yang lemah. Jika Anda memiliki kode R yang mengimplementasikan prosedur yang Anda sarankan, saya akan dengan senang hati mengilustrasikan jawaban saya dengan sebuah contoh. jika tidak, Anda bisa mendapatkan penjelasan lebih lanjut di bab 2 buku teks ini .
user603
Prosedur yang saya sarankan pada dasarnya dijelaskan di sini: indiana.edu/~kruschke/BEST termasuk kode R. Saya harus memikirkan solusi Anda! Namun, tampaknya bukan Bayesian dalam arti bahwa ia tidak memodelkan semua data, hanya subset yang "bertahan" langkah 2.
Rasmus Bååth
1
Saya sekarang telah melakukannya!
Rasmus Bååth
1

Dalam analisis bayesian menggunakan distribusi Gamma terbalik sebagai prior untuk presisi (kebalikan dari varians) adalah pilihan umum. Atau distribusi Wishart terbalik untuk model multivarian. Menambahkan sebelumnya pada varians meningkatkan ketahanan terhadap outlier.

Ada sebuah makalah yang bagus dari Andrew Gelman: "Distribusi sebelumnya untuk parameter varians dalam model hierarkis" di mana ia membahas apa pilihan yang baik untuk prior pada varian dapat.

jpmuc
sumber
4
Maaf tapi saya gagal melihat bagaimana ini menjawab pertanyaan. Saya tidak meminta yang kuat sebelumnya, tetapi untuk model yang kuat .
Rasmus Bååth
0

μNσ2μtN

σD

D|μ,σN(μ,σ2)
D(d1,,dN)
p(D|μ,σ2)=1(2πσ)Nexp(N2σ2((mμ2)+s2))
ms2
m=1Ni=1Ndis2=1Ni=1Ndi2m2
p(μ,σ2|D)p(D|μ,σ2)p(μ,σ2)
(μ,σ2)p(μ,σ2|D)p(σ2|D)
σ2|DIG(α+N/2,2β+Ns2)α,β>0
σ2αβtμ
yannick
sumber
1
σ2
1
Itu semua tergantung pada apa yang Anda maksud dengan robust. Apa yang Anda katakan saat ini adalah bahwa Anda ingin ketahanan data. Apa yang saya usulkan adalah kekokohan model wrt spesifikasi. Keduanya merupakan jenis ketahanan yang berbeda.
yannick
2
Saya akan mengatakan bahwa contoh-contoh yang saya berikan, MAD dan gunakan pada distribusi sebagai distribusi untuk data adalah contoh-contoh ketahanan sehubungan dengan data.
Rasmus Bååth
Saya akan mengatakan Rasmus benar dan begitu juga Gelman er al dalam BDA3, seperti halnya pemahaman dasar bahwa distribusi ini memiliki ekor yang lebih gemuk daripada normal untuk parameter lokasi yang sama
Brash Equilibrium
0

Saya telah mengikuti diskusi dari pertanyaan awal. Rasmus ketika Anda mengatakan ketahanan saya yakin maksud Anda dalam data (outlier, bukan miss-spesifikasi distribusi). Saya akan mengambil distribusi data menjadi distribusi Laplace alih-alih t-distribusi, kemudian seperti dalam regresi normal di mana kita memodelkan mean, di sini kita akan memodelkan median (sangat kuat) alias regresi median (kita semua tahu). Biarkan model menjadi:

ϵ ( 0 , σ 2 )Y=βX+ϵ , memiliki laplace .ϵ(0,σ2)

Tentu saja tujuan kami adalah memperkirakan parameter model. Kami berharap prior kami tidak jelas untuk memiliki model yang objektif. Model yang ada memiliki posterior dari bentuk . Memberikan sebelum normal dengan varians yang besar membuat seperti sebelum kabur dan sebelum chis-squared dengan kecil derajat kebebasan untuk meniru sebuah jeffrey ini sebelum (samar-samar sebelum) diberikan kepada keβ σ 2f(β,σ,Y,X)βσ2. Dengan sampler Gibbs apa yang terjadi? normal sebelum + laplace likehood = ???? kami tahu. Juga chi-square sebelum + kemungkinan laplace = ??? kita tidak tahu distribusinya. Untungnya bagi kita ada teorema dalam (Aslan, 2010) yang mengubah kemungkinan laplace menjadi skala campuran dari distribusi normal yang kemudian memungkinkan kita untuk menikmati sifat konjugasi dari prior kita. Saya pikir seluruh proses yang dijelaskan sepenuhnya kuat dalam hal outlier. Dalam pengaturan multivariat, chi-square menjadi distribusi wishart, dan kami menggunakan multivariat laplace dan distribusi normal.

Chamberlain Foncha
sumber
2
Solusi Anda tampaknya berfokus pada estimasi lokasi yang kuat (rata-rata / median). Pertanyaan saya agak tentang estimasi skala dengan properti konsistensi sehubungan dengan mengambil SD ketika distribusi data menghasilkan sebenarnya normal.
Rasmus Bååth
Dengan perkiraan lokasi yang kuat, skala sebagai fungsi lokasi langsung mendapat manfaat dari kekokohan lokasi. Tidak ada cara lain untuk membuat skalanya kuat.
Chamberlain Foncha
Pokoknya saya harus mengatakan bahwa saya sedang menunggu untuk melihat bagaimana masalah ini akan ditangani terutama dengan distribusi normal seperti yang Anda tekankan.
Chamberlain Foncha
0

Misalkan Anda memiliki grup dan Anda ingin memodelkan distribusi varians sampel mereka, mungkin terkait dengan beberapa kovariat . Yaitu, anggap titik data Anda untuk grup adalah . Pertanyaannya di sini adalah, "Apa model yang kuat untuk kemungkinan varians sampel?" Salah satu cara untuk pendekatan ini adalah untuk model data ditransformasikan berasal darix k 1 ... K Var ( y k ) [ 0 , ) ln [ Var ( y k ) ] tKxk1KVar(yk)[0,)ln[Var(yk)]tn, maka Anda dapat memilih distribusi probabilitas dengan dukungan nyata positif yang diketahui memiliki ekor berat dibandingkan dengan distribusi lain dengan lokasi yang sama. Sebagai contoh, ada jawaban baru-baru ini untuk pertanyaan di Cross Divalidasi tentang apakah distribusi lognormal atau gamma memiliki ekor yang lebih berat, dan ternyata distribusi lognormal melakukannya (terima kasih kepada @Glen_b untuk kontribusinya). Selain itu, Anda bisa menjelajahi keluarga setengah-Cauchy.

Alasan yang sama berlaku jika Anda menetapkan distribusi sebelumnya atas parameter skala untuk distribusi normal. Intinya, distribusi lognormal dan invers-gamma tidak disarankan jika Anda ingin membentuk batas menghindari sebelum untuk keperluan perkiraan mode posterior karena mereka memuncak tajam jika Anda men-parameterkannya sehingga mode mendekati nol. Lihat BDA3 bab 13 untuk diskusi. Jadi selain mengidentifikasi model yang kuat dalam hal ketebalan ekor, perlu diingat bahwa kurtosis juga penting bagi kesimpulan Anda.

Saya harap ini membantu Anda sebanyak jawaban Anda untuk salah satu pertanyaan terakhir saya membantu saya.

Keseimbangan kurang ajar
sumber
1
Pertanyaan saya adalah tentang situasi ketika Anda memiliki satu kelompok dan bagaimana memperkirakan dengan kuat skala kelompok itu. Dalam kasus outlier saya tidak percaya varians sampel dianggap kuat.
Rasmus Bååth
Jika Anda memiliki satu grup, dan Anda memperkirakan distribusi normalnya, maka pertanyaan Anda berlaku untuk bentuk parameter sebelumnya melebihi parameter skalanya. Seperti jawaban saya, Anda dapat menggunakan distribusi pada transformasi log-nya atau memilih distribusi tailed fat dengan dukungan nyata positif, berhati-hati tentang aspek-aspek lain dari distribusi itu seperti kurtosis-nya. Intinya, jika Anda menginginkan model yang kuat untuk parameter skala, gunakan pada distribusi melalui log-nya atau distribusi lemak lainnya.
Brash Equilibrium