fdtr.c
Go to the documentation of this file.
1 /* fdtr.c
2  *
3  * F distribution
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double df1, df2;
10  * double x, y, fdtr();
11  *
12  * y = fdtr( df1, df2, x );
13  *
14  * DESCRIPTION:
15  *
16  * Returns the area from zero to x under the F density
17  * function (also known as Snedcor's density or the
18  * variance ratio density). This is the density
19  * of x = (u1/df1)/(u2/df2), where u1 and u2 are random
20  * variables having Chi square distributions with df1
21  * and df2 degrees of freedom, respectively.
22  *
23  * The incomplete beta integral is used, according to the
24  * formula
25  *
26  * P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
27  *
28  *
29  * The arguments a and b are greater than zero, and x is
30  * nonnegative.
31  *
32  * ACCURACY:
33  *
34  * Tested at random points (a,b,x).
35  *
36  * x a,b Relative error:
37  * arithmetic domain domain # trials peak rms
38  * IEEE 0,1 0,100 100000 9.8e-15 1.7e-15
39  * IEEE 1,5 0,100 100000 6.5e-15 3.5e-16
40  * IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12
41  * IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13
42  * See also incbet.c.
43  *
44  *
45  * ERROR MESSAGES:
46  *
47  * message condition value returned
48  * fdtr domain a<0, b<0, x<0 0.0
49  *
50  */
51 
52 /* fdtrc()
53  *
54  * Complemented F distribution
55  *
56  *
57  *
58  * SYNOPSIS:
59  *
60  * double df1, df2;
61  * double x, y, fdtrc();
62  *
63  * y = fdtrc( df1, df2, x );
64  *
65  * DESCRIPTION:
66  *
67  * Returns the area from x to infinity under the F density
68  * function (also known as Snedcor's density or the
69  * variance ratio density).
70  *
71  *
72  * inf.
73  * -
74  * 1 | | a-1 b-1
75  * 1-P(x) = ------ | t (1-t) dt
76  * B(a,b) | |
77  * -
78  * x
79  *
80  *
81  * The incomplete beta integral is used, according to the
82  * formula
83  *
84  * P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
85  *
86  *
87  * ACCURACY:
88  *
89  * Tested at random points (a,b,x) in the indicated intervals.
90  * x a,b Relative error:
91  * arithmetic domain domain # trials peak rms
92  * IEEE 0,1 1,100 100000 3.7e-14 5.9e-16
93  * IEEE 1,5 1,100 100000 8.0e-15 1.6e-15
94  * IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13
95  * IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12
96  * See also incbet.c.
97  *
98  * ERROR MESSAGES:
99  *
100  * message condition value returned
101  * fdtrc domain a<0, b<0, x<0 0.0
102  *
103  */
104 
105 /* fdtri()
106  *
107  * Inverse of F distribution
108  *
109  *
110  *
111  * SYNOPSIS:
112  *
113  * double df1, df2;
114  * double x, p, fdtri();
115  *
116  * x = fdtri( df1, df2, p );
117  *
118  * DESCRIPTION:
119  *
120  * Finds the F density argument x such that the integral
121  * from -infinity to x of the F density is equal to the
122  * given probability p.
123  *
124  * This is accomplished using the inverse beta integral
125  * function and the relations
126  *
127  * z = incbi( df2/2, df1/2, p )
128  * x = df2 (1-z) / (df1 z).
129  *
130  * Note: the following relations hold for the inverse of
131  * the uncomplemented F distribution:
132  *
133  * z = incbi( df1/2, df2/2, p )
134  * x = df2 z / (df1 (1-z)).
135  *
136  * ACCURACY:
137  *
138  * Tested at random points (a,b,p).
139  *
140  * a,b Relative error:
141  * arithmetic domain # trials peak rms
142  * For p between .001 and 1:
143  * IEEE 1,100 100000 8.3e-15 4.7e-16
144  * IEEE 1,10000 100000 2.1e-11 1.4e-13
145  * For p between 10^-6 and 10^-3:
146  * IEEE 1,100 50000 1.3e-12 8.4e-15
147  * IEEE 1,10000 50000 3.0e-12 4.8e-14
148  * See also fdtrc.c.
149  *
150  * ERROR MESSAGES:
151  *
152  * message condition value returned
153  * fdtri domain p <= 0 or p > 1 NaN
154  * v < 1
155  *
156  */
157 
158 /*
159  * Cephes Math Library Release 2.3: March, 1995
160  * Copyright 1984, 1987, 1995 by Stephen L. Moshier
161  */
162 
163 
164 #include "mconf.h"
165 
166 
167 double fdtrc(double a, double b, double x)
168 {
169  double w;
170 
171  if ((a <= 0.0) || (b <= 0.0) || (x < 0.0)) {
172  sf_error("fdtrc", SF_ERROR_DOMAIN, NULL);
173  return NAN;
174  }
175  w = b / (b + a * x);
176  return incbet(0.5 * b, 0.5 * a, w);
177 }
178 
179 
180 double fdtr(double a, double b, double x)
181 {
182  double w;
183 
184  if ((a <= 0.0) || (b <= 0.0) || (x < 0.0)) {
185  sf_error("fdtr", SF_ERROR_DOMAIN, NULL);
186  return NAN;
187  }
188  w = a * x;
189  w = w / (b + w);
190  return incbet(0.5 * a, 0.5 * b, w);
191 }
192 
193 
194 double fdtri(double a, double b, double y)
195 {
196  double w, x;
197 
198  if ((a <= 0.0) || (b <= 0.0) || (y <= 0.0) || (y > 1.0)) {
199  sf_error("fdtri", SF_ERROR_DOMAIN, NULL);
200  return NAN;
201  }
202  y = 1.0 - y;
203  /* Compute probability for x = 0.5. */
204  w = incbet(0.5 * b, 0.5 * a, 0.5);
205  /* If that is greater than y, then the solution w < .5.
206  * Otherwise, solve at 1-y to remove cancellation in (b - b*w). */
207  if (w > y || y < 0.001) {
208  w = incbi(0.5 * b, 0.5 * a, y);
209  x = (b - b * w) / (a * w);
210  }
211  else {
212  w = incbi(0.5 * a, 0.5 * b, 1.0 - y);
213  x = b * w / (a * (1.0 - w));
214  }
215  return x;
216 }
incbi
double incbi(double aa, double bb, double yy0)
Definition: incbi.c:51
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
b
Scalar * b
Definition: benchVecAdd.cpp:17
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
fdtrc
double fdtrc(double a, double b, double x)
Definition: fdtr.c:167
incbet
double incbet(double aa, double bb, double xx)
Definition: incbet.c:283
fdtri
double fdtri(double a, double b, double y)
Definition: fdtr.c:194
y
Scalar * y
Definition: level1_cplx_impl.h:124
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
fdtr
double fdtr(double a, double b, double x)
Definition: fdtr.c:180
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
NULL
#define NULL
Definition: ccolamd.c:609
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:32