If I make N trials, each independent with p chance of success, what is the probability that X or more of them will be successful?

Question

In my problem, $p$, $X$, and $N$ are fixed and I want to know the probability. Currently, I am calculating this recursively, as for any given $p$, $X$, $N$, the probability is:

function Prob(p,X,N) {
  return p*Prob(p,X-1,N-1) + (1-p)*Prob(p,X,N-1)
}

The intuition behind this is that if we were successful (first term) we have $X-1$ more successes to find and $N-1$ more trials. If we weren't successful $(1-p)$, we have $X$ more successes to find and $N-1$ more trials. There are also base cases that I left out where $X=0$ (probability is 1 because we've gotten all the successes) and $N$ less than $X$ (probability is zero because we have no more trials)
The recursion is really killing the speed of computation, so I'm looking for a closed form or approximation or just a point in the right direction.
Thanks

Answer 1

Most programming languages have "Cumulative Distribution Functions" (CDF) for various distributions. As a result, you shouldn't have to individually compute each probability and then sum up to $X$. For example, in R, to obtain the CDF of a binomial distribution with $p=.20$, $n=50$, and $X=10$, you don't need to calculate:

> sum(dbinom(x=0:10, size=50, prob=0.20))
[1] 0.5835594

Instead, you could simply use:

> pbinom(q=10, size=50, prob=0.20)
[1] 0.5835594

Comparing the times, you see that the built in CDF function is generally faster:

> Sys.time()->start;
> sum(dbinom(0:54, size=100, prob=.50) )
[1] 0.8158992
> print(Sys.time()-start);
Time difference of 0.002012014 secs
> 
> Sys.time()->start;
> pbinom(q=54, size=100, prob=0.50)
[1] 0.8158992
> print(Sys.time()-start);
Time difference of 0.0009989738 secs

Also, if $n$ is large and/or $p$ is close to 0.50, you can use a normal approximation to the binomial. For example, suppose you wanted to calculate the Probability of at least 41 successes in 80 tries when $p=0.5$. You could calculate:

n=80
x=41
p=.5
q=1-p

pbinom(n,size=n,prob=p)-pbinom(x,size=n,prob=p)
[1] 0.3687772

or you could use the normal approximation:

1-pnorm(x,mean=n*p,sd=sqrt(n*p*q))
0.4115316

and with continuity correction, you get:

1-pnorm(x+.5,mean=n*p,sd=sqrt(n*p*q))
0.3686578