The R Project SVN R

Rev

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

Rev 5458 Rev 7671
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) 2000 The R Development Core Team
4
 *
5
 *
5
 *  This program is free software; you can redistribute it and/or modify
6
 *  This program is free software; you can redistribute it and/or modify
6
 *  it under the terms of the GNU General Public License as published by
7
 *  it under the terms of the GNU General Public License as published by
7
 *  the Free Software Foundation; either version 2 of the License, or
8
 *  the Free Software Foundation; either version 2 of the License, or
8
 *  (at your option) any later version.
9
 *  (at your option) any later version.
Line 14... Line 15...
14
 *
15
 *
15
 *  You should have received a copy of the GNU General Public License
16
 *  You should have received a copy of the GNU General Public License
16
 *  along with this program; if not, write to the Free Software
17
 *  along with this program; if not, write to the Free Software
17
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
18
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
18
 *
19
 *
19
 *  SYNOPSIS
-
 
20
 *
-
 
21
 *    #include "Mathlib.h"
-
 
22
 *    double dhyper(double x, double NR, double NB, double n);
-
 
23
 *
-
 
24
 *  DESCRIPTION
20
 *  DESCRIPTION
25
 *
21
 *
26
 *    The density of the hypergeometric distribution.
22
 *    The density of the hypergeometric distribution.
27
 */
23
 */
28
 
24
 
29
#include "Mathlib.h"
25
#include "Mathlib.h"
30
 
26
 
31
double dhyper(double x, double NR, double NB, double n)
27
double dhyper(double x, double NR, double NB, double n, int give_log)
32
{
28
{
33
    double N;
29
    double N;
34
#ifdef IEEE_754
30
#ifdef IEEE_754
35
    if (ISNAN(x) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
31
    if (ISNAN(x) || ISNAN(NR) || ISNAN(NB) || ISNAN(n))
36
	return x + NR + NB + n;
32
	return x + NR + NB + n;
Line 39... Line 35...
39
    NR = floor(NR + 0.5);
35
    NR = floor(NR + 0.5);
40
    NB = floor(NB + 0.5);
36
    NB = floor(NB + 0.5);
41
    N = NR + NB;
37
    N = NR + NB;
42
    n = floor(n + 0.5);
38
    n = floor(n + 0.5);
43
    if (NR < 0 || NB < 0 || n < 0 || n > N) {
39
    if (NR < 0 || NB < 0 || n < 0 || n > N) {
44
	ML_ERROR(ME_DOMAIN);
40
	ML_ERR_return_NAN;
45
	return ML_NAN;
-
 
46
    }
41
    }
47
    if (x < fmax2(0, n - NB) || x > fmin2(n, NR))
42
    if (x < fmax2(0, n - NB) || x > fmin2(n, NR))
48
	return 0;
43
	return R_D__0;
49
    return exp(lfastchoose(NR, x) + lfastchoose(NB, n - x)
44
    return R_D_exp(lfastchoose(NR, x) + lfastchoose(NB, n - x)
50
	       - lfastchoose(N, n));
45
		   - lfastchoose(N, n));
51
}
46
}