;$Id: binomial.pro,v 1.12 2004/01/21 15:54:48 scottm Exp $ ; ; Copyright (c) 1994-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;+ ; NAME: ; BINOMIAL ; ; PURPOSE: ; This function computes the probabilty (bp) such that: ; Probability(X => v) = bp ; where X is a random variable from the cumulative binomial distribution ; (Bernouli distribution). ; ; CATEGORY: ; Statistics. ; ; CALLING SEQUENCE: ; Result = Binomial(V, N, P) ; ; INPUTS: ; V: A non-negative integer specifying the minimal number of ; times an event E occurs in (N) independent performances. ; ; N: A non-negative integer specifying the number of performances. ; ; P: A non-negative scalar or array, in the interval [0.0, 1.0], ; of type float or double that specifies the probability of ; occurrence or success of a single independent performance. ; ; KEYWORDS: ; ; DOUBLE = Set this keyword to force the computation to be done in ; double-precision arithmetic. ; ; GAUSSIAN = Set this keyword to force the computation to be done using ; the Gaussian approximation. ; ; EXAMPLES: ; Compute the probability of obtaining at least two 6s in rolling a ; die four times. The result should be 0.131944 ; result = binomial(2, 4, 1./6.) ; ; Compute the probability of obtaining exactly two 6s in rolling a ; die four times. The result should be 0.115741 ; result = binomial(2, 4, 1./6.) - binomial(3, 4, 1./6.) ; ; Compute the probability of obtaining three or fewer 6s in rolling ; a die four times. The result should be 0.999228 ; result = (binomial(0, 4, 1./6.) - binomial(1, 4, 1./6.)) + $ ; (binomial(1, 4, 1./6.) - binomial(2, 4, 1./6.)) + $ ; (binomial(2, 4, 1./6.) - binomial(3, 4, 1./6.)) + $ ; (binomial(3, 4, 1./6.) - binomial(4, 4, 1./6.)) ; ; PROCEDURE: ; BINOMIAL computes the probability that an event E occurs at least ; (V) times in (N) independent performances. The event E is assumed ; to have a probability of occurance or success (P) in a single ; performance. ; ; If an overflow occurs during computation, then the Gaussian ; distribution is used to approximate the cumulative binomial ; distribution. ; ; REFERENCE: ; ADVANCED ENGINEERING MATHEMATICS (seventh edition) ; Erwin Kreyszig ; ISBN 0-471-55380-8 ; ; Schaum's Outline of Theory and Problems of Probability and Statistics, ; M.R. Spiegel, McGraw-Hill, 1975. ; ; MODIFICATION HISTORY: ; Modified by: GGS, RSI, July 1994 ; Minor changes to code. Rewrote documentation header. ; CT, RSI, June 2000: Added keywords DOUBLE, GAUSSIAN; ; changed algorithm to use LNGAMMA function; ; doesn't use Gaussian unless overflow or /GAUSSIAN; ; now allows array input for P. ;- function binomial, v, n, pIn, $ DOUBLE=double, $ GAUSSIAN=gaussian on_error, 2 ;Return to caller if error occurs. if MIN(pIn) lt 0. or MAX(pIn) gt 1. then message, $ 'P must be in the interval [0.0, 1.0]' if v lt 0 then message, 'V must be nonnegative.' if n lt 0 then message, 'N must be nonnegative.' double = (N_ELEMENTS(double) GT 0) ? KEYWORD_SET(double) : $ (SIZE(pIn,/TNAME) EQ 'DOUBLE') p = double ? DOUBLE(pIn) : FLOAT(pIn) sum = p*0 IF (v EQ 0) THEN RETURN, sum+1 ; one IF (v GT n) THEN RETURN, sum ; zero q = 1 - p nn = fix(n) vv = fix(v) nn = double ? DOUBLE(nn) : FLOAT(nn) vv = double ? DOUBLE(vv) : FLOAT(vv) ; for accuracy use logarithms to evaluate factorials logFact = LNGAMMA(nn+1)-LNGAMMA(vv+1)-LNGAMMA(nn-vv+1) ; if we are going to overflow, then switch to Gaussian approx overflow = logFact GT (double ? 700 : 80) do_gauss = (N_ELEMENTS(gaussian) EQ 0) AND overflow IF KEYWORD_SET(gaussian) OR do_gauss THEN BEGIN ; use Gaussian sum = 1 - GAUSS_PDF((v-0.5-n*p)/SQRT(n*p*q)) ENDIF ELSE BEGIN ; use the binomial summation previous = (sum = EXP(logFact) * p^vv * q^(nn-vv)) zero = WHERE(q EQ 0, nZero) IF (nZero GT 0) THEN q[zero] = 1 ; avoid divide by zero FOR x = vv+1, nn DO BEGIN ; ratio of next term in sequence to previous previous = previous*(nn-x+1)/x*(p/q) sum = sum + previous ENDFOR IF (nZero GT 0) THEN sum[zero] = 1 ; if p=1 then sum=1 ENDELSE RETURN, sum end