bdtr.c
Go to the documentation of this file.
1 /* bdtr.c
2  *
3  * Binomial distribution
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * int k, n;
10  * double p, y, bdtr();
11  *
12  * y = bdtr( k, n, p );
13  *
14  * DESCRIPTION:
15  *
16  * Returns the sum of the terms 0 through k of the Binomial
17  * probability density:
18  *
19  * k
20  * -- ( n ) j n-j
21  * > ( ) p (1-p)
22  * -- ( j )
23  * j=0
24  *
25  * The terms are not summed directly; instead the incomplete
26  * beta integral is employed, according to the formula
27  *
28  * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
29  *
30  * The arguments must be positive, with p ranging from 0 to 1.
31  *
32  * ACCURACY:
33  *
34  * Tested at random points (a,b,p), with p between 0 and 1.
35  *
36  * a,b Relative error:
37  * arithmetic domain # trials peak rms
38  * For p between 0.001 and 1:
39  * IEEE 0,100 100000 4.3e-15 2.6e-16
40  * See also incbet.c.
41  *
42  * ERROR MESSAGES:
43  *
44  * message condition value returned
45  * bdtr domain k < 0 0.0
46  * n < k
47  * x < 0, x > 1
48  */
49 /* bdtrc()
50  *
51  * Complemented binomial distribution
52  *
53  *
54  *
55  * SYNOPSIS:
56  *
57  * int k, n;
58  * double p, y, bdtrc();
59  *
60  * y = bdtrc( k, n, p );
61  *
62  * DESCRIPTION:
63  *
64  * Returns the sum of the terms k+1 through n of the Binomial
65  * probability density:
66  *
67  * n
68  * -- ( n ) j n-j
69  * > ( ) p (1-p)
70  * -- ( j )
71  * j=k+1
72  *
73  * The terms are not summed directly; instead the incomplete
74  * beta integral is employed, according to the formula
75  *
76  * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
77  *
78  * The arguments must be positive, with p ranging from 0 to 1.
79  *
80  * ACCURACY:
81  *
82  * Tested at random points (a,b,p).
83  *
84  * a,b Relative error:
85  * arithmetic domain # trials peak rms
86  * For p between 0.001 and 1:
87  * IEEE 0,100 100000 6.7e-15 8.2e-16
88  * For p between 0 and .001:
89  * IEEE 0,100 100000 1.5e-13 2.7e-15
90  *
91  * ERROR MESSAGES:
92  *
93  * message condition value returned
94  * bdtrc domain x<0, x>1, n<k 0.0
95  */
96 /* bdtri()
97  *
98  * Inverse binomial distribution
99  *
100  *
101  *
102  * SYNOPSIS:
103  *
104  * int k, n;
105  * double p, y, bdtri();
106  *
107  * p = bdtri( k, n, y );
108  *
109  * DESCRIPTION:
110  *
111  * Finds the event probability p such that the sum of the
112  * terms 0 through k of the Binomial probability density
113  * is equal to the given cumulative probability y.
114  *
115  * This is accomplished using the inverse beta integral
116  * function and the relation
117  *
118  * 1 - p = incbi( n-k, k+1, y ).
119  *
120  * ACCURACY:
121  *
122  * Tested at random points (a,b,p).
123  *
124  * a,b Relative error:
125  * arithmetic domain # trials peak rms
126  * For p between 0.001 and 1:
127  * IEEE 0,100 100000 2.3e-14 6.4e-16
128  * IEEE 0,10000 100000 6.6e-12 1.2e-13
129  * For p between 10^-6 and 0.001:
130  * IEEE 0,100 100000 2.0e-12 1.3e-14
131  * IEEE 0,10000 100000 1.5e-12 3.2e-14
132  * See also incbi.c.
133  *
134  * ERROR MESSAGES:
135  *
136  * message condition value returned
137  * bdtri domain k < 0, n <= k 0.0
138  * x < 0, x > 1
139  */
140 
141 /* bdtr() */
142 
143 /*
144  * Cephes Math Library Release 2.3: March, 1995
145  * Copyright 1984, 1987, 1995 by Stephen L. Moshier
146  */
147 
148 #include "mconf.h"
149 
150 double bdtrc(double k, int n, double p) {
151  double dk, dn;
152  double fk = floor(k);
153 
154  if (isnan(p) || isnan(k)) {
155  return NAN;
156  }
157 
158  if (p < 0.0 || p > 1.0 || n < fk) {
159  sf_error("bdtrc", SF_ERROR_DOMAIN, NULL);
160  return NAN;
161  }
162 
163  if (fk < 0) {
164  return 1.0;
165  }
166 
167  if (fk == n) {
168  return 0.0;
169  }
170 
171  dn = n - fk;
172  if (k == 0) {
173  if (p < .01)
174  dk = -expm1(dn * log1p(-p));
175  else
176  dk = 1.0 - pow(1.0 - p, dn);
177  } else {
178  dk = fk + 1;
179  dk = incbet(dk, dn, p);
180  }
181  return dk;
182 }
183 
184 double bdtr(double k, int n, double p) {
185  double dk, dn;
186  double fk = floor(k);
187 
188  if (isnan(p) || isnan(k)) {
189  return NAN;
190  }
191 
192  if (p < 0.0 || p > 1.0 || fk < 0 || n < fk) {
193  sf_error("bdtr", SF_ERROR_DOMAIN, NULL);
194  return NAN;
195  }
196 
197  if (fk == n) return 1.0;
198 
199  dn = n - fk;
200  if (fk == 0) {
201  dk = pow(1.0 - p, dn);
202  } else {
203  dk = fk + 1.;
204  dk = incbet(dn, dk, 1.0 - p);
205  }
206  return dk;
207 }
208 
209 double bdtri(double k, int n, double y) {
210  double p, dn, dk;
211  double fk = floor(k);
212 
213  if (isnan(k)) {
214  return NAN;
215  }
216 
217  if (y < 0.0 || y > 1.0 || fk < 0.0 || n <= fk) {
218  sf_error("bdtri", SF_ERROR_DOMAIN, NULL);
219  return NAN;
220  }
221 
222  dn = n - fk;
223 
224  if (fk == n) return 1.0;
225 
226  if (fk == 0) {
227  if (y > 0.8) {
228  p = -expm1(log1p(y - 1.0) / dn);
229  } else {
230  p = 1.0 - pow(y, 1.0 / dn);
231  }
232  } else {
233  dk = fk + 1;
234  p = incbet(dn, dk, 0.5);
235  if (p > 0.5)
236  p = incbi(dk, dn, 1.0 - y);
237  else
238  p = 1.0 - incbi(dn, dk, y);
239  }
240  return p;
241 }
incbi
double incbi(double aa, double bb, double yy0)
Definition: incbi.c:51
isnan
#define isnan(X)
Definition: main.h:93
n
int n
Definition: BiCGSTAB_simple.cpp:1
bdtr
double bdtr(double k, int n, double p)
Definition: bdtr.c:184
expm1
double expm1(double x)
Definition: unity.c:106
incbet
double incbet(double aa, double bb, double xx)
Definition: incbet.c:283
log1p
double log1p(double x)
Definition: unity.c:49
y
Scalar * y
Definition: level1_cplx_impl.h:124
bdtri
double bdtri(double k, int n, double y)
Definition: bdtr.c:209
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
mconf.h
p
float * p
Definition: Tutorial_Map_using.cpp:9
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
bdtrc
double bdtrc(double k, int n, double p)
Definition: bdtr.c:150
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481


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