; $Id: binomial.pro,v 1.1 1993/04/02 18:54:39 idl Exp $

;  Copyright (c) 1991-1993, Research Systems Inc.  All rights
;  reserved. Unauthorized reproduction prohibited.


 function factorial,n,min,fac1
; if min and fac1 are undefined, then factorial returns n!.
; Otherwise, fac1 * min * (min+1)....n is returned.

  if N_Elements(min) eq 0 then min = 2
  if N_Elements(fac1) EQ 0 THEN fac = 1.00  $
    ELSE fac = fac1

  if min GT n then return, fac

  n1 =  n < 10

 if min lt 11 then  $
  for i = min, n1 DO fac= i*fac

 if (n lt 11.0 ) then return,fac

 n1 = 11 > min         ;use logs to preserve precision

 sum = alog(findgen(n - n1 +1) + n1)
 sum = Total(sum)
  return, fac * exp(sum)
end


function binomial,x,n,p 
;+
; NAME:
;	BINOMIAL
;
; PURPOSE: 
;	Binomial implements the cumulative binomial distribution.
;
; CATEGORY:
;	Statistics.
;
; CALLING SEQUENCE:
;	Result = BINOMIAL(X, N, P) 
;
; INPUTS:
;	X:	The number of successes.
;
;	N:	The number of trials.
;
;	P:	The success probability.
;
; OUTPUT:
;	Binomial returns the probability of x or more successes in a binomial
;	experiment with n trials and success probability p. If n exceeds 25,
;	the Gaussian approximation is computed.
; 
;-
 if p lt 0 or p gt 1 THEN BEGIN
    print,'binomial -- must have 0<=p<=1'
    return,-1
 ENDIF

 if n lt 0 or x lt 0 then BEGIN
   printf,'binomial -- number of successes and number of' 
   printf, 'trials must be nonnegative'
   return, -1
 ENDIF

 if x EQ 0 THEN return,1
 if n gt 25 then return, 1 - gaussint((x - n*p)/sqrt(n*p*(1-p)))
 if x gt n then return,0

 nn = fix(n)
 xx = fix(x)
 n2 = xx < (nn-xx)
 n3 = xx > (nn-xx)
n1=factorial(nn,n3+1,p)
n1 =  n1/factorial(n2)
sum = n1 * p^(xx-1)*(1-p)^(nn-xx)
for i = xx+1,nn do BEGIN
 n1 = (nn-i+1) * n1/float(i)
 sum = sum + n1	 *p^(i-1) * (1-p)^(nn-i)
ENDFOR
return, sum
end