LAL  7.1.7.1-37cf38b
XLALChisq.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005,2007,2008,2010,2015 Kipp Cannon
3  *
4  * This program is free software; you can redistribute it and/or modify it
5  * under the terms of the GNU General Public License as published by the
6  * Free Software Foundation; either version 2 of the License, or (at your
7  * option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful, but
10  * WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
12  * Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License along
15  * with this program; if not, write to the Free Software Foundation, Inc.,
16  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 
20 #include <math.h>
21 
22 
23 #include <gsl/gsl_cdf.h>
24 #include <gsl/gsl_sf.h>
25 
26 
27 #include <lal/LALChisq.h>
28 #include <lal/LALConstants.h>
29 #include <lal/XLALGSL.h>
30 #include <lal/XLALError.h>
31 
32 
33 /**
34  * Compute the natural logarithm of the complementary cumulative
35  * probability function of the \f$\chi^{2}\f$ distribution.
36  *
37  * Returns the natural logarithm of the probability that \f$x_{1}^{2} +
38  * \cdots + x_{\mathrm{dof}}^{2} \geq \chi^{2}\f$, where \f$x_{1}, \ldots,
39  * x_{\mathrm{dof}}\f$ are independent zero mean unit variance Gaussian
40  * random variables. The integral expression is \f$\ln Q = \ln
41  * \int_{\chi^{2}/2}^{\infty} x^{\frac{n}{2}-1} \mathrm{e}^{-x} /
42  * \Gamma(n/2) \, \mathrm{d}x\f$, where \f$n = \mathrm{dof} =\f$ number of
43  * degrees of freedom.
44  *
45  * Results agree with Mathematica to 15 digits or more where the two have
46  * been compared and where Mathematica is able to compute a result. For
47  * example, it does not seem to be possible to obtain a numerical result
48  * from Mathematica for the equivalent of XLALLogChisqCCDF(10000.0, 8.5)
49  * using any number of digits in the intermediate calculations, though the
50  * two implementations agree to better than 15 digits at
51  * XLALLogChisqCCDF(10000.0, 8.0)
52  */
53 
55  double chi2,
56  double dof
57 )
58 {
59  /*
60  * Notes:
61  *
62  * Q(chi^2, dof) = Gamma(dof/2, chi^2/2) / Gamma(dof/2)
63  *
64  * The incomplete gamma function, Gamma(a, x), can by approximated by
65  * computing the first few terms of
66  *
67  * Gamma(a, x) = e^-x x^(a-1) [ 1 + (a-1)/x + (a-1)(a-2)/x^2 + ... ]
68  *
69  * See Abramowitz and Stegun, (6.5.32). The natural logarithm of
70  * this is
71  *
72  * ln Gamma(a, x) = (a - 1) * ln(x) - x + ln[ 1 + (a-1)/x +
73  * (a-1)(a-2)/x^2 + ... ]
74  *
75  * and we can use the log1p() function to evaluate the log of the
76  * sum.
77  */
78 
79  double ln_prob;
80 
81  if((chi2 < 0.0) || (dof <= 0.0))
83 
84  /* start with a technique for intermediate probabilities */
85  XLAL_CALLGSL(ln_prob = log(gsl_cdf_chisq_Q(chi2, dof)));
86 
87  if(ln_prob > -LAL_LN2) {
88  /* the result turned out to be large, so recompute using
89  * the CDF */
90  XLAL_CALLGSL(ln_prob = log1p(-gsl_cdf_chisq_P(chi2, dof)));
91  } else if(ln_prob < -600.) {
92  /* the result turned out to be very small, so recompute
93  * using an asymptotic approximation */
94  const double a = dof / 2.;
95  const double x = chi2 / 2.;
96  double sum;
97  double term;
98  double i;
99 
100  /* calculate the log of the numerator. if this aymptotic
101  * form doesn't work out, we can choose to report that as
102  * an error or we can fall back to whatever result we got
103  * from GSL. this is left as a compile-time choice. if
104  * the choice is to report this as an error, it can either
105  * be reported as the failure of the series to converge or
106  * as a general loss of precision in the result (because
107  * the only option left is to fall back to the GSL result).
108  * this is also left as a compile-time choice. */
109  for(i = 1., sum = 0., term = 1.; fabs(term/sum) > 1e-17; i++) {
110  double factor = (a - i) / x;
111 #if 0
112 #if 1
113  if(fabs(factor) >= 1.)
115  if(i > 100.)
117 #else
118  if(fabs(factor) >= 1. || i > 100.)
120 #endif
121 #else
122  if(fabs(factor) >= 1. || i > 100.)
123  goto done;
124 #endif
125  term *= factor;
126  sum += term;
127  }
128  ln_prob = (a - 1.) * log(x) - x + log1p(sum);
129 
130  /* subtract the log of the denominator */
131  XLAL_CALLGSL(ln_prob -= gsl_sf_lngamma(a));
132  }
133 
134  /* check that the final answer is the log of a legal probability */
135 done:
136  if(ln_prob > 0.0)
138 
139  return ln_prob;
140 }
#define XLAL_CALLGSL(statement)
Definition: XLALGSL.h:33
#define LAL_LN2
natural log of 2, ln(2)
Definition: LALConstants.h:172
Series is diverging.
Definition: XLALError.h:457
Output range error.
Definition: XLALError.h:412
static const INT4 a
Definition: Random.c:79
double XLALLogChisqCCDF(double chi2, double dof)
Compute the natural logarithm of the complementary cumulative probability function of the distributi...
Definition: XLALChisq.c:54
Input domain error.
Definition: XLALError.h:411
#define XLAL_ERROR_REAL8(...)
Macro to invoke a failure from a XLAL routine returning a REAL8.
Definition: XLALError.h:753
Exceeded maximum number of iterations.
Definition: XLALError.h:456
Loss of accuracy.
Definition: XLALError.h:460