The R Project SVN R

Rev

Rev 77685 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 77685 Rev 78687
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) 1999-2020  The R Core Team
4 
 *  Copyright (C) 1999-2014  The R Core Team
 4 
 *  Copyright (C) 1998       Ross Ihaka
5 
 *  Copyright (C) 2004	     Morten Welinder
 5 
 *  Copyright (C) 2004	     Morten Welinder
6 
 *  Copyright (C) 2004	     The R Foundation
 6 
 *  Copyright (C) 2004	     The R Foundation
7 
 *
 7 
 *
8 
 *  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
9 
 *  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
Line 67... Line 67...
67 
    while (x > 0 && term >= DBL_EPSILON * sum) {
 67 
    while (x > 0 && term >= DBL_EPSILON * sum) {
68 
	term *= x * (NB - n + x) / (n + 1 - x) / (NR + 1 - x);
 68 
	term *= x * (NB - n + x) / (n + 1 - x) / (NR + 1 - x);
69 
	sum += term;
 69 
	sum += term;
70 
	x--;
 70 
	x--;
71 
    }
 71 
    }
72
72 
 
73 
    double ss = (double) sum;
 73 
    double ss = (double) sum;
74 
    return log_p ? log1p(ss) : 1 + ss;
 74 
    return log_p ? log1p(ss) : 1 + ss;
75 
}
 75 
}
76 
 
 76 
 
77 
 
 77 
 
78 
/* FIXME: The old phyper() code was basically used in ./qhyper.c as well
 78 
/* FIXME: The old phyper() code was basically used in ./qhyper.c as well
79 
 * -----  We need to sync this again!
 79 
 * -----  We need to sync this again!
80 
*/
 - 
 
- 
 
 80 
                      q         m           n         k   */
81 
double phyper (double x, double NR, double NB, double n,
 81 
double phyper (double x, double NR, double NB, double n,
82 
	       int lower_tail, int log_p)
 82 
	       int lower_tail, int log_p)
83 
{
 83 
{
84 
/* Sample of  n balls from  NR red  and	 NB black ones;	 x are red */
 84 
/* Sample of  n balls from  NR red  and	 NB black ones;	 x are red */
85 
 
 85 
 
Line 105... Line 105...
105 
	NR = oldNB;
 105 
	NR = oldNB;
106 
	x = n - x - 1;
 106 
	x = n - x - 1;
107 
	lower_tail = !lower_tail;
 107 
	lower_tail = !lower_tail;
108 
    }
 108 
    }
109 
 
 109 
 
- 
 
 110 
    /* support of dhyper() as a function of its parameters
- 
 
 111 
     * R:  .suppHyper <- function(m,n,k) max(0, k-n) : min(k, m)
- 
 
 112 
     * --  where R's (m,n, k) == (NR,NB, n)  here */
110 
    if (x < 0)
 113 
    if (x < 0 || x < n - NB)
111 
	return R_DT_0;
 114 
	return R_DT_0;
112 
    if (x >= NR || x >= n)
 115 
    if (x >= NR || x >= n)
113 
	return R_DT_1;
 116 
	return R_DT_1;
114 
 
 - 
 
115 
    d  = dhyper (x, NR, NB, n, log_p);
 117 
    d  = dhyper (x, NR, NB, n, log_p);
- 
 
 118 
    // dhyper(.., log_p=FALSE) > 0 mathematically, but not always numerically :
- 
 
 119 
    if((!log_p && d == 0.) ||
- 
 
 120 
        (log_p && d == ML_NEGINF))
- 
 
 121 
	return R_DT_0;
116 
    pd = pdhyper(x, NR, NB, n, log_p);
 122 
    pd = pdhyper(x, NR, NB, n, log_p);
117 
 
 123 
 
118 
    return log_p ? R_DT_Log(d + pd) : R_D_Lval(d * pd);
 124 
    return log_p ? R_DT_Log(d + pd) : R_D_Lval(d * pd);
119 
}
 125 
}
- 
 
 126 
 
- 
 
 127 
// NB: MM has code for  AS 152 (Lund, 1980) >> R_77 (Shea, 1989) >> R_86 (Berger, 1991)