21 |
/* Modified for cdfs by Brian Gough, June 2003 */ |
/* Modified for cdfs by Brian Gough, June 2003 */ |
22 |
|
|
23 |
static double |
static double |
24 |
beta_cont_frac (const double a, const double b, const double x) |
beta_cont_frac (const double a, const double b, const double x, |
25 |
|
const double epsabs) |
26 |
{ |
{ |
27 |
const unsigned int max_iter = 512; /* control iterations */ |
const unsigned int max_iter = 512; /* control iterations */ |
28 |
const double cutoff = 2.0 * GSL_DBL_MIN; /* control the zero cutoff */ |
const double cutoff = 2.0 * GSL_DBL_MIN; /* control the zero cutoff */ |
34 |
double den_term = 1.0 - (a + b) * x / (a + 1.0); |
double den_term = 1.0 - (a + b) * x / (a + 1.0); |
35 |
|
|
36 |
if (fabs (den_term) < cutoff) |
if (fabs (den_term) < cutoff) |
37 |
den_term = cutoff; |
den_term = GSL_NAN; |
38 |
|
|
39 |
den_term = 1.0 / den_term; |
den_term = 1.0 / den_term; |
40 |
cf = den_term; |
cf = den_term; |
48 |
/* first step */ |
/* first step */ |
49 |
den_term = 1.0 + coeff * den_term; |
den_term = 1.0 + coeff * den_term; |
50 |
num_term = 1.0 + coeff / num_term; |
num_term = 1.0 + coeff / num_term; |
51 |
|
|
52 |
if (fabs (den_term) < cutoff) |
if (fabs (den_term) < cutoff) |
53 |
den_term = cutoff; |
den_term = GSL_NAN; |
54 |
|
|
55 |
if (fabs (num_term) < cutoff) |
if (fabs (num_term) < cutoff) |
56 |
num_term = cutoff; |
num_term = GSL_NAN; |
57 |
|
|
58 |
den_term = 1.0 / den_term; |
den_term = 1.0 / den_term; |
59 |
|
|
60 |
delta_frac = den_term * num_term; |
delta_frac = den_term * num_term; |
65 |
/* second step */ |
/* second step */ |
66 |
den_term = 1.0 + coeff * den_term; |
den_term = 1.0 + coeff * den_term; |
67 |
num_term = 1.0 + coeff / num_term; |
num_term = 1.0 + coeff / num_term; |
68 |
|
|
69 |
if (fabs (den_term) < cutoff) |
if (fabs (den_term) < cutoff) |
70 |
den_term = cutoff; |
den_term = GSL_NAN; |
71 |
|
|
72 |
if (fabs (num_term) < cutoff) |
if (fabs (num_term) < cutoff) |
73 |
num_term = cutoff; |
num_term = GSL_NAN; |
74 |
|
|
75 |
den_term = 1.0 / den_term; |
den_term = 1.0 / den_term; |
76 |
|
|
80 |
if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON) |
if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON) |
81 |
break; |
break; |
82 |
|
|
83 |
|
if (cf * fabs (delta_frac - 1.0) < epsabs) |
84 |
|
break; |
85 |
|
|
86 |
++iter_count; |
++iter_count; |
87 |
} |
} |
88 |
|
|
92 |
return cf; |
return cf; |
93 |
} |
} |
94 |
|
|
95 |
|
/* The function beta_inc_AXPY(A,Y,a,b,x) computes A * beta_inc(a,b,x) |
96 |
|
+ Y taking account of possible cancellations when using the |
97 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
hypergeometric transformation beta_inc(a,b,x)=1-beta(b,a,1-x). |
98 |
|
|
99 |
|
It also adjusts the accuracy of beta_inc() to fit the overall |
100 |
|
absolute error when A*beta_inc is added to Y. (e.g. if Y >> |
101 |
|
A*beta_inc then the accuracy of beta_inc can be reduced) */ |
102 |
|
|
103 |
static double |
static double |
104 |
beta_inc (const double a, const double b, const double x) |
beta_inc_AXPY (const double A, const double Y, |
105 |
|
const double a, const double b, const double x) |
106 |
{ |
{ |
107 |
if (x == 0.0) |
if (x == 0.0) |
108 |
{ |
{ |
109 |
return 0; |
return A * 0 + Y; |
110 |
} |
} |
111 |
else if (x == 1.0) |
else if (x == 1.0) |
112 |
{ |
{ |
113 |
return 1; |
return A * 1 + Y; |
114 |
} |
} |
115 |
else |
else |
116 |
{ |
{ |
117 |
double ln_beta = gsl_sf_lnbeta (a, b); |
double ln_beta = gsl_sf_lnbeta (a, b); |
118 |
double ln_pre = -ln_beta + a * log (x) + b * log1p (-x); |
double ln_pre = -ln_beta + a * log (x) + b * log1p (-x); |
119 |
|
|
120 |
double prefactor = exp(ln_pre); |
double prefactor = exp (ln_pre); |
121 |
|
|
122 |
if (x < (a + 1.0) / (a + b + 2.0)) |
if (x < (a + 1.0) / (a + b + 2.0)) |
123 |
{ |
{ |
124 |
/* Apply continued fraction directly. */ |
/* Apply continued fraction directly. */ |
125 |
double cf = beta_cont_frac (a, b, x); |
double epsabs = fabs (Y / (A * prefactor / a)) * GSL_DBL_EPSILON; |
126 |
|
|
127 |
return prefactor * cf / a; |
double cf = beta_cont_frac (a, b, x, epsabs); |
128 |
|
|
129 |
|
return A * (prefactor * cf / a) + Y; |
130 |
} |
} |
131 |
else |
else |
132 |
{ |
{ |
133 |
/* Apply continued fraction after hypergeometric transformation. */ |
/* Apply continued fraction after hypergeometric transformation. */ |
134 |
double cf = beta_cont_frac (b, a, 1.0 - x); |
double epsabs = |
135 |
|
fabs ((A + Y) / (A * prefactor / b)) * GSL_DBL_EPSILON; |
136 |
|
double cf = beta_cont_frac (b, a, 1.0 - x, epsabs); |
137 |
double term = prefactor * cf / b; |
double term = prefactor * cf / b; |
138 |
|
|
139 |
return 1 - term; |
if (A == -Y) |
140 |
|
{ |
141 |
|
return -A * term; |
142 |
|
} |
143 |
|
else |
144 |
|
{ |
145 |
|
return A * (1 - term) + Y; |
146 |
|
} |
147 |
} |
} |
148 |
} |
} |
149 |
} |
} |
150 |
|
|
151 |
|
/* Direct series evaluation for testing purposes only */ |
152 |
|
|
153 |
|
#if 0 |
154 |
|
static double |
155 |
|
beta_series (const double a, const double b, const double x, |
156 |
|
const double epsabs) |
157 |
|
{ |
158 |
|
double f = x / (1 - x); |
159 |
|
double c = (b - 1) / (a + 1) * f; |
160 |
|
double s = 1; |
161 |
|
double n = 0; |
162 |
|
|
163 |
|
s += c; |
164 |
|
|
165 |
|
do |
166 |
|
{ |
167 |
|
n++; |
168 |
|
c *= -f * (2 + n - b) / (2 + n + a); |
169 |
|
s += c; |
170 |
|
} |
171 |
|
while (n < 512 && fabs (c) > GSL_DBL_EPSILON * fabs (s) + epsabs); |
172 |
|
|
173 |
|
s /= (1 - x); |
174 |
|
|
175 |
|
return s; |
176 |
|
} |
177 |
|
#endif |
178 |
|
|