Interval prediksi untuk variabel acak binomial

14

Apa rumus (perkiraan atau tepat) untuk interval prediksi untuk variabel acak Binomial?

Asumsikan , dan kami mengamati (diambil dari ). The dikenal.YBinom(n,p)yYn

Tujuan kami adalah untuk mendapatkan interval prediksi 95% untuk hasil imbang baru dari .Y

Estimasi titik adalah , di mana . Sebuah selang kepercayaan untuk \ hat {p} sangat mudah, tapi saya tidak dapat menemukan formula untuk interval prediksi untuk Y . Jika kita tahu p (daripada \ hat {p} ), maka interval prediksi 95% hanya melibatkan menemukan kuantil binomial. Apakah ada sesuatu yang jelas saya abaikan?np^p Yp pp^=ynp^Ypp^

Statseeker
sumber
1
Lihat Apa metode non-Bayesian yang ada untuk inferensi prediktif? . Dalam hal ini metode menggunakan pivot tidak tersedia (saya tidak berpikir) tetapi Anda bisa menggunakan salah satu kemungkinan prediksi. Atau tentu saja, pendekatan Bayesian.
Scortchi
1
Hai teman-teman, saya ingin mengambil waktu sejenak untuk mengatasi masalah yang muncul. - Mengenai kepercayaan untuk p: Saya tidak tertarik untuk itu. - tentang prediksi menjadi 95% dari distribusi: ya, itulah interval prediksi apa pun konteksnya (dalam regresi Anda harus mengasumsikan kesalahan normal, sedangkan interval kepercayaan bergantung pada CLT - ya, contoh prediksi jumlah head di flip koin sudah benar. Apa yang membuat masalah ini sulit adalah bahwa kita tidak sekarang "p", tetapi memiliki perkiraan.
Statseeker
3
@Adison membaca buku Interval Statistik oleh G. Hahn dan W. Meeker. Mereka menjelaskan perbedaan antara interval kepercayaan, interval prediksi, interval toleransi dan interval kredibel Bayesian. Interval prediksi 95% tidak mengandung 95% dari distribusi. Itu melakukan apa yang interval paling sering lakukan. Jika Anda berulang kali mengambil sampel dari B (n, p) dan menggunakan metode yang sama setiap kali untuk menghasilkan interval prediksi 95% untuk p maka 95% interval prediksi Anda akan berisi nilai sebenarnya dari p. Jika Anda ingin mencakup 95% dari distribusi, membangun interval toleransi.
Michael R. Chernick
Interval toleransi mencakup persentase dari distribusi. Untuk interval toleransi 95% untuk 90% dari distribusi Anda mengulangi proses berkali-kali dan menggunakan metode yang sama untuk menghasilkan interval setiap kali kemudian dalam sekitar 95% dari kasus setidaknya 90% dari distribusi akan jatuh dalam interval. dan 5% dari waktu kurang dari 90% dari distribusi akan terkandung dalam interval.
Michael R. Chernick
3
Lawless & Fredette (2005), "Interval Prediksi Frekuensi dan Distribusi Prediktif", Biometrika , 92 , 3 adalah referensi bagus lainnya, selain yang ada di tautan yang saya berikan.
Scortchi

Jawaban:

24

Ok, ayo coba ini. Saya akan memberikan dua jawaban - yang Bayesian, yang menurut saya sederhana dan alami, dan salah satu yang mungkin sering.

Solusi Bayesian

Kita asumsikan Beta sebelumnya pada , i, e., P ~ B e t a ( α , β ) , karena model Beta-Binomial adalah konjugat, yang berarti bahwa distribusi posterior juga distribusi Beta dengan parameter α = α + k , β = β + n - k , (saya menggunakan k untuk menunjukkan jumlah keberhasilan dalam n percobaan, bukan y ). Dengan demikian, kesimpulan sangat disederhanakan. Sekarang, jika Anda memiliki pengetahuan sebelumnya tentang nilai kemungkinanppBeta(α,β)α^=α+k,β^=β+nkkny , Anda dapat menggunakannya untuk mengatur nilai α dan β , yaitu, untuk menentukan Beta Anda sebelumnya, jika tidak, Anda dapat menganggap seragam (noninformatif) sebelumnya, dengan α = β = 1 , atau prior noninformatif lainnya (lihat contohdi sini). Bagaimanapun, posterior Andahalαβα=β=1

Pr(p|n,k)=Beta(α+k,β+nk)

Dalam inferensi Bayesian, semua yang penting adalah probabilitas posterior, yang berarti bahwa setelah Anda tahu itu, Anda dapat membuat kesimpulan untuk semua kuantitas lain dalam model Anda. Anda ingin membuat inferensi pada diamati : khususnya, pada vektor hasil baru y = y 1 , ... , y m , di mana m belum tentu sama untuk n . Khususnya, untuk setiap j = 0 , ... , m , kami ingin menghitung probabilitas untuk mendapatkan keberhasilan j secara tepat dalam percobaan m berikutnya , mengingat bahwa kami mendapat kyy=y1,,ymmnj=0,,mjmkkeberhasilan dalam sebelumnya uji coba; fungsi massa prediktif posterior:n

Pr(j|m,y)=Pr(j|m,n,k)=01Pr(j,p|m,n,k)dp=01Pr(j|p,m,n,k)Pr(p|n,k)dp

Namun, model Binomial kami untuk berarti bahwa, kondisional pada p memiliki nilai tertentu, kemungkinan memiliki j keberhasilan dalam m percobaan tidak tergantung pada hasil masa lalu: itu hanyaYpjm

f(j|m,p)=(jm)pj(1p)j

Jadi ungkapan itu menjadi

Pr(j|m,n,k)=01(jm)pj(1p)jPr(p|n,k)dp=01(jm)pj(1p)jBeta(α+k,β+nk)dp

Hasil integral ini adalah distribusi terkenal yang disebut distribusi Beta-Binomial: melewatkan bagian-bagian, kita mendapatkan ekspresi yang mengerikan

Pr(j|m,n,k)=m!j!(mj)!Γ(α+β+n)Γ(α+k)Γ(β+nk)Γ(α+k+j)Γ(β+n+mkj)Γ(α+β+n+m)

Estimasi titik kami untuk , mengingat kerugian kuadratik, tentu saja adalah rata-rata dari distribusi ini, yaitu,j

μ=m(α+k)(α+β+n)

Sekarang, mari kita cari interval prediksi. Karena ini adalah distribusi diskrit, kita tidak memiliki ekspresi bentuk tertutup untuk , sehingga P r ( j 1j j 2 ) = 0,95 . Alasannya adalah bahwa, tergantung pada bagaimana Anda mendefinisikan suatu kuantil, untuk distribusi diskrit, fungsi kuantil bukanlah fungsi atau fungsi diskontinyu. Tapi ini bukan masalah besar: untuk m kecil , Anda bisa menuliskan probabilitas m P r ( j = 0[j1,j2]Pr(j1jj2)=0.95mm dan dari sini temukan j 1 , j 2 sedemikian rupa sehinggaPr(j=0|m,n,k),Pr(j1|m,n,k),,Pr(jm1|m,n,k)j1,j2

Pr(j1jj2)=Pr(jj2|m,n,k)Pr(j<j1|m,n,k)0.95

Tentu saja Anda akan menemukan lebih dari satu pasangan, sehingga Anda akan idealnya mencari terkecil seperti yang di atas adalah puas. Catat itu[j1,j2]

Pr(j=0|m,n,k)=p0,Pr(j1|m,n,k)=p1,,Pr(jm1|m,n,k)=pm1

are just the values of the CMF (Cumulative Mass Function) of the Beta-Binomial distribution, and as such there is a closed form expression, but this is in terms of the generalized hypergeometric function and thus is quite complicated. I'd rather just install the R package extraDistr and call pbbinom to compute the CMF of the Beta-Binomial distribution. Specifically, if you want to compute all the probabilities p0,,pm1 in one go, just write:

library(extraDistr)  
jvec <- seq(0, m-1, by = 1) 
probs <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)

di mana alphadan betaadalah nilai-nilai parameter Beta Anda sebelumnya, yaitu, dan β (dengan demikian 1 jika Anda menggunakan seragam sebelum lebih dari p ). Tentu saja semua akan jauh lebih sederhana jika R menyediakan fungsi kuantil untuk distribusi Beta-Binomial, tetapi sayangnya tidak.αβp

Contoh praktis dengan solusi Bayesian

Misalkan , k = 70 (dengan demikian kami awalnya mengamati 70 keberhasilan dalam 100 percobaan). Kami menginginkan estimasi titik dan interval prediksi-95% untuk jumlah keberhasilan j dalam m = 20 percobaan berikutnya. Kemudiann=100k=70jm=20

n <- 100
k <- 70
m <- 20
alpha <- 1
beta  <- 1

p

bayesian_point_estimate <- m * (alpha + k)/(alpha + beta + n) #13.92157

j

jvec <- seq(0, m-1, by = 1)
library(extraDistr)
probabilities <- pbbinom(jvec, m, alpha = alpha + k, beta = beta + n - k)

Peluangnya adalah

> probabilities
 [1] 1.335244e-09 3.925617e-08 5.686014e-07 5.398876e-06
 [5] 3.772061e-05 2.063557e-04 9.183707e-04 3.410423e-03
 [9] 1.075618e-02 2.917888e-02 6.872028e-02 1.415124e-01
[13] 2.563000e-01 4.105894e-01 5.857286e-01 7.511380e-01
[17] 8.781487e-01 9.546188e-01 9.886056e-01 9.985556e-01

j2Pr(jj2|m,n,k)0.975j1 such that Pr(j<j1|m,n,k)=Pr(jj11|m,n,k)0.025. This way, we will have

Pr(j1jj2|m,n,k)=Pr(jj2|m,n,k)Pr(j<j1|m,n,k)0.9750.025=0.95

Thus, by looking at the above probabilities, we see that j2=18 and j1=9. The probability of this Bayesian prediction interval is 0.9778494, which is larger than 0.95. We could find shorter intervals such that Pr(j1jj2|m,n,k)0.95, but in that case at least one of the two inequalities for the tail probabilities wouldn't be satisfied.

Frequentist solution

I'll follow the treatment of Krishnamoorthy and Peng, 2011. Let YBinom(m,p) and XBinom(n,p) be independently Binominally distributed. We want a 12αprediction interval for Y, based on a observation of X. In other words we look for I=[L(X;n,m,α),U(X;n,m,α)] such that:

PrX,Y(YI)=PrX,Y(L(X;n,m,α)YU(X;n,m,α)]12α

The "12α" is due to the fact that we are dealing with a discrete random variable, and thus we cannot expect to get exact coverage...but we can look for an interval which has always at least the nominal coverage, thus a conservative interval. Now, it can be proved that the conditional distribution of X given X+Y=k+j=s is hypergeometric with sample size s, number of successes in the population n and population size n+m. Thus the conditional pmf is

Pr(X=k|X+Y=s,n,n+m)=(nk)(msk)(m+ns)

The conditional CDF of X given X+Y=s is thus

Pr(Xk|s,n,n+m)=H(k;s,n,n+m)=i=0k(ni)(msi)(m+ns)

The first great thing about this CDF is that it doesn't depend on p, which we don't know. The second great thing is that it allows to easily find our PI: as a matter of fact, if we observed a value k of X, then the 1α lower prediction limit is the smallest integer L such that

Pr(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α

correspondingly, the the 1α upper prediction limit is the largest integer such that

Pr(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>α

Thus, [L,U] is a prediction interval for Y of coverage at least 12α. Note that when p is close to 0 or 1, this interval is conservative even for large n, m, i.e., its coverage is quite larger than 12α.

Practical example with the Frequentist solution

Same setting as before, but we don't need to specify α and β (there are no priors in the Frequentist framework):

n <- 100
k <- 70
m <- 20

The point estimate is now obtained using the MLE estimate for the probability of successes, p^=kn, which in turns leads to the following estimate for the number of successes in m trials:

frequentist_point_estimate <- m * k/n #14

For the prediction interval, the procedure is a bit different. We look for the largest U such that Pr(Xk|k+U,n,n+m)=H(k;k+U,n,n+m)>α, thus let's compute the above expression for all U in [0,m]:

jvec <- seq(0, m, by = 1)
probabilities <- phyper(k,n,m,k+jvec)

We can see that the largest U such that the probability is still larger than 0.025 is

jvec[which.min(probabilities > 0.025) - 1] # 18

Same as for the Bayesian approach. The lower prediction bound L is the smallest integer such that Pr(Xk|k+L,n,n+m)=1H(k1;k+L,n,n+m)>α, thus

probabilities <- 1-phyper(k-1,n,m,k+jvec)
jvec[which.max(probabilities > 0.025) - 1] # 8

Thus our frequentist "exact" prediction interval is [L,U]=[8,18].

DeltaIV
sumber