3 #include <gsl/gsl_roots.h>
4 #include <gsl/gsl_sf.h>
5 #include <libpspp/pool.h>
8 const struct substring empty_string = {NULL, 0};
11 expr_error (void *aux UNUSED, const char *format, ...)
16 /* FIXME: we can do better about saying where the error
18 m.category = MSG_SYNTAX;
19 m.severity = MSG_ERROR;
20 msg_location (&m.where);
21 va_start (args, format);
22 m.text = xvasprintf (format, args);
29 expr_ymd_to_ofs (double year, double month, double day)
35 if (y != year || m != month || d != day)
37 msg (SE, _("One of the arguments to a DATE function is not an integer. "
38 "The result will be system-missing."));
42 return calendar_gregorian_to_offset (y, m, d, expr_error, NULL);
46 expr_ymd_to_date (double year, double month, double day)
48 double ofs = expr_ymd_to_ofs (year, month, day);
49 return ofs != SYSMIS ? ofs * DAY_S : SYSMIS;
53 expr_wkyr_to_date (double week, double year)
59 msg (SE, _("The week argument to DATE.WKYR is not an integer. "
60 "The result will be system-missing."));
63 else if (w < 1 || w > 53)
65 msg (SE, _("The week argument to DATE.WKYR is outside the acceptable "
67 "The result will be system-missing."));
72 double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
74 return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1));
81 expr_yrday_to_date (double year, double yday)
87 msg (SE, _("The day argument to DATE.YRDAY is not an integer. "
88 "The result will be system-missing."));
91 else if (yd < 1 || yd > 366)
93 msg (SE, _("The day argument to DATE.YRDAY is outside the acceptable "
95 "The result will be system-missing."));
100 double yr_1_1 = expr_ymd_to_ofs (year, 1, 1);
101 if (yr_1_1 != SYSMIS)
102 return DAY_S * (yr_1_1 + yd - 1.);
109 expr_yrmoda (double year, double month, double day)
111 if (year >= 0 && year <= 99)
113 else if (year != (int) year && year > 47516)
115 msg (SE, _("The year argument to YRMODA is greater than 47516. "
116 "The result will be system-missing."));
120 return expr_ymd_to_ofs (year, month, day);
124 compare_string (const struct substring *a, const struct substring *b)
128 for (i = 0; i < a->length && i < b->length; i++)
129 if (a->string[i] != b->string[i])
130 return a->string[i] < b->string[i] ? -1 : 1;
131 for (; i < a->length; i++)
132 if (a->string[i] != ' ')
134 for (; i < b->length; i++)
135 if (b->string[i] != ' ')
141 count_valid (double *d, size_t d_cnt)
147 for (i = 0; i < d_cnt; i++)
148 valid_cnt += is_valid (d[i]);
153 alloc_string (struct expression *e, size_t length)
157 s.string = pool_alloc (e->eval_pool, length);
162 copy_string (struct expression *e, const char *old, size_t length)
164 struct substring s = alloc_string (e, length);
165 memcpy (s.string, old, length);
169 /* Returns the noncentral beta cumulative distribution function
170 value for the given arguments.
172 FIXME: The accuracy of this function is not entirely
173 satisfactory. We only match the example values given in AS
174 310 to the first 5 significant digits. */
176 ncdf_beta (double x, double a, double b, double lambda)
180 if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
186 /* Algorithm AS 226. */
187 double x0, a0, beta, temp, gx, q, ax, sumq, sum;
188 double err_max = 2 * DBL_EPSILON;
193 x0 = floor (c - 5.0 * sqrt (c));
197 beta = (gsl_sf_lngamma (a0)
199 - gsl_sf_lngamma (a0 + b));
200 temp = gsl_sf_beta_inc (a0, b, x);
201 gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
203 q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
215 gx = x * (a + b + iter - 1.) * gx / (a + iter);
221 err_bound = (temp - gx) * sumq;
223 while (iter < iter_max && err_bound > err_max);
229 /* Algorithm AS 310. */
231 int iter, iter_lower, iter_upper, iter1, iter2, j;
232 double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
234 double err_max = 2 * DBL_EPSILON;
240 iter_lower = m - 5. * m_sqrt;
241 iter_upper = m + 5. * m_sqrt;
243 t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
247 beta = (gsl_sf_lngamma (a + m)
249 - gsl_sf_lngamma (a + m + b));
250 s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
252 ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
257 while (iter1 >= iter_lower && q >= err_max)
261 gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
268 t0 = (gsl_sf_lngamma (a + b)
269 - gsl_sf_lngamma (a + 1.)
270 - gsl_sf_lngamma (b));
271 s0 = a * log (x) + b * log (1. - x);
274 for (j = 0; j < iter1; j++)
277 s += exp (t0 + s0 + j * log (x));
278 t1 = log (a + b + j) - log (a + 1. + j) + t0;
282 err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
289 double ebd = err_bound + (1. - psum) * temp;
290 if (ebd < err_max || iter >= iter_upper)
298 gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
307 cdf_bvnor (double x0, double x1, double r)
309 double z = x0 * x0 - 2. * r * x0 * x1 + x1 * x1;
310 return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
314 idf_fdist (double P, double df1, double df2)
316 double temp = gslextras_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
317 return temp * df2 / ((1. - temp) * df1);
321 * Mathlib : A C Library of Special Functions
322 * Copyright (C) 1998 Ross Ihaka
323 * Copyright (C) 2000 The R Development Core Team
325 * This program is free software; you can redistribute it and/or
327 * it under the terms of the GNU General Public License as
329 * the Free Software Foundation; either version 2 of the
331 * (at your option) any later version.
333 * This program is distributed in the hope that it will be
335 * but WITHOUT ANY WARRANTY; without even the implied warranty
337 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
338 * GNU General Public License for more details.
340 * You should have received a copy of the GNU General Public
342 * along with this program; if not, write to the Free Software
343 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
347 /* Returns the density of the noncentral beta distribution with
348 noncentrality parameter LAMBDA. */
350 npdf_beta (double x, double a, double b, double lambda)
352 if (lambda < 0. || a <= 0. || b <= 0.)
354 else if (lambda == 0.)
355 return gsl_ran_beta_pdf (x, a, b);
358 double max_error = 2 * DBL_EPSILON;
360 double term = gsl_ran_beta_pdf (x, a, b);
361 double lambda2 = 0.5 * lambda;
362 double weight = exp (-lambda2);
363 double sum = weight * term;
364 double psum = weight;
366 for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
367 weight *= lambda2 / k;
368 term *= x * (a + b) / a;
369 sum += weight * term;