The R Project SVN R

Rev

Rev 12976 | Rev 42302 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 12976 Rev 30227
Line 25... Line 25...
25 
 *
 25 
 *
26 
 *    Given a sequence of r successes and b failures, we sample n (\le b+r)
 26 
 *    Given a sequence of r successes and b failures, we sample n (\le b+r)
27 
 *    items without replacement. The hypergeometric probability is the
 27 
 *    items without replacement. The hypergeometric probability is the
28 
 *    probability of x successes:
 28 
 *    probability of x successes:
29 
 *
 29 
 *
- 
 
 30 
 *		       choose(r, x) * choose(b, n-x)
- 
 
 31 
 *	p(x; r,b,n) =  -----------------------------  =
- 
 
 32 
 *			       choose(r+b, n)
- 
 
 33 
 *
30 
 *                   dbinom(x,r,p) * dbinom(n-x,b,p)
 34 
 *		      dbinom(x,r,p) * dbinom(n-x,b,p)
31 
 *      p(x;r,b,n) = ---------------------------------
 35 
 *		    = --------------------------------
32 
 *                             dbinom(n,r+b,p)
 36 
 *			       dbinom(n,r+b,p)
33 
 *
 37 
 *
34 
 *    for any p. For numerical stability, we take p=n/(r+b); with this choice,
 38 
 *    for any p. For numerical stability, we take p=n/(r+b); with this choice,
35 
 *    the denominator is not exponentially small.
 39 
 *    the denominator is not exponentially small.
36 
 */
 40 
 */
37 
 
 41 
 
38 
#include "nmath.h"
 42 
#include "nmath.h"
39 
#include "dpq.h"
 43 
#include "dpq.h"
40 
 
 44 
 
41 
double dhyper(double x, double r, double b, double n, int give_log)
 45 
double dhyper(double x, double r, double b, double n, int give_log)
42 
{
46 
{
43 
    double p, q, p1, p2, p3;
 47 
    double p, q, p1, p2, p3;
44 
 
 48 
 
45 
#ifdef IEEE_754
 49 
#ifdef IEEE_754
46 
    if (ISNAN(x) || ISNAN(r) || ISNAN(b) || ISNAN(n))
 50 
    if (ISNAN(x) || ISNAN(r) || ISNAN(b) || ISNAN(n))
47 
        return x + r + b + n;
 51 
	return x + r + b + n;
48 
#endif
 52 
#endif
49 
 
 53 
 
50 
    if (R_D_notnnegint(r) || R_D_notnnegint(b) || R_D_notnnegint(n) || n > r+b)
 54 
    if (R_D_negInonint(r) || R_D_negInonint(b) || R_D_negInonint(n) || n > r+b)
51 
	ML_ERR_return_NAN;
 55 
	ML_ERR_return_NAN;
- 
 
 56 
    if (R_D_negInonint(x))
- 
 
 57 
	return(R_D__0);
52 
 
 58 
 
53 
    if (R_D_notnnegint(x)) return(R_D__0);
 - 
 
54 
    x = R_D_forceint(x);
 59 
    x = R_D_forceint(x);
55 
    r = R_D_forceint(r);
 60 
    r = R_D_forceint(r);
56 
    b = R_D_forceint(b);
 61 
    b = R_D_forceint(b);
57 
    n = R_D_forceint(n);
 62 
    n = R_D_forceint(n);
58 
 
 63 
 
Line 60... Line 65...
60 
    if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
 65 
    if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
61 
 
 66 
 
62 
    p = ((double)n)/((double)(r+b));
 67 
    p = ((double)n)/((double)(r+b));
63 
    q = ((double)(r+b-n))/((double)(r+b));
 68 
    q = ((double)(r+b-n))/((double)(r+b));
64 
 
 69 
 
65 
    p1 = dbinom_raw(x,  r, p,q,give_log);
 70 
    p1 = dbinom_raw(x,	r, p,q,give_log);
66 
    p2 = dbinom_raw(n-x,b, p,q,give_log);
 71 
    p2 = dbinom_raw(n-x,b, p,q,give_log);
67 
    p3 = dbinom_raw(n,r+b, p,q,give_log);
 72 
    p3 = dbinom_raw(n,r+b, p,q,give_log);
68 
 
 73 
 
69 
    return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
 74 
    return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
70 
}
 75 
}