1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2008, 2010, 2011, 2015, 2016 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include "math/distributions.h"
21 #include <gsl/gsl_cdf.h>
22 #include <gsl/gsl_randist.h>
23 #include <gsl/gsl_sf.h>
24 #include <libpspp/misc.h>
27 #include "data/val-type.h"
29 /* Returns the noncentral beta cumulative distribution function
30 value for the given arguments.
32 FIXME: The accuracy of this function is not entirely
33 satisfactory. We only match the example values given in AS
34 310 to the first 5 significant digits. */
36 ncdf_beta (double x, double a, double b, double lambda)
40 if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
46 /* Algorithm AS 226. */
47 double x0, a0, beta, temp, gx, q, ax, sumq, sum;
48 double err_max = 2 * DBL_EPSILON;
53 x0 = floor (c - 5.0 * sqrt (c));
57 beta = (gsl_sf_lngamma (a0)
59 - gsl_sf_lngamma (a0 + b));
60 temp = gsl_sf_beta_inc (a0, b, x);
61 gx = exp (a0 * log (x) + b * log (1. - x) - beta - log (a0));
63 q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
75 gx = x * (a + b + iter - 1.) * gx / (a + iter);
81 err_bound = (temp - gx) * sumq;
83 while (iter < iter_max && err_bound > err_max);
89 /* Algorithm AS 310. */
91 int iter, iter_lower, iter_upper, iter1, iter2, j;
92 double t, q, r, psum, beta, s1, gx, fx, temp, ftemp, t0, s0, sum, s;
94 double err_max = 2 * DBL_EPSILON;
100 iter_lower = m - 5. * m_sqrt;
101 iter_upper = m + 5. * m_sqrt;
103 t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
107 beta = (gsl_sf_lngamma (a + m)
109 - gsl_sf_lngamma (a + m + b));
110 s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
112 ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
117 while (iter1 >= iter_lower && q >= err_max)
121 gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
128 t0 = (gsl_sf_lngamma (a + b)
129 - gsl_sf_lngamma (a + 1.)
130 - gsl_sf_lngamma (b));
131 s0 = a * log (x) + b * log (1. - x);
134 for (j = 0; j < iter1; j++)
137 s += exp (t0 + s0 + j * log (x));
138 t1 = log (a + b + j) - log (a + 1. + j) + t0;
142 err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
149 double ebd = err_bound + (1. - psum) * temp;
150 if (ebd < err_max || iter >= iter_upper)
158 gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
167 cdf_bvnor (double x0, double x1, double r)
169 double z = pow2 (x0) - 2. * r * x0 * x1 + pow2 (x1);
170 return exp (-z / (2. * (1 - r * r))) * (2. * M_PI * sqrt (1 - r * r));
174 idf_fdist (double P, double df1, double df2)
176 double temp = gsl_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
177 return temp * df2 / ((1. - temp) * df1);
181 * Mathlib : A C Library of Special Functions
182 * Copyright (C) 1998 Ross Ihaka
183 * Copyright (C) 2000 The R Development Core Team
185 * This program is free software; you can redistribute it and/or
187 * it under the terms of the GNU General Public License as
189 * the Free Software Foundation; either version 2 of the
191 * (at your option) any later version.
193 * This program is distributed in the hope that it will be
195 * but WITHOUT ANY WARRANTY; without even the implied warranty
197 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
198 * GNU General Public License for more details.
200 * You should have received a copy of the GNU General Public
202 * along with this program; if not, write to the Free Software
203 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
207 /* Returns the density of the noncentral beta distribution with
208 noncentrality parameter LAMBDA. */
210 npdf_beta (double x, double a, double b, double lambda)
212 if (lambda < 0. || a <= 0. || b <= 0.)
214 else if (lambda == 0.)
215 return gsl_ran_beta_pdf (x, a, b);
218 double max_error = 2 * DBL_EPSILON;
220 double term = gsl_ran_beta_pdf (x, a, b);
221 double lambda2 = 0.5 * lambda;
222 double weight = exp (-lambda2);
223 double sum = weight * term;
224 double psum = weight;
226 for (k = 1; k <= max_iter && 1 - psum < max_error; k++) {
227 weight *= lambda2 / k;
228 term *= x * (a + b) / a;
229 sum += weight * term;