#include #include "helpers.h" #include #include #include "pool.h" #include "private.h" const struct fixed_string empty_string = {NULL, 0}; static void expr_error (void *aux UNUSED, const char *format, ...) { struct error e; va_list args; /* FIXME: we can do better about saying where the error occurred. */ e.class = SE; err_location (&e.where); e.title = NULL; va_start (args, format); err_vmsg (&e, format, args); va_end (args); } double expr_ymd_to_ofs (double year, double month, double day) { int y = year; int m = month; int d = day; if (y != year || m != month || d != day) { msg (SE, _("One of the arguments to a DATE function is not an integer. " "The result will be system-missing.")); return SYSMIS; } return calendar_gregorian_to_offset (y, m, d, expr_error, NULL); } double expr_ymd_to_date (double year, double month, double day) { double ofs = expr_ymd_to_ofs (year, month, day); return ofs != SYSMIS ? ofs * DAY_S : SYSMIS; } double expr_wkyr_to_date (double week, double year) { int w = week; if (w != week) { msg (SE, _("The week argument to DATE.WKYR is not an integer. " "The result will be system-missing.")); return SYSMIS; } else if (w < 1 || w > 53) { msg (SE, _("The week argument to DATE.WKYR is outside the acceptable " "range of 1 to 53. " "The result will be system-missing.")); return SYSMIS; } else { double yr_1_1 = expr_ymd_to_ofs (year, 1, 1); if (yr_1_1 != SYSMIS) return DAY_S * (yr_1_1 + WEEK_DAY * (w - 1)); else return SYSMIS; } } double expr_yrday_to_date (double year, double yday) { int yd = yday; if (yd != yday) { msg (SE, _("The day argument to DATE.YRDAY is not an integer. " "The result will be system-missing.")); return SYSMIS; } else if (yd < 1 || yd > 366) { msg (SE, _("The day argument to DATE.YRDAY is outside the acceptable " "range of 1 to 366. " "The result will be system-missing.")); return SYSMIS; } else { double yr_1_1 = expr_ymd_to_ofs (year, 1, 1); if (yr_1_1 != SYSMIS) return DAY_S * (yr_1_1 + yd - 1.); else return SYSMIS; } } double expr_yrmoda (double year, double month, double day) { if (year >= 0 && year <= 99) year += 1900; else if (year != (int) year && year > 47516) { msg (SE, _("The year argument to YRMODA is greater than 47516. " "The result will be system-missing.")); return SYSMIS; } return expr_ymd_to_ofs (year, month, day); } int compare_string (const struct fixed_string *a, const struct fixed_string *b) { size_t i; for (i = 0; i < a->length && i < b->length; i++) if (a->string[i] != b->string[i]) return a->string[i] < b->string[i] ? -1 : 1; for (; i < a->length; i++) if (a->string[i] != ' ') return 1; for (; i < b->length; i++) if (b->string[i] != ' ') return -1; return 0; } size_t count_valid (double *d, size_t d_cnt) { size_t valid_cnt; size_t i; valid_cnt = 0; for (i = 0; i < d_cnt; i++) valid_cnt += is_valid (d[i]); return valid_cnt; } struct fixed_string alloc_string (struct expression *e, size_t length) { struct fixed_string s; s.length = length; s.string = pool_alloc (e->eval_pool, length); return s; } struct fixed_string copy_string (struct expression *e, const char *old, size_t length) { struct fixed_string s = alloc_string (e, length); memcpy (s.string, old, length); return s; } struct func_params { double Ptarget; double a; double b; }; static double generalized_idf (double P, double a, double b, double (*cdf) (double, void *)) { struct func_params params; gsl_function f; gsl_root_fsolver *fsolver; int iter; int status; double root; if (P < 0 || P > 1) return SYSMIS; params.Ptarget = P; params.a = a; params.b = b; f.function = cdf; f.params = ¶ms; fsolver = gsl_root_fsolver_alloc (gsl_root_fsolver_brent); gsl_root_fsolver_set (fsolver, &f, 0, 1); iter = 0; do { double x_lower, x_upper; status = gsl_root_fsolver_iterate (fsolver); if (status != 0) return SYSMIS; x_lower = gsl_root_fsolver_x_lower (fsolver); x_upper = gsl_root_fsolver_x_upper (fsolver); status = gsl_root_test_interval (x_lower, x_upper, 0, 2 * DBL_EPSILON); } while (status == GSL_CONTINUE && ++iter < 100); root = gsl_root_fsolver_root (fsolver); gsl_root_fsolver_free (fsolver); return root; } static double cdf_beta (double x, void *params_) { struct func_params *params = params_; return gsl_cdf_beta_P (x, params->a, params->b) - params->Ptarget; } double idf_beta (double P, double a, double b) { #if 1 return generalized_idf (P, a, b, cdf_beta); #else double x = a / (a + b); double dx = 1.; while (fabs (dx) > 2 * DBL_EPSILON) { dx = (gsl_sf_beta_inc (a, b, x) - P) / gsl_ran_beta_pdf (x, a, b); x -= dx; if (x < 0) x += (dx - x) / 2; } return x; #endif } /* Returns the noncentral beta cumulative distribution function value for the given arguments. FIXME: The accuracy of this function is not entirely satisfactory. We only match the example values given in AS 310 to the first 5 significant digits. */ double ncdf_beta (double x, double a, double b, double lambda) { double c; if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.) return SYSMIS; c = lambda / 2.; if (lambda < 54.) { /* Algorithm AS 226. */ double x0, a0, beta, temp, gx, q, ax, sumq, sum; double err_max = 2 * DBL_EPSILON; double err_bound; int iter_max = 100; int iter; x0 = floor (c - 5.0 * sqrt (c)); if (x0 < 0.) x0 = 0.; a0 = a + x0; beta = (gsl_sf_lngamma (a0) + gsl_sf_lngamma (b) - gsl_sf_lngamma (a0 + b)); temp = gsl_sf_beta_inc (a0, b, x); gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0)); if (a0 >= a) q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.); else q = exp (-c); ax = q * temp; sumq = 1. - q; sum = ax; iter = 0; do { iter++; temp -= gx; gx = x * (a + b + iter - 1.) * gx / (a + iter); q *= c / iter; sumq -= q; ax = temp * q; sum += ax; err_bound = (temp - gx) * sumq; } while (iter < iter_max && err_bound > err_max); return sum; } else { /* Algorithm AS 310. */ double m, m_sqrt; int iter, iter_lower, iter_upper, iter1, iter2, j; double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s; double err_bound; double err_max = 2 * DBL_EPSILON; iter = 0; m = floor (c + .5); m_sqrt = sqrt (m); iter_lower = m - 5. * m_sqrt; iter_upper = m + 5. * m_sqrt; t = -c + m * log (c) - gsl_sf_lngamma (m + 1.); q = exp (t); r = q; psum = q; beta = (gsl_sf_lngamma (a + m) + gsl_sf_lngamma (b) - gsl_sf_lngamma (a + m + b)); s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta; fx = gx = exp (s1); ftemp = temp = gsl_sf_beta_inc (a + m, b, x); iter++; sum = q * temp; iter1 = m; while (iter1 >= iter_lower && q >= err_max) { q = q * iter1 / c; iter++; gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx; iter1--; temp += gx; psum += q; sum += q * temp; } t0 = (gsl_sf_lngamma (a + b) - gsl_sf_lngamma (a + 1.) - gsl_sf_lngamma (b)); s0 = a * log (x) + b * log (1. - x); s = 0.; for (j = 0; j < iter1; j++) { double t1; s += exp (t0 + s0 + j * log (x)); t1 = log (a + b + j) - log (a + 1. + j) + t0; t0 = t1; } err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s); q = r; temp = ftemp; gx = fx; iter2 = m; for (;;) { double ebd = err_bound + (1. - psum) * temp; if (ebd < err_max || iter >= iter_upper) break; iter2++; iter++; q = q * c / iter2; psum += q; temp -= gx; gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx; sum += q * temp; } return sum; } } double cdf_bvnor (double x0, double x1, double r) { double z = x0 * x0 - 2. * r * x0 * x1 + x1 * x1; return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r)); } double idf_fdist (double P, double df1, double df2) { double temp = idf_beta (P, df1 / 2, df2 / 2); return temp * df2 / ((1. - temp) * df1); } /* * Mathlib : A C Library of Special Functions * Copyright (C) 1998 Ross Ihaka * Copyright (C) 2000 The R Development Core Team * * This program is free software; you can redistribute it and/or * modify * it under the terms of the GNU General Public License as * published by * the Free Software Foundation; either version 2 of the * License, or * (at your option) any later version. * * This program is distributed in the hope that it will be * useful, * but WITHOUT ANY WARRANTY; without even the implied warranty * of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public * License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA * 02111-1307 USA. */ /* Returns the density of the noncentral beta distribution with noncentrality parameter LAMBDA. */ double npdf_beta (double x, double a, double b, double lambda) { if (lambda < 0. || a <= 0. || b <= 0.) return SYSMIS; else if (lambda == 0.) return gsl_ran_beta_pdf (x, a, b); else { double max_error = 2 * DBL_EPSILON; int max_iter = 200; double term = gsl_ran_beta_pdf (x, a, b); double lambda2 = 0.5 * lambda; double weight = exp (-lambda2); double sum = weight * term; double psum = weight; int k; for (k = 1; k <= max_iter && 1 - psum < max_error; k++) { weight *= lambda2 / k; term *= x * (a + b) / a; sum += weight * term; psum += weight; a += 1; } return sum; } }