Cumulative binomial distribution

| Forums | Site Guide |

Although the text wraps on this page, if you copy and paste it into an editor it should come out with correct formatting

/**
 * Program to compute cumulative values in the binomial probability distribution
 * Written by Matt Squire (June 2008) :: http://www.InsideReality.net
 * Source code made available for educational purposes only
 * No guarantees are made about the accuracy or correctness of this program.
 * In particular, I'm not sure what the maximum inputs which it can handle are
 *
 * Usage: binomial n k p
 *
 * ----------------------------------------------------------------------------------------
 *
 * For n independent trials each with probability p of success, and q = 1 - p of failure,
 * the probability of k successes is given by
 *
 *	Pr(K = k) =	nCk * p^k * q^(n-k)			(for integer k)
 *
 * 		     	     n!
 *	Where nCk =	-----------  (the number of ways to select k successes from n trials)
 *			 (n-k)!k!
 *
 * The cumulative distribution is the sum of all the probabilities K = 0 up to K = k
 *
 * 			  k
 * 			 ---
 *	Pr(K <= k) =	 \	nCj * p^j * q^(n-j)		(for integer k)
 *			 /
 *			 ---
 *			j = 0
 */
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
typedef unsigned long long int_64;

int_64 combinations(unsigned, unsigned);
double probability(unsigned, unsigned, double);

int main(int argc, char* argv[])
{
	if(argc < 4)
	{
		printf("Usage: binomial n k p\n");
		return 1;
	}
		
	int n = atoi(argv[1]);
	if(n < 1)
	{
		printf("n should be at least 1\n");
		return 1;
	}
		
	int k = atoi(argv[2]);
	if(k < 0 || k > n)
	{
		printf("k should be at least 1 and at most n\n");
		return 1;
	}
		
	double p = atof(argv[3]);
	if(p < 0 || p > 1)
	{
		printf("p should be between 0 and 1 (inclusive)\n");
		//return 1;
	}
	
	//Input is good, so time to crunch some numbers
	printf("Input:\n");
	printf("\tn = %d\n\tk = %d\n\tp = %0.3f\n", n, k, p);
	printf("---------\n");
	printf("Pr(K = %d) = %0.4f\n", k, probability(n, k, p));
	
	//Sum probabilities over desired range
	double sum = 0;
	for(int i = 0; i <= k; i++)
		sum += probability(n, i, p);
		
	printf("Pr(K <= %d) = %0.4f\n", k, sum);
}

/**
 * Calculates Pr(K = k) = nCk * p^k * q^(n-k)
 */
double probability(unsigned n, unsigned k, double p)
{
	if(p < 0 || p > 1)
		return -1;
	
	return combinations(n, k) * pow(p, k) * pow(1 - p, n - k);
}

/**
 * Calculates nCk = n! / (n-k)!k!
 */
 int_64 combinations(unsigned n, unsigned k)
 {
 	if(k == 0 || k == n)
 		return 1;
 		
 	if(k > n)
 		return 0;
 	
 	//After a couple of attempts at an efficient solution, I just decided to use this,
 	//from http://en.wikipedia.org/wiki/Binomial_coefficient
 	//Much faster than recursion and doesn't overflow easily
	if (k > n/2)
		k = n-k;

	long double p = 1;
	for (unsigned i = 1; i <= k; i++)
		p = p * (n-k+i) / i;

	return p + 0.5;
 }