The R Project SVN R

Rev

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

Rev 8628 Rev 30227
Line 26... Line 26...
26 
#include "dpq.h"
 26 
#include "dpq.h"
27 
 
 27 
 
28 
double qhyper(double p, double NR, double NB, double n,
 28 
double qhyper(double p, double NR, double NB, double n,
29 
	      int lower_tail, int log_p)
 29 
	      int lower_tail, int log_p)
30 
{
 30 
{
31 
/* This is basically the same code as  ./phyper.c -- keep in sync! */
 31 
/* This is basically the same code as  ./phyper.c  *used* to be --> FIXME! */
32 
    double N, xstart, xend, xr, xb, sum, term;
 32 
    double N, xstart, xend, xr, xb, sum, term;
33 
    int small_N;
 33 
    int small_N;
34 
 
 - 
 
35 
#ifdef IEEE_754
 34 
#ifdef IEEE_754
36 
    if (ISNAN(p) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
 35 
    if (ISNAN(p) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
37 
	return p + NR + NB + n;
 36 
	return p + NR + NB + n;
38 
#endif
 37 
#endif
39 
    if(!R_FINITE(p) || !R_FINITE(NR) || !R_FINITE(NB) || !R_FINITE(n))
 38 
    if(!R_FINITE(p) || !R_FINITE(NR) || !R_FINITE(NB) || !R_FINITE(n))
Line 42... Line 41...
42 
 
 41 
 
43 
    NR = floor(NR + 0.5);
 42 
    NR = floor(NR + 0.5);
44 
    NB = floor(NB + 0.5);
 43 
    NB = floor(NB + 0.5);
45 
    N = NR + NB;
 44 
    N = NR + NB;
46 
    n = floor(n + 0.5);
 45 
    n = floor(n + 0.5);
47 
    if (NR < 0 || NR < 0 || n < 0 || n > N)
 46 
    if (NR < 0 || NB < 0 || n < 0 || n > N)
48 
	ML_ERR_return_NAN;
 47 
	ML_ERR_return_NAN;
49 
 
 48 
 
- 
 
 49 
    /* Goal:  Find  xr (= #{red balls in sample}) such that
- 
 
 50 
     *   phyper(xr,  NR,NB, n) >= p > phyper(xr - 1,  NR,NB, n)
- 
 
 51 
     */
- 
 
 52 
 
50 
    xstart = fmax2(0, n - NB);
 53 
    xstart = fmax2(0, n - NB);
51 
    xend = fmin2(n, NR);
 54 
    xend = fmin2(n, NR);
52 
    if(p == R_DT_0) return xstart;
 55 
    if(p == R_DT_0) return xstart;
53 
    if(p == R_DT_1) return xend;
 56 
    if(p == R_DT_1) return xend;
54 
 
 57 
 
55 
    xr = xstart;
 58 
    xr = xstart;
56 
    xb = n - xr;
 59 
    xb = n - xr;/* always ( = #{black balls in sample} ) */
57 
 
 60 
 
58 
    small_N = (N < 1000); /* won't have underflow in product below */
 61 
    small_N = (N < 1000); /* won't have underflow in product below */
59 
    /* if N is small,  term := product.ratio( bin.coef );
 62 
    /* if N is small,  term := product.ratio( bin.coef );
60 
       otherwise work with its logarithm to protect against underflow */
 63 
       otherwise work with its logarithm to protect against underflow */
61 
    term = lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n);
 64 
    term = lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n);
62 
    if(small_N) term = exp(term);
 65 
    if(small_N) term = exp(term);
63 
    NR -= xr;
 66 
    NR -= xr;
64 
    NB -= xb;
 67 
    NB -= xb;
65 
 
 68 
 
66 
    sum = term;
 - 
 
67 
    if(!lower_tail || log_p) {
 69 
    if(!lower_tail || log_p) {
68 
	p = R_DT_qIv(p);
 70 
	p = R_DT_qIv(p);
69 
    }
 71 
    }
70 
    p *= 1 - 64*DBL_EPSILON;
 72 
    p *= 1 - 64*DBL_EPSILON;
71 
    sum = small_N ? term : exp(term);
 73 
    sum = small_N ? term : exp(term);