Интервал прогнозирования для биномиальной случайной величины

14

Какова формула (приблизительная или точная) для интервала предсказания для биномиальной случайной величины?

Предположим, что YBinom(n,p) , и мы наблюдаем y (взятый из Y ). n известно.

Наша цель состоит в том, чтобы получить интервал прогнозирования на 95% для нового розыгрыша от Y .

Оценка точка np^ , где р = уp^=yn . Доверительный интервал для р является простым, но я не могу найти формулу для интервала прогнозирования дляY. Если бы мы зналир(а не р ), то интервал прогнозирования 95% просто предполагает нахождение квантилей бинома. Есть ли что-то очевидное, что я пропускаю?p^Ypp^

Statseeker
источник
1
См. Какие не байесовские методы существуют для прогнозного вывода? , В этом случае метод с использованием сводок недоступен (я не думаю), но вы можете использовать одну из предиктивных вероятностей. Или, конечно, байесовский подход.
Scortchi - Восстановить Монику
1
Привет, ребята, я хотел бы воспользоваться моментом для решения проблем, которые были подняты. - относительно доверия для p: я не заинтересован в этом. - в отношении прогнозирования, составляющего 95% распределения: да, это именно то, что интервалы прогнозирования являются независимыми от контекста (в регрессии вы должны допускать нормальные ошибки, где в качестве доверительных интервалов используется CLT - да, пример прогнозирования количества голов в бросок монеты - это правильно. Что делает эту проблему трудной, так это то, что у нас нет «р», а есть оценка.
Statseeker
3
@ Аддисон. Прочитайте книгу «Статистические интервалы» Дж. Хана и В. Микера. Они объясняют разницу между доверительными интервалами, интервалами прогнозирования, интервалами допуска и байесовскими достоверными интервалами. Интервал прогнозирования 95% не содержит 95% распределения. Это делает то, что делают самые частые интервалы. Если вы неоднократно производите выборку из B (n, p) и каждый раз используете один и тот же метод для получения 95% интервала прогнозирования для p, то в 95% интервалов прогнозирования вы будете содержать истинное значение p. Если вы хотите покрыть 95% распределения, создайте интервал допуска.
Майкл Р. Черник
Интервалы допусков покрывают процент от распределения. Для интервала допуска 95% для 90% распределения вы снова повторяете процесс много раз и используете один и тот же метод для генерации интервала каждый раз, тогда примерно в 95% случаев по меньшей мере 90% распределения попадет в интервал и в 5% случаев в интервале будет содержаться менее 90% распределения.
Майкл Р. Черник
3
Lawless & Fredette (2005), «Интервалы прогнозирования для частоты и прогнозирующие распределения», Biometrika , 92 , 3 - еще одна хорошая ссылка, в дополнение к ссылкам, которые я дал.
Scortchi - Восстановить Монику

Ответы:

24

Хорошо, давайте попробуем это. Я дам два ответа - байесовский, который, на мой взгляд, простой и естественный, и один из возможных частых.

Байесовское решение

Мы предполагаем , бета перед на , я, е., Р ~ Б е т в ( & alpha ; , & beta ) , так как модель Бета-биномиальное сопряжена, что означает , что задняя распределение также бета - распределение с параметрами α = α + к , β = β + п - к , (я использую K , чтобы обозначить число успехов в п испытаниях, вместо того , у ). Таким образом, вывод значительно упрощается. Теперь, если у вас есть некоторые предварительные знания о вероятных значенияхppBeta(α,β)α^=α+k,β^=β+nkknyp, you could use it to set the values of α and β, i.e., to define your Beta prior, otherwise you could assume a uniform (noninformative) prior, with α=β=1, or other noninformative priors (see for example here). In any case, your posterior is

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

В байесовском умозаключении все, что имеет значение, это апостериорная вероятность, означающая, что, как только вы это знаете, вы можете делать выводы для всех других величин в вашей модели. Вы хотите сделать вывод о наблюдаемых : в частности, о векторе новых результатов y = y 1 , , y m , где m не обязательно равно n . В частности, для каждого j = 0 , , m мы хотим вычислить вероятность достижения именно j успехов в следующих m испытаниях, учитывая, что мы получили kyy=y1,,ymmnj=0,,mjmkуспехи в предыдущих испытаниях; задняя предиктивная функция массы: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

However, our Binomial model for Y means that, conditionally on p having a certain value, the probability of having j successes in m trials doesn't depend on past results: it's simply

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

Thus the expression becomes

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

The result of this integral is a well-known distribution called the Beta-Binomial distribution: skipping the passages, we get the horrible expression

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

Our point estimate for j, given quadratic loss, is of course the mean of this distribution, i.e.,

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

Now, let's look for a prediction interval. Since this is a discrete distribution, we don't have a closed form expression for [j1,j2], such that Pr(j1jj2)=0.95. The reason is that, depending on how you define a quantile, for a discrete distribution the quantile function is either not a function or is a discontinuous function. But this is not a big problem: for small m, you can just write down the m probabilities Pr(j=0|m,n,k),Pr(j1|m,n,k),,Pr(jm1|m,n,k) and from here find j1,j2 such that

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

Of course you would find more than one couple, so you would ideally look for the smallest [j1,j2] such that the above is satisfied. Note that

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)

where alpha and beta are the values of the parameters of your Beta prior, i.e., α and β (thus 1 if you're using a uniform prior over p). Of course it would all be much simpler if R provided a quantile function for the Beta-Binomial distribution, but unfortunately it doesn't.

Practical example with the Bayesian solution

Let n=100, k=70 (thus we initially observed 70 successes in 100 trials). We want a point estimate and a 95%-prediction interval for the number of successes j in the next m=20 trials. Then

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

where I assumed a uniform prior on p: depending on the prior knowledge for your specific application, this may or may not be a good prior. Thus

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

Clearly a non-integer estimate for j doesn't make sense, so we could just round to the nearest integer (14). Then, for the prediction interval:

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

The probabilities are

> 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

For an equal-tail probabilities interval, we want the smallest j2 such that Pr(jj2|m,n,k)0.975 and the largest j1 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
источник