The R Project SVN R

Rev

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

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 1 
/*
- 
 
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 1998 Ross Ihaka
- 
 
 4 
 *  Copyright (C) 1999-2000  The R Development Core Team
- 
 
 5 
 *  Copyright (C) 2004	     Morten Welinder
- 
 
 6 
 *  Copyright (C) 2004	     The R Foundation
- 
 
 7 
 *
- 
 
 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
- 
 
 10 
 *  the Free Software Foundation; either version 2 of the License, or
- 
 
 11 
 *  (at your option) any later version.
- 
 
 12 
 *
- 
 
 13 
 *  This program is distributed in the hope that it will be useful,
- 
 
 14 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- 
 
 15 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- 
 
 16 
 *  GNU General Public License for more details.
- 
 
 17 
 *
- 
 
 18 
 *  You should have received a copy of the GNU General Public License
- 
 
 19 
 *  along with this program; if not, a copy is available at
- 
 
 20 
 *  http://www.r-project.org/Licenses/
- 
 
 21 
 *
- 
 
 22 
 *  DESCRIPTION
- 
 
 23 
 *
- 
 
 24 
 *	The distribution function of the hypergeometric distribution.
- 
 
 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 
 
- 
 
 47 
#include "nmath.h"
- 
 
 48 
#include "dpq.h"
- 
 
 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 
*/
- 
 
 80 
double phyper (double x, double NR, double NB, double n,
- 
 
 81 
	       int lower_tail, int log_p)
- 
 
 82 
{
- 
 
 83 
/* Sample of  n balls from  NR red  and	 NB black ones;	 x are red */
- 
 
 84 
 
- 
 
 85 
    double d, pd;
- 
 
 86 
 
- 
 
 87 
#ifdef IEEE_754
- 
 
 88 
    if(ISNAN(x) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
- 
 
 89 
	return x + NR + NB + n;
- 
 
 90 
#endif
- 
 
 91 
 
- 
 
 92 
    x = floor (x + 1e-7);
- 
 
 93 
    NR = R_D_forceint(NR);
- 
 
 94 
    NB = R_D_forceint(NB);
- 
 
 95 
    n  = R_D_forceint(n);
- 
 
 96 
 
- 
 
 97 
    if (NR < 0 || NB < 0 || !R_FINITE(NR + NB) || n < 0 || n > NR + NB)
- 
 
 98 
	ML_ERR_return_NAN;
- 
 
 99 
 
- 
 
 100 
    if (x * (NR + NB) > n * NR) {
- 
 
 101 
	/* Swap tails.	*/
- 
 
 102 
	double oldNB = NB;
- 
 
 103 
	NB = NR;
- 
 
 104 
	NR = oldNB;
- 
 
 105 
	x = n - x - 1;
- 
 
 106 
	lower_tail = !lower_tail;
- 
 
 107 
    }
- 
 
 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 
}