Qual é a fórmula (aproximada ou exata) para um intervalo de previsão para uma variável aleatória binomial?
Suponha , e observa- (desenhado a partir de ). O é conhecido.Y n
O nosso objectivo é a obtenção de um intervalo de previsão de 95% para um novo desenho de .
A estimativa pontual é , onde . Um intervalo de confiança para é simples, mas não consigo encontrar uma fórmula para um intervalo de previsão para . Se soubéssemos (em vez de ), um intervalo de previsão de 95% envolve apenas encontrar os quantis de um binômio. Existe algo óbvio que estou ignorando?p = yP YP P
confidence-interval
binomial
prediction-interval
Statseeker
fonte
fonte
Respostas:
Ok, vamos tentar isso. Darei duas respostas - a bayesiana, que em minha opinião é simples e natural, e uma das possíveis freqüentadoras.
Solução Bayesiana
Assumimos um antes Beta emp , i, e., p∼Beta(α,β) , porque o modelo de Beta-binomial é conjugado, o que significa que a distribuição a posteriori é também uma distribuição beta com parâmetros α = α + k , β = β + n - k , (estou usando k para indicar o número de sucessos em n ensaios, em vez de y ). Assim, a inferência é bastante simplificada. Agora, se você tem algum conhecimento prévio sobre os valores prováveis deα^=α+k,β^=β+n−k k n y p , você pode usá-lo para definir os valores deα eβ , ou seja, para definir seu Beta anterior, caso contrário, você pode assumir um anterior uniforme (não informativo), comα=β=1 ou outros anteriores não informativos (veja, por exemplo,aqui)) Em qualquer caso, seu posterior é
Na inferência bayesiana, tudo o que importa é a probabilidade posterior, significando que, uma vez que você saiba disso, poderá fazer inferências para todas as outras quantidades em seu modelo. Você deseja fazer inferência nos observáveisy : em particular, em um vetor de novos resultados y=y1,…,ym , em que m não é necessariamente igual a n . Especificamente, para cada j=0,…,m , queremos calcular a probabilidade de ter exatamente j sucessos nos próximos m testes, considerando que obtivemos k sucessos nas n tentativas anteriores ; a função de massa preditiva posterior:
Entretanto, nosso modelo binomial para significa que, condicionalmente em p ter um certo valor, a probabilidade de ter j sucessos em m ensaios não depende de resultados passados: é simplesmenteY p j m
Assim, a expressão se torna
O resultado dessa integral é uma distribuição bem conhecida chamada distribuição beta-binomial: pulando as passagens, obtemos a expressão horrível
Nossa estimativa pontual para , dada a perda quadrática, é obviamente a média dessa distribuição, ou seja,j
Agora, vamos procurar um intervalo de previsão. Como essa é uma distribuição discreta, não temos uma expressão de forma fechada para , de modo que P r ( j 1 ≤ j ≤ j 2 ) = 0,95 . O motivo é que, dependendo de como você define um quantil, para uma distribuição discreta, a função quantil não é uma função ou é uma função descontínua. Mas esse não é um grande problema: para m pequeno , basta escrever as m probabilidades P r ( j = 0[j1,j2] Pr(j1≤j≤j2)=0.95 m m Pr(j=0|m,n,k),Pr(j≤1|m,n,k),…,Pr(j≤m−1|m,n,k) and from here find j1,j2 such that
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
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 packagep0,…,pm−1 in one go, just write:
extraDistr
and callpbbinom
to compute the CMF of the Beta-Binomial distribution. Specifically, if you want to compute all the probabilitieswhereα 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.
alpha
andbeta
are the values of the parameters of your Beta prior, i.e.,Practical example with the Bayesian solution
Letn=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
where I assumed a uniform prior onp : depending on the prior knowledge for your specific application, this may or may not be a good prior. Thus
Clearly a non-integer estimate forj doesn't make sense, so we could just round to the nearest integer (14). Then, for the prediction interval:
The probabilities are
For an equal-tail probabilities interval, we want the smallestj2 such that Pr(j≤j2|m,n,k)≥0.975 and the largest j1 such that Pr(j<j1|m,n,k)=Pr(j≤j1−1|m,n,k)≤0.025 . This way, we will have
Thus, by looking at the above probabilities, we see thatj2=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(j1≤j≤j2|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. LetY∼Binom(m,p) and X∼Binom(n,p) be independently Binominally distributed. We want a 1−2α− 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:
The "≥1−2α " 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
The conditional CDF ofX given X+Y=s is thus
The first great thing about this CDF is that it doesn't depend onp , 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
correspondingly, the the1−α upper prediction limit is the largest integer such that
Thus,[L,U] is a prediction interval for Y of coverage at least 1−2α . 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 1−2α .
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):
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:
For the prediction interval, the procedure is a bit different. We look for the largestU such that Pr(X≤k|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] :
We can see that the largestU such that the probability is still larger than 0.025 is
Same as for the Bayesian approach. The lower prediction boundL is the smallest integer such that Pr(X≥k|k+L,n,n+m)=1−H(k−1;k+L,n,n+m)>α , thus
Thus our frequentist "exact" prediction interval is[L,U]=[8,18] .
fonte