The R Project SVN R

Rev

Rev 11723 | Rev 36820 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
574 ihaka 1 
/*
11723 maechler 2 
 *  AUTHOR
3 
 *    Catherine Loader, catherine@research.bell-labs.com.
4 
 *    October 23, 2000.
574 ihaka 5 
 *
11723 maechler 6 
 *  Merge in to R:
12976 pd 7 
 *	Copyright (C) 2000, 2001 The R Core Development Team
11723 maechler 8 
 *
574 ihaka 9 
 *  This program is free software; you can redistribute it and/or modify
10 
 *  it under the terms of the GNU General Public License as published by
11 
 *  the Free Software Foundation; either version 2 of the License, or
12 
 *  (at your option) any later version.
13 
 *
14 
 *  This program is distributed in the hope that it will be useful,
15 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17 
 *  GNU General Public License for more details.
18 
 *
19 
 *  You should have received a copy of the GNU General Public License
20 
 *  along with this program; if not, write to the Free Software
5458 ripley 21 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
574 ihaka 22 
 *
23 
 *
11723 maechler 24 
 * DESCRIPTION
25 
 *
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
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.
574 ihaka 36 
 */
37 
 
8431 ripley 38 
#include "nmath.h"
8074 maechler 39 
#include "dpq.h"
574 ihaka 40 
 
11723 maechler 41 
double dhyper(double x, double r, double b, double n, int give_log)
42 
{
43 
    double p, q, p1, p2, p3;
44 
 
574 ihaka 45 
#ifdef IEEE_754
11723 maechler 46 
    if (ISNAN(x) || ISNAN(r) || ISNAN(b) || ISNAN(n))
47 
        return x + r + b + n;
574 ihaka 48 
#endif
11723 maechler 49 
 
50 
    if (R_D_notnnegint(r) || R_D_notnnegint(b) || R_D_notnnegint(n) || n > r+b)
7671 maechler 51 
	ML_ERR_return_NAN;
11723 maechler 52 
 
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);
57 
    n = R_D_forceint(n);
58 
 
12976 pd 59 
    if (n < x || r < x || n - x > b) return(R_D__0);
11723 maechler 60 
    if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
61 
 
62 
    p = ((double)n)/((double)(r+b));
63 
    q = ((double)(r+b-n))/((double)(r+b));
64 
 
12976 pd 65 
    p1 = dbinom_raw(x,  r, p,q,give_log);
66 
    p2 = dbinom_raw(n-x,b, p,q,give_log);
67 
    p3 = dbinom_raw(n,r+b, p,q,give_log);
11723 maechler 68 
 
69 
    return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
574 ihaka 70 
}