Hitung interval kepercayaan untuk rata-rata distribusi beta

12

Pertimbangkan distribusi beta untuk satu set peringkat di [0,1]. Setelah menghitung rata-rata:

μ=αα+β

Apakah ada cara untuk memberikan interval kepercayaan sekitar rata-rata ini?

Dominikus
sumber
1
dominic - Anda telah menentukan rata- rata populasi . Interval kepercayaan akan didasarkan pada beberapa estimasi rata-rata itu. Statistik sampel apa yang Anda gunakan?
Glen_b -Reinstate Monica
Glen_b - Hai, saya menggunakan seperangkat peringkat yang dinormalisasi (dari suatu produk) dalam interval [0,1]. Apa yang saya cari adalah perkiraan interval di sekitar rata-rata (untuk tingkat kepercayaan tertentu), misalnya: mean + - 0,02
dominic
2
dominic: Biarkan saya coba lagi. Anda tidak tahu artinya populasi . Jika Anda ingin perkiraan berada di tengah-tengah interval Anda ( estimasi setengah-lebar , seperti dalam komentar Anda), Anda akan memerlukan beberapa penaksir untuk jumlah itu di urutan tengah untuk menempatkan interval di sekitarnya. Apa yang kamu gunakan untuk itu? Kemungkinan maksimum? Metode saat? sesuatu yang lain? ±
Glen_b -Reinstate Monica
Glen_b - terima kasih atas kesabaran Anda. Saya akan menggunakan MLE
dominic
2
Dominikus; dalam hal itu, untuk yang besar akan menggunakan sifat asimptotik dari penduga kemungkinan maksimum; estimasi ML dari akan didistribusikan secara asimptotik secara normal dengan mean dan kesalahan standar yang dapat dihitung dari Informasi Fisher . Dalam sampel kecil seseorang kadang-kadang dapat menghitung distribusi MLE (meskipun dalam kasus beta saya ingat bahwa sulit); alternatifnya adalah mensimulasikan distribusi pada ukuran sampel Anda untuk memahami perilakunya di sana. μ μnμμ
Glen_b -Reinstate Monica

Jawaban:

22

Meskipun ada metode khusus untuk menghitung interval kepercayaan untuk parameter dalam distribusi beta, saya akan menjelaskan beberapa metode umum, yang dapat digunakan untuk (hampir) semua jenis distribusi , termasuk distribusi beta, dan mudah diimplementasikan dalam R .

Interval kepercayaan kemungkinan profil

Mari kita mulai dengan estimasi kemungkinan maksimum dengan interval kepercayaan kemungkinan profil yang sesuai. Pertama kita perlu beberapa data sampel:

# Sample size
n = 10

# Parameters of the beta distribution
alpha = 10
beta = 1.4

# Simulate some data
set.seed(1)
x = rbeta(n, alpha, beta)

# Note that the distribution is not symmetrical
curve(dbeta(x,alpha,beta))

Fungsi kepadatan probabilitas untuk distribusi beta.

Maksud sebenarnya / teoritis adalah

> alpha/(alpha+beta)
0.877193

Sekarang kita harus membuat fungsi untuk menghitung fungsi kemungkinan log negatif untuk sampel dari distribusi beta, dengan mean sebagai salah satu parameter. Kita dapat menggunakan dbeta()fungsi, tetapi karena ini tidak menggunakan parametrisation yang melibatkan mean, kita harus menyatakan parameternya ( α dan β ) sebagai fungsi dari mean dan beberapa parameter lainnya (seperti standar deviasi):

# Negative log likelihood for the beta distribution
nloglikbeta = function(mu, sig) {
  alpha = mu^2*(1-mu)/sig^2-mu
  beta = alpha*(1/mu-1)
  -sum(dbeta(x, alpha, beta, log=TRUE))
}

Untuk menemukan perkiraan kemungkinan maksimum, kita dapat menggunakan mle()fungsi di stats4perpustakaan:

library(stats4)
est = mle(nloglikbeta, start=list(mu=mean(x), sig=sd(x)))

Abaikan saja peringatan untuk saat ini. Itu disebabkan oleh algoritma pengoptimalan yang mencoba nilai yang tidak valid untuk parameter, memberikan nilai negatif untuk α dan / atau β . (Untuk menghindari peringatan, Anda dapat menambahkan lowerargumen dan mengubah optimasi yang methoddigunakan.)

Sekarang kami memiliki taksiran dan interval kepercayaan untuk dua parameter kami:

> est
Call:
mle(minuslogl = nloglikbeta, start = list(mu = mean(x), sig = sd(x)))

Coefficients:
        mu        sig 
0.87304148 0.07129112

> confint(est)
Profiling...
         2.5 %    97.5 %
mu  0.81336555 0.9120350
sig 0.04679421 0.1276783

Perhatikan bahwa, seperti yang diharapkan, interval kepercayaan tidak simetris:

par(mfrow=c(1,2))
plot(profile(est)) # Profile likelihood plot

Alur kemungkinan profil untuk distribusi beta.

(Garis magenta kedua-luar menunjukkan interval kepercayaan 95%.)

Juga perhatikan bahwa bahkan hanya dengan 10 pengamatan, kami mendapatkan perkiraan yang sangat baik (interval kepercayaan yang sempit).

Sebagai alternatif mle(), Anda dapat menggunakan fitdistr()fungsi dari MASSpaket. Ini juga menghitung penaksir kemungkinan maksimum, dan memiliki keuntungan bahwa Anda hanya perlu menyediakan kerapatan, bukan kemungkinan log negatif, tetapi tidak memberi Anda interval kepercayaan kemungkinan profil, hanya interval kepercayaan asimtotik (simetris).

Pilihan yang lebih baik adalah mle2()(dan fungsi terkait) dari bbmlepaket, yang agak lebih fleksibel dan kuat daripada mle(), dan memberikan plot yang sedikit lebih bagus.

Interval kepercayaan bootstrap

Pilihan lain adalah menggunakan bootstrap. Ini sangat mudah digunakan dalam R, dan Anda bahkan tidak perlu menyediakan fungsi kerapatan:

> library(simpleboot)
> x.boot = one.boot(x, mean, R=10^4)
> hist(x.boot)                # Looks good
> boot.ci(x.boot, type="bca") # Confidence interval
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = x.boot, type = "bca")

Intervals : 
Level       BCa          
95%   ( 0.8246,  0.9132 )  
Calculations and Intervals on Original Scale

Bootstrap memiliki keunggulan tambahan yang berfungsi meskipun data Anda tidak berasal dari distribusi beta.

Interval kepercayaan asimptotik

Untuk interval kepercayaan pada mean, jangan lupakan interval kepercayaan asimptotik lama yang baik berdasarkan teorema limit pusat (dan distribusi- t ). Selama kita memiliki ukuran sampel yang besar (sehingga CLT berlaku dan distribusi rata-rata sampel mendekati normal) atau nilai-nilai besar dari α dan β (sehingga distribusi beta sendiri kira-kira normal), ia bekerja dengan baik. Di sini kita tidak memiliki keduanya, tetapi interval kepercayaan masih tidak terlalu buruk:

> t.test(x)$conf.int
[1] 0.8190565 0.9268349

Untuk nilai n yang sedikit lebih besar (dan nilai yang tidak terlalu ekstrim dari kedua parameter), interval kepercayaan asimptotik bekerja sangat baik.

Karl Ove Hufthammer
sumber
Terima kasih Karl. Pertanyaan cepat: bagaimana Anda menentukan alfa dan beta Anda? Saya menggunakan varians dan mean sampel untuk mendapatkan alfa dan beta, tapi saya pikir saya mungkin telah mengacaukan mean sampel dengan mean populasi, jadi saya tidak yakin saya melakukannya dengan cara yang benar ... lihat komentar Glen_b di atas .
Dominic
Untuk menentukan α dan β sebagai fungsi dari mean dan standar deviasi, saya hanya membalikkan fungsi untuk mean dan standar deviasi sebagai fungsi α dan β (tapi saya yakin Anda juga dapat mencarinya di internet).
Karl Ove Hufthammer
α,β
0

Lihat regresi Beta. Pengantar yang baik untuk melakukannya menggunakan R dapat ditemukan di sini:

http://cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf

Cara lain (sangat mudah) untuk membangun interval kepercayaan adalah dengan menggunakan pendekatan boostrap non-parametrik. Wikipedia memiliki info bagus:

http://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29

Juga video yang bagus di sini:

http://www.youtube.com/watch?v=ZCXg64l9R_4

Rasmus Bååth
sumber