The R Project SVN R

Rev

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

Rev 8431 Rev 11723
Line 1... Line 1...
1 
/*
 1 
/*
- 
 
 2 
 *  AUTHOR
2 
 *  Mathlib : A C Library of Special Functions
 3 
 *    Catherine Loader, catherine@research.bell-labs.com.
- 
 
 4 
 *    October 23, 2000.
- 
 
 5 
 *
3 
 *  Copyright (C) 1998 Ross Ihaka
 6 
 *  Merge in to R:
4 
 *  Copyright (C) 2000 The R Development Core Team
 7 
 *	Copyright (C) 2000, The R Core Development Team
5 
 *
 8 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 9 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 10 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 11 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 12 
 *  (at your option) any later version.
Line 15... Line 18...
15 
 *
 18 
 *
16 
 *  You should have received a copy of the GNU General Public License
 19 
 *  You should have received a copy of the GNU General Public License
17 
 *  along with this program; if not, write to the Free Software
 20 
 *  along with this program; if not, write to the Free Software
18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
 21 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
19 
 *
 22 
 *
20 
 *  DESCRIPTION
 - 
 
21 
 *
 23 
 *
- 
 
 24 
 * DESCRIPTION
- 
 
 25 
 *
- 
 
 26 
 *    Given a sequence of r successes and b failures, we sample n (\le b+r)
22 
 *    The density of the hypergeometric distribution.
 27 
 *    items without replacement. The hypergeometric probability is the
- 
 
 28 
 *    probability of x successes:
- 
 
 29 
 *
- 
 
 30 
 *                   dbinom(x,r,p) * dbinom(n-x,b,p)
- 
 
 31 
 *      p(x;r,b,n) = ---------------------------------
- 
 
 32 
 *                             dbinom(n,r+b,p)
- 
 
 33 
 *
- 
 
 34 
 *    for any p. For numerical stability, we take p=n/(r+b); with this choice,
- 
 
 35 
 *    the denominator is not exponentially small.
23 
 */
 36 
 */
24 
 
 37 
 
25 
#include "nmath.h"
 38 
#include "nmath.h"
26 
#include "dpq.h"
 39 
#include "dpq.h"
27 
 
 40 
 
28 
double dhyper(double x, double NR, double NB, double n, int give_log)
 41 
double dhyper(double x, double r, double b, double n, int give_log)
29 
{
 42 
{
30 
    double N;
 43 
    double p, q, p1, p2, p3;
- 
 
 44 
 
31 
#ifdef IEEE_754
 45 
#ifdef IEEE_754
32 
    if (ISNAN(x) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
 46 
    if (ISNAN(x) || ISNAN(r) || ISNAN(b) || ISNAN(n))
33 
	return x + NR + NB + n;
 47 
        return x + r + b + n;
34 
#endif
 48 
#endif
35 
    x = floor(x + 0.5);
 - 
 
36 
    NR = floor(NR + 0.5);
 - 
 
37 
    NB = floor(NB + 0.5);
 - 
 
38 
    N = NR + NB;
 - 
 
39 
    n = floor(n + 0.5);
 - 
 
- 
 
 49 
 
40 
    if (NR < 0 || NB < 0 || n < 0 || n > N) {
 50 
    if (R_D_notnnegint(r) || R_D_notnnegint(b) || R_D_notnnegint(n) || n > r+b)
41 
	ML_ERR_return_NAN;
 51 
	ML_ERR_return_NAN;
42 
    }
 52 
 
43 
    if (x < fmax2(0, n - NB) || x > fmin2(n, NR))
 53 
    if (R_D_notnnegint(x)) return(R_D__0);
- 
 
 54 
    x = R_D_forceint(x);
- 
 
 55 
    r = R_D_forceint(r);
- 
 
 56 
    b = R_D_forceint(b);
44 
	return R_D__0;
 57 
    n = R_D_forceint(n);
- 
 
 58 
 
45 
    return R_D_exp(lfastchoose(NR, x) + lfastchoose(NB, n - x)
 59 
    if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
- 
 
 60 
 
46 
		   - lfastchoose(N, n));
 61 
    p = ((double)n)/((double)(r+b));
- 
 
 62 
    q = ((double)(r+b-n))/((double)(r+b));
- 
 
 63 
 
- 
 
 64 
    p1 = dbinom_raw(x,r,p,q,give_log);
- 
 
 65 
    p2 = dbinom_raw(n-x,b,p,q,give_log);
- 
 
 66 
    p3 = dbinom_raw(n,r+b,p,q,give_log);
- 
 
 67 
 
- 
 
 68 
    return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
47 
}
 69 
}