How to find a confidence interval for the total number of events

Question

I have detector which will detect an event with some probability p. If the detector says that an event occured, then that is always the case, so there are not false-positives. After I run it for some time, I get k events detected. I would like to calculate what the total number of events that occured was, detected or otherwise, with some confidence, say 95%.

So for example, let's say I get 13 events detected. I would like to be able to calculate that there were between 13 and 19 events with 95% confidence based on p.

Here's what I've tried so far:

The probability of detecting k events if there were n total is:

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

The sum of that over n from k to infinity is:

1/p

Which means, that the probability of there being n events total is:

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

So if I want to be 95% sure I should find the first partial sum f(k) + f(k+1) + f(k+2) ... + f(k+m) which is at least 0.95 and the answer is [k, k+m]. Is this the correct approach? Also is there a closed formula for the answer?

Answer 1

I would choose to use the negative binomial distribution, which returns the probability that there will be X failures before the k_th success, when the constant probability of a success is p.

Using an example

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

the mean and sd for the failures are given by

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

The distribution of the failures X, will have approximately that shape

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

So, the number of failures will be (with 95% confidence) approximately between

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

and

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

So you inerval would be [k+qnbinom(.025,k,p),k+qnbinom(.975,k,p)] (using the example's numbers [21,38] )

Answer 2

Assuming you want to pick a distribution for n, p(n) you can apply Bayes law.

You know that the probability of k events occuring given that n have actually occured is governed by a binomial distribtion

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

The thing you really want to know is the probability of n events having actually occured, given that you observed k. By Bayes lay:

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

By applying the theorem of total probability, we can write:

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

So without further information, about the distribution of $p(n)$ you can't really go any further.

However, if you want to pick a distribution for $p(n)$ for which there is a value $n$ greater than which $p(n) = 0$, or sufficiently close to zero, then you can do a bit better. For example, assume that the distribution of $n$ is uniform in the range $[0,n_{max}]$. this case:

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

The Bayesian formulation simplifies to:

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

As for the final part of the problem, I agree that the best approach is to perform a cumulative summation over $p(n|k)$, to generate the cummulative probability distribution function, and iterate until the 0.95 limit is reached.

Given that this question migrated from SO, toy sample code in python is attached below

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]

Answer 3

If you measure $k$ events and know your detection efficiency is $p$ you can automatically correct your measured result up to the "true" count $k_\mathrm{true} = k/p$.

Your question is then about finding the range of $k_\mathrm{true}$ where 95% of the observations will fall. You can use the Feldman-Cousins method to estimate this interval. If you have access to ROOT there is a class to do this calculation for you.

You would calculate the upper and lower limits with Feldman-Cousins from the uncorrected number of events $k$ and then scale them up to 100% with $1/p$. This way the actual number of measurements determines your uncertainty, not some scaled number that wasn't measured.

{
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;
}

Answer 4

I think you misunderstood the purpose of confidence intervals. Confidence intervals allow you to assess where the true value of the parameter is located. So, in your case, you can construct a confidence interval for $p$. It does not make sense to construct an interval for the data.

Having said that, once you have an estimate of $p$ you can calculate the probability that you will observe different realizations such as 14, 15 etc using the binomial pdf.