Как найти доверительный интервал для общего количества событий

У меня есть детектор, который обнаружит событие с некоторой вероятностью р . Если детектор сообщает, что событие произошло, то это всегда так, поэтому ложных срабатываний нет. Через некоторое время я получаю k обнаруженных событий. Я хотел бы подсчитать, каково было общее количество событий, которые произошли, обнаружены или нет, с некоторой достоверностью, скажем, 95%.

Например, допустим, я обнаружил 13 событий. Я хотел бы иметь возможность рассчитать, что было от 13 до 19 событий с 95% достоверностью на основе p .

Вот что я пробовал до сих пор:

Вероятность обнаружения k событий, если их было всего n, равна:

binomial(n, k) * p^k * (1 - p)^(n - k)

Сумма этого по n от k до бесконечности равна:

1/p

Это означает, что вероятность того, что всего будет n событий, равна:

f(n) = binomial(n, k) * p^(k + 1) * (1 - p)^(n - k)

Так что, если я хочу быть 95% уверен , что я должен найти первую частичную сумму , f(k) + f(k+1) + f(k+2) ... + f(k+m)которая составляет , по меньшей мере , 0,95 и ответ [k, k+m]. Это правильный подход? Также есть закрытая формула для ответа?

Я бы предпочел использовать отрицательное биномиальное распределение , которое возвращает вероятность того, что будет X неудач до k-го успеха, когда постоянная вероятность успеха равна p.

Используя пример

k=17 # number of successes
p=.6 # constant probability of success

среднее значение и сд для сбоев задаются

mean.X <- k*(1-p)/p
sd.X <- sqrt(k*(1-p)/p^2) 

Распределение сбоев X будет иметь примерно такую ​​форму

plot(dnbinom(0:(mean.X + 3 * sd.X),k,p),type='l')

Таким образом, количество отказов будет (с доверительной вероятностью 95%) примерно между

qnbinom(.025,k,p)
[1] 4

а также

qnbinom(.975,k,p)
[1] 21

Таким образом, вы будете иметь значение [k + qnbinom (.025, k, p), k + qnbinom (.975, k, p)] (используя числа примера [21,38])

Предполагая, что вы хотите выбрать распределение для n, p (n), вы можете применить закон Байеса.

Вы знаете, что вероятность k событий, происходящих при условии, что n действительно произошло, определяется биномиальным распределением

$p(k|n) = {n \choose k} p^k (1-p)^{(n-k)}$

То, что вы действительно хотите знать, это вероятность того, что n событий действительно произошло, учитывая, что вы наблюдали k. По Байесу выложу

$p(n|k) = \frac{p(k|n)p(n)}{p(k)}$

Применяя теорему полной вероятности, мы можем написать:

$p(n|k) = \frac{p(k|n)p(n)}{\sum_{n'} p(k|n')p(n')}$

Так что без дополнительной информации о распределении вы не сможете пойти дальше.$p(n)$

Однако, если вы хотите выбрать распределение для для которого есть значение большее, чем , или достаточно близкое к нулю, вы можете сделать это немного лучше. Например, предположим, что распределение является равномерным в диапазоне . этот случай:$p(n)$$n$$p(n) = 0$$n$$[0,n_{max}]$

$p(n) = \frac{1}{n_{max}}$

Байесовская формулировка упрощает:

$p(n|k) = \frac{p(k|n)}{\sum_{n'} p(k|n')}$

Что касается последней части проблемы, я согласен с тем, что наилучшим подходом является выполнение кумулятивного суммирования по , генерация кумулятивной функции распределения вероятностей и итерация до достижения предела 0,95.$p(n|k)$

Учитывая, что этот вопрос перенесен из SO, пример кода игрушки на python прилагается ниже

import numpy.random

p = 0.8
nmax = 200

def factorial(n):
    if n == 0:
        return 1
    return reduce( lambda a,b : a*b, xrange(1,n+1), 1 )

def ncr(n,r):
    return factorial(n) / (factorial(r) * factorial(n-r))

def binomProbability(n, k, p):
    p1 = ncr(n,k)
    p2 = p**k
    p3 = (1-p)**(n-k)
    return p1*p2*p3

def posterior( n, k, p ):
    def p_k_given_n( n, k ):
        return binomProbability(n, k, p)
    def p_n( n ):
        return 1./nmax
    def p_k( k ):
        return sum( [ p_n(nd)*p_k_given_n(nd,k) for nd in range(k,nmax) ] )
    return (p_k_given_n(n,k) * p_n(n)) / p_k(k)


observed_k   = 80
p_n_given_k  = [ posterior( n, observed_k, p ) for n in range(0,nmax) ]
cp_n_given_k = numpy.cumsum(p_n_given_k)
for n in xrange(0,nmax):
    print n, p_n_given_k[n], cp_n_given_k[n]

Если вы измеряете событий и знаете, что ваша эффективность обнаружения равна вы можете автоматически скорректировать свой измеренный результат до «истинного» счета .$k$$p$$k_\mathrm{true} = k/p$

Ваш вопрос заключается в том, чтобы найти диапазон в который попадет 95% наблюдений. Вы можете использовать метод Фельдмана-Казинса, чтобы оценить этот интервал. Если у вас есть доступ к ROOT, есть класс для этого расчета.$k_\mathrm{true}$

Вы можете рассчитать верхний и нижний пределы с Фельдманом-кузенами из нескорректированного числа событий а затем масштабировать их до 100% с помощью . Таким образом, фактическое количество измерений определяет вашу неопределенность, а не какое-то масштабированное число, которое не было измерено.$k$$1/p$

{
gSystem->Load("libPhysics");

const double lvl = 0.95;
TFeldmanCousins f(lvl);

const double p = 0.95;
const double k = 13;
const double k_true = k/p;

const double k_bg = 0;

const double upper = f.CalculateUperLimit(k, k_bg) / p;
const double lower = f.GetLowerLimit() / p;

std::cout << "["
  lower <<"..."<<
  k_true <<"..."<<
  upper <<
  "]" << std::endl;
}

Я думаю, что вы неправильно поняли цель доверительных интервалов. Доверительные интервалы позволяют оценить, где находится истинное значение параметра. Итак, в вашем случае вы можете построить доверительный интервал для . Не имеет смысла строить интервал для данных.$p$

Сказав, что, получив оценку вы можете вычислить вероятность того, что вы будете наблюдать различные реализации, такие как 14, 15 и т. Д., Используя биномиальный pdf.$p$