pintos-os.org Git - pspp/blob - src/math/distributions.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2010, 2011, 2015, 2016 Free Software Foundation, Inc.
   3
   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.
   8
   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.
  13
  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/>. */
  16
  17 #include <config.h>
  18
  19 #include "math/distributions.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <gsl/gsl_randist.h>
  23 #include <gsl/gsl_sf.h>
  24 #include <libpspp/misc.h>
  25 #include <math.h>
  26
  27 #include "data/val-type.h"
  28
  29 /* Returns the noncentral beta cumulative distribution function
  30    value for the given arguments.
  31
  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. */
  35 double
  36 ncdf_beta (double x, double a, double b, double lambda)
  37 {
  38   double c;
  39
  40   if (x <= 0. || x >= 1. || a <= 0. || b <= 0. || lambda <= 0.)
  41     return SYSMIS;
  42
  43   c = lambda / 2.;
  44   if (lambda < 54.)
  45     {
  46       /* Algorithm AS 226. */
  47       double x0, a0, beta, temp, gx, q, ax, sumq, sum;
  48       double err_max = 2 * DBL_EPSILON;
  49       double err_bound;
  50       int iter_max = 100;
  51       int iter;
  52
  53       x0 = floor (c - 5.0 * sqrt (c));
  54       if (x0 < 0.)
  55         x0 = 0.;
  56       a0 = a + x0;
  57       beta = (gsl_sf_lngamma (a0)
  58               + gsl_sf_lngamma (b)
  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));
  62       if (a0 >= a)
  63         q = exp (-c + x0 * log (c)) - gsl_sf_lngamma (x0 + 1.);
  64       else
  65         q = exp (-c);
  66       ax = q * temp;
  67       sumq = 1. - q;
  68       sum = ax;
  69
  70       iter = 0;
  71       do
  72         {
  73           iter++;
  74           temp -= gx;
  75           gx = x * (a + b + iter - 1.) * gx / (a + iter);
  76           q *= c / iter;
  77           sumq -= q;
  78           ax = temp * q;
  79           sum += ax;
  80
  81           err_bound = (temp - gx) * sumq;
  82         }
  83       while (iter < iter_max && err_bound > err_max);
  84
  85       return sum;
  86     }
  87   else
  88     {
  89       /* Algorithm AS 310. */
  90       double m, m_sqrt;
  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;
  93       double err_bound;
  94       double err_max = 2 * DBL_EPSILON;
  95
  96       iter = 0;
  97
  98       m = floor (c + .5);
  99       m_sqrt = sqrt (m);
 100       iter_lower = m - 5. * m_sqrt;
 101       iter_upper = m + 5. * m_sqrt;
 102
 103       t = -c + m * log (c) - gsl_sf_lngamma (m + 1.);
 104       q = exp (t);
 105       r = q;
 106       psum = q;
 107       beta = (gsl_sf_lngamma (a + m)
 108               + gsl_sf_lngamma (b)
 109               - gsl_sf_lngamma (a + m + b));
 110       s1 = (a + m) * log (x) + b * log (1. - x) - log (a + m) - beta;
 111       fx = gx = exp (s1);
 112       ftemp = temp = gsl_sf_beta_inc (a + m, b, x);
 113       iter++;
 114       sum = q * temp;
 115       iter1 = m;
 116
 117       while (iter1 >= iter_lower && q >= err_max)
 118         {
 119           q = q * iter1 / c;
 120           iter++;
 121           gx = (a + iter1) / (x * (a + b + iter1 - 1.)) * gx;
 122           iter1--;
 123           temp += gx;
 124           psum += q;
 125           sum += q * temp;
 126         }
 127
 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);
 132
 133       s = 0.;
 134       for (j = 0; j < iter1; j++)
 135         {
 136           double t1;
 137           s += exp (t0 + s0 + j * log (x));
 138           t1 = log (a + b + j) - log (a + 1. + j) + t0;
 139           t0 = t1;
 140         }
 141
 142       err_bound = (1. - gsl_sf_gamma_inc_P (iter1, c)) * (temp + s);
 143       q = r;
 144       temp = ftemp;
 145       gx = fx;
 146       iter2 = m;
 147       for (;;)
 148         {
 149           double ebd = err_bound + (1. - psum) * temp;
 150           if (ebd < err_max || iter >= iter_upper)
 151             break;
 152
 153           iter2++;
 154           iter++;
 155           q = q * c / iter2;
 156           psum += q;
 157           temp -= gx;
 158           gx = x * (a + b + iter2 - 1.) / (a + iter2) * gx;
 159           sum += q * temp;
 160         }
 161
 162       return sum;
 163     }
 164 }
 165
 166 double
 167 cdf_bvnor (double x0, double x1, double r)
 168 {
 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));
 171 }
 172
 173 double
 174 idf_fdist (double P, double df1, double df2)
 175 {
 176   double temp = gsl_cdf_beta_Pinv (P, df1 / 2, df2 / 2);
 177   return temp * df2 / ((1. - temp) * df1);
 178 }
 179
 180 /*
 181  *  Mathlib : A C Library of Special Functions
 182  *  Copyright (C) 1998 Ross Ihaka
 183  *  Copyright (C) 2000 The R Development Core Team
 184  *
 185  *  This program is free software; you can redistribute it and/or
 186  *  modify
 187  *  it under the terms of the GNU General Public License as
 188  *  published by
 189  *  the Free Software Foundation; either version 2 of the
 190  *  License, or
 191  *  (at your option) any later version.
 192  *
 193  *  This program is distributed in the hope that it will be
 194  *  useful,
 195  *  but WITHOUT ANY WARRANTY; without even the implied warranty
 196  *  of
 197  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 198  *  GNU General Public License for more details.
 199  *
 200  *  You should have received a copy of the GNU General Public
 201  *  License
 202  *  along with this program; if not, write to the Free Software
 203  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 204  *  02110-1301 USA.
 205  */
 206
 207 /* Returns the density of the noncentral beta distribution with
 208    noncentrality parameter LAMBDA. */
 209 double
 210 npdf_beta (double x, double a, double b, double lambda)
 211 {
 212   if (lambda < 0. || a <= 0. || b <= 0.)
 213     return SYSMIS;
 214   else if (lambda == 0.)
 215     return gsl_ran_beta_pdf (x, a, b);
 216   else
 217     {
 218       double max_error = 2 * DBL_EPSILON;
 219       int max_iter = 200;
 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;
 225       int k;
 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;
 230         psum += weight;
 231         a += 1;
 232       }
 233       return sum;
 234     }
 235 }
 236