The R Project SVN R

Rev

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

Rev 8431 Rev 11723
Line 1... Line 1...
1 
/*
 1 
/*
- 
 
 2 
 *  AUTHOR
2 
 *  Mathlib : A C Library of Special Functions
 3 
 *    Catherine Loader, catherine@research.bell-labs.com.
- 
 
 4 
 *    October 23, 2000.
- 
 
 5 
 *
3 
 *  Copyright (C) 1998 Ross Ihaka
 6 
 *  Merge in to R:
4 
 *  Copyright (C) 2000 The R Development Core Team
 7 
 *	Copyright (C) 2000, The R Core Development Team
5 
 *
 8 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 9 
 *  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
 10 
 *  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
 11 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 12 
 *  (at your option) any later version.
Line 13... Line 16...
13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 16 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 
 *  GNU General Public License for more details.
 17 
 *  GNU General Public License for more details.
15 
 *
 18 
 *
16 
 *  You should have received a copy of the GNU General Public License
 19 
 *  You should have received a copy of the GNU General Public License
17 
 *  along with this program; if not, write to the Free Software
 20 
 *  along with this program; if not, write to the Free Software
18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 21 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
- 
 
 22 
 *
19 
 *
 23 
 *
20 
 *  DESCRIPTION
 24 
 *  DESCRIPTION
21 
 *
 25 
 *
22 
 *    The density function of the F distribution.
 26 
 *    The density function of the F distribution.
- 
 
 27 
 *    To evaluate it, write it as a Binomial probability with p = x*m/(n+x*m).
- 
 
 28 
 *    For m >= 2, we use the simplest conversion.
- 
 
 29 
 *    For m < 2, (m-2)/2 < 0 so the conversion will not work, and we must use
- 
 
 30 
 *               a second conversion.
- 
 
 31 
 *    Note the division by p; this seems unavoidable
- 
 
 32 
 *    for m < 2, since the F density has a singularity as x (or p) -> 0.
23 
 */
 33 
 */
24 
 
 34 
 
25 
#include "nmath.h"
 35 
#include "nmath.h"
26 
#include "dpq.h"
 36 
#include "dpq.h"
27 
 
 37 
 
28 
double df(double x, double n1, double n2, int give_log)
 38 
double df(double x, double m, double n, int give_log)
29 
{
 39 
{
30 
    double a;
 40 
    double p, q, f, dens;
- 
 
 41 
 
31 
#ifdef IEEE_754
 42 
#ifdef IEEE_754
32 
    if (ISNAN(x) || ISNAN(n1) || ISNAN(n2))
 43 
    if (ISNAN(x) || ISNAN(m) || ISNAN(n))
33 
	return x + n1 + n2;
 44 
	return x + m + n;
34 
#endif
 45 
#endif
35 
    if (n1 <= 0 || n2 <= 0) ML_ERR_return_NAN;
 46 
    if (m <= 0 || n <= 0) ML_ERR_return_NAN;
- 
 
 47 
    if (x <= 0.) return(R_D__0);
- 
 
 48 
 
- 
 
 49 
    f = 1./(n+x*m);
- 
 
 50 
    q = n*f;
- 
 
 51 
    p = x*m*f;
36 
 
 52 
 
37 
    if (x <= 0.0)
 53 
    if (m >= 2) {
38 
	return R_D__0;
 54 
	f = m*q/2;
39 
    a = (n1 / n2) * x;
 55 
	dens = dbinom_raw((m-2)/2, (m+n-2)/2, p, q, give_log);
40 
    n1 /= 2;
 56 
    }
41 
    n2 /= 2;
 57 
    else {
42 
    return give_log ?
 58 
	f = m*m*q / (2*p*(m+n));
43 
	log(a) * n1  + log(1. + a)*(- (n1 + n2)) - log(x) - lbeta(n1, n2) :
 59 
	dens = dbinom_raw(m/2, (m+n)/2, p, q, give_log);
- 
 
 60 
    }
44 
	pow(a  , n1) * pow(1. + a , - (n1 + n2)) /  (  x  *  beta(n1, n2));
 61 
    return(give_log ? log(f)+dens : f*dens);
45 
}
 62 
}