The R Project SVN R

Rev

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

Rev 8431 Rev 30227
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
 3 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 1999-2000  The R Development Core Team
 4 
 *  Copyright (C) 1999-2000  The R Development Core Team
- 
 
 5 
 *  Copyright (C) 2004	     Morten Welinder
- 
 
 6 
 *  Copyright (C) 2004	     The R Foundation
5 
 *
 7 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 8 
 *  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
 9 
 *  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
 10 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 11 
 *  (at your option) any later version.
Line 18... Line 20...
18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
 20 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
19 
 *
 21 
 *
20 
 *  DESCRIPTION
 22 
 *  DESCRIPTION
21 
 *
 23 
 *
22 
 *	The distribution function of the hypergeometric distribution.
 24 
 *	The distribution function of the hypergeometric distribution.
23 
 */
 25 
 *
- 
 
 26 
 * Current implementation based on posting
- 
 
 27 
 * From: Morten Welinder <terra@gnome.org>
- 
 
 28 
 * Cc: R-bugs@biostat.ku.dk
- 
 
 29 
 * Subject: [Rd] phyper accuracy and efficiency (PR#6772)
- 
 
 30 
 * Date: Thu, 15 Apr 2004 18:06:37 +0200 (CEST)
- 
 
 31 
 ......
- 
 
 32 
 
- 
 
 33 
 The current version has very serious cancellation issues.  For example,
- 
 
 34 
 if you ask for a small right-tail you are likely to get total cancellation.
- 
 
 35 
 For example,  phyper(59, 150, 150, 60, FALSE, FALSE) gives 6.372680161e-14.
- 
 
 36 
 The right answer is dhyper(0, 150, 150, 60, FALSE) which is 5.111204798e-22.
- 
 
 37 
 
- 
 
 38 
 phyper is also really slow for large arguments.
- 
 
 39 
 
- 
 
 40 
 Therefore, I suggest using the code below. This is a sniplet from Gnumeric ...
- 
 
 41 
 The code isn't perfect.  In fact, if  x*(NR+NB)  is close to	n*NR,
- 
 
 42 
 then this code can take a while. Not longer than the old code, though.
- 
 
 43 
 
- 
 
 44 
 -- Thanks to Ian Smith for ideas.
- 
 
 45 
*/
- 
 
 46 
 
24 
#include "nmath.h"
 47 
#include "nmath.h"
25 
#include "dpq.h"
 48 
#include "dpq.h"
26 
 
 49 
 
- 
 
 50 
static double pdhyper (double x, double NR, double NB, double n, int log_p)
- 
 
 51 
{
- 
 
 52 
/*
- 
 
 53 
 * Calculate
- 
 
 54 
 *
- 
 
 55 
 *	    phyper (x, NR, NB, n, TRUE, FALSE)
- 
 
 56 
 *   [log]  ----------------------------------
- 
 
 57 
 *	       dhyper (x, NR, NB, n, FALSE)
- 
 
 58 
 *
- 
 
 59 
 * without actually calling phyper.  This assumes that
- 
 
 60 
 *
- 
 
 61 
 *     x * (NR + NB) <= n * NR
- 
 
 62 
 *
- 
 
 63 
 */
- 
 
 64 
    double sum = 0;
- 
 
 65 
    double term = 1;
- 
 
 66 
 
- 
 
 67 
    while (x > 0 && term >= DBL_EPSILON * sum) {
- 
 
 68 
	term *= x * (NB - n + x) / (n + 1 - x) / (NR + 1 - x);
- 
 
 69 
	sum += term;
- 
 
 70 
	x--;
- 
 
 71 
    }
- 
 
 72 
 
- 
 
 73 
    return log_p ? log1p(sum) : 1 + sum;
- 
 
 74 
}
- 
 
 75 
 
- 
 
 76 
 
- 
 
 77 
/* FIXME: The old phyper() code was basically used in ./qhyper.c as well
- 
 
 78 
 * -----  We need to sync this again!
- 
 
 79 
*/
27 
double phyper(double x, double NR, double NB, double n,
 80 
double phyper (double x, double NR, double NB, double n,
28 
	      int lower_tail, int log_p)
 81 
	       int lower_tail, int log_p)
29 
{
 82 
{
30 
/* Sample of  n balls from  NR red  and	 NB black ones;	 x are red */
 83 
/* Sample of  n balls from  NR red  and	 NB black ones;	 x are red */
31 
 
 84 
 
32 
/* basically the same code is used also in  ./qhyper.c -- keep in sync! */
 - 
 
33 
    double N, xstart, xend, xr, xb, sum, term;
 - 
 
34 
    int small_N;
 85 
    double d, pd;
- 
 
 86 
 
35 
#ifdef IEEE_754
 87 
#ifdef IEEE_754
36 
    if(ISNAN(x) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
 88 
    if(ISNAN(x) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
37 
	return x + NR + NB + n;
 89 
	return x + NR + NB + n;
38 
    if(!R_FINITE(x) || !R_FINITE(NR) || !R_FINITE(NB) || !R_FINITE(n))
 - 
 
39 
	ML_ERR_return_NAN;
 - 
 
40 
#endif
 90 
#endif
41 
 
 91 
 
42 
    x = floor(x + 1e-7);
 92 
    x = floor (x + 1e-7);
43 
    NR = floor(NR + 0.5);
 93 
    NR = R_D_forceint(NR);
44 
    NB = floor(NB + 0.5);
 94 
    NB = R_D_forceint(NB);
45 
    N = NR + NB;
 - 
 
46 
    n = floor(n + 0.5);
 95 
    n  = R_D_forceint(n);
- 
 
 96 
 
47 
    if (NR < 0 || NB < 0 || n < 0 || n > N)
 97 
    if (NR < 0 || NB < 0 || !R_FINITE(NR + NB) || n < 0 || n > NR + NB)
48 
	ML_ERR_return_NAN;
 98 
	ML_ERR_return_NAN;
49 
 
 99 
 
50 
    xstart = fmax2(0, n - NB);
 - 
 
51 
    xend = fmin2(n, NR);
 - 
 
52 
    if(x < xstart) return R_DT_0;
 - 
 
53 
    if(x >= xend)  return R_DT_1;
 100 
    if (x * (NR + NB) > n * NR) {
54 
 
 - 
 
55 
    xr = xstart;
 101 
	/* Swap tails.	*/
56 
    xb = n - xr;
 102 
	double oldNB = NB;
57 
 
 - 
 
58 
    small_N = (N < 1000); /* won't have underflow in product below */
 - 
 
59 
    /* if N is small,  term := product.ratio( bin.coef );
 - 
 
60 
       otherwise work with its logarithm to protect against underflow */
 - 
 
61 
    term = lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n);
 - 
 
62 
    if(small_N) term = exp(term);
 - 
 
63 
    NR -= xr;
 103 
	NB = NR;
64 
    NB -= xb;
 104 
	NR = oldNB;
65 
    sum = 0.0;
 105 
	x = n - x - 1;
66 
    while(xr <= x) {
 - 
 
67 
	sum += (small_N ? term : exp(term));
 106 
	lower_tail = !lower_tail;
68 
	xr++;
 - 
 
69 
	NB++;
 - 
 
70 
	if(small_N) term *= (NR / xr) * (xb / NB);
 - 
 
71 
	else	term += log((NR / xr) * (xb / NB));
 - 
 
72 
	xb--;
 - 
 
73 
	NR--;
 - 
 
74 
    }
 107 
    }
75 
    return R_DT_val(sum);
 - 
 
76 
}
 - 
 
77 
 
 108 
 
- 
 
 109 
    if (x < 0)
- 
 
 110 
	return R_DT_0;
- 
 
 111 
 
- 
 
 112 
    d  = dhyper (x, NR, NB, n, log_p);
- 
 
 113 
    pd = pdhyper(x, NR, NB, n, log_p);
- 
 
 114 
 
- 
 
 115 
    return log_p ? R_DT_Log(d + pd) : R_D_Lval(d * pd);
- 
 
 116 
}