1 /*
2 * AUTHOR
3 * Catherine Loader, catherine@research.bell-labs.com.
4 * October 23, 2000.
5 *
6 * Merge in to R:
7 * Copyright (C) 2000, 2005 The R Core Team
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, a copy is available at
21 * https://www.R-project.org/Licenses/
22 *
23 *
24 * DESCRIPTION
25 *
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.
33 */
34
35 #include "nmath.h"
36 #include "dpq.h"
37
df(double x,double m,double n,int give_log)38 double df(double x, double m, double n, int give_log)
39 {
40 double p, q, f, dens;
41
42 #ifdef IEEE_754
43 if (ISNAN(x) || ISNAN(m) || ISNAN(n))
44 return x + m + n;
45 #endif
46 if (m <= 0 || n <= 0) ML_WARN_return_NAN;
47 if (x < 0.) return(R_D__0);
48 if (x == 0.) return(m > 2 ? R_D__0 : (m == 2 ? R_D__1 : ML_POSINF));
49 if (!R_FINITE(m) && !R_FINITE(n)) { /* both +Inf */
50 if(x == 1.) return ML_POSINF; else return R_D__0;
51 }
52 if (!R_FINITE(n)) /* must be +Inf by now */
53 return(dgamma(x, m/2, 2./m, give_log));
54 if (m > 1e14) {/* includes +Inf: code below is inaccurate there */
55 dens = dgamma(1./x, n/2, 2./n, give_log);
56 return give_log ? dens - 2*log(x): dens/(x*x);
57 }
58
59 f = 1./(n+x*m);
60 q = n*f;
61 p = x*m*f;
62
63 if (m >= 2) {
64 f = m*q/2;
65 dens = dbinom_raw((m-2)/2, (m+n-2)/2, p, q, give_log);
66 }
67 else {
68 f = m*m*q / (2*p*(m+n));
69 dens = dbinom_raw(m/2, (m+n)/2, p, q, give_log);
70 }
71 return(give_log ? log(f)+dens : f*dens);
72 }
73