LAL  7.5.0.1-8083555
XLALMarcumQ.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2014 Kipp Cannon
3  *
4  * Implementation of the algorithm described in
5  *
6  * A. Gil, J. Segura, and N. M. Temme. Algorithm 939: Computation of the
7  * Marcum Q-Function. ACM Transactions on Mathematical Software (TOMS),
8  * Volume 40 Issue 3, April 2014, Article No. 20. arXiv:1311.0681
9  *
10  * with a few modifications. In particular, a different expression is used
11  * here for the 0th term in the asymptotic expansion for large (xy) than
12  * shown in Section 4.1 such that it remains numerically stable in the x ~=
13  * y limit, and a special case of the recursion implemented that remains
14  * valid when x = y.
15  *
16  * This program is free software; you can redistribute it and/or modify it
17  * under the terms of the GNU General Public License as published by the
18  * Free Software Foundation; either version 2 of the License, or (at your
19  * option) any later version.
20  *
21  * This program is distributed in the hope that it will be useful, but
22  * WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
24  * Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License along
27  * with this program; if not, write to the Free Software Foundation, Inc.,
28  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
29  */
30 
31 /*
32  * ============================================================================
33  *
34  * Preamble
35  *
36  * ============================================================================
37  */
38 
39 
40 /*
41  * stuff from the C library
42  */
43 
44 
45 #include <math.h>
46 
47 
48 /*
49  * stuff from GSL
50  */
51 
52 
53 #include <gsl/gsl_integration.h>
54 #include <gsl/gsl_sf_bessel.h>
55 #include <gsl/gsl_sf_erf.h>
56 #include <gsl/gsl_sf_gamma.h>
57 
58 
59 /*
60  * stuff from LAL
61  */
62 
63 
64 #include <lal/LALConstants.h>
65 #include <lal/XLALError.h>
66 #include <lal/LALMarcumQ.h>
67 
68 
69 /*
70  * ============================================================================
71  *
72  * Support Code
73  *
74  * ============================================================================
75  */
76 
77 
78 /*
79  * equation (8), normalized incomplete Gamma function,
80  *
81  * Q_{\mu}(x) = \Gamma(\mu, x) / \Gamma(\mu).
82  *
83  * we use the GSL implementations, but Gil et al.'s tests of their Marcum Q
84  * function implemention rely on their own implementation of these
85  * functions. GSL's implementations might not work as well as theirs and
86  * so the Marcum Q implemented here might not meet their accuracy claims.
87  * on the other hand, GSL might be using exactly their aglorithm, or it
88  * might be using an even better one, who knows.
89  *
90  * their algorithm is described in
91  *
92  * Gil , A., Segura , J., and Temme , N. M. 2012. Efficient and accurate
93  * algorithms for the computation and inversion of the incomplete gamma
94  * function ratios. SIAM J. Sci. Comput. 34(6), A2965–A2981.
95  * arXiv:1306.1754.
96  */
97 
98 
99 static double Q(double mu, double x)
100 {
101  return gsl_sf_gamma_inc_Q(mu, x);
102 }
103 
104 
105 /*
106  * \sqrt{x / y} * equation (14). see
107  *
108  * W. Gautschi, J. Slavik, On the Computation of Modified Bessel Function
109  * Ratios, Mathematics of Computation, Vol. 32, No. 143 (Jul., 1978), pp.
110  * 865-875.
111  *
112  * for information on using continued fractions to compute the ratio
113  * directly.
114  */
115 
116 
117 static double c_mu(double mu, double xi)
118 {
119  return gsl_sf_bessel_Inu_scaled(mu, xi) / gsl_sf_bessel_Inu_scaled(mu - 1., xi);
120 }
121 
122 
123 /*
124  * log of equation (32). computed using log \Gamma functions from GSL.
125  */
126 
127 
128 static double lnA(int n, double mu)
129 {
130  mu += 0.5;
131  return gsl_sf_lngamma(mu + n) - gsl_sf_lngamma(mu - n) - n * log(2.) - gsl_sf_lngamma(n + 1);
132 }
133 
134 
135 /*
136  * equation (84). \zeta^{2} / 2 = phi(xi) - phi(z0), with special case for
137  * y ~= x + 1 given by expansion in (85) and (86).
138  */
139 
140 
141 static double half_zeta_2(double x, double y)
142 {
143  if(fabs(y - x - 1.) < 1e-3) {
144  /*
145  * we do this a little differently than implied by what's
146  * in Gil et al. we rewrite (85) as
147  *
148  * -(2x + 1)^(3/2) * (y - x - 1) / (2x+1)^2 * \sum_{k} ...
149  * = -(2x + 1)^(3/2) * \sum_{k} c_{k} z^{k + 1}
150  *
151  * but then we're returning \zeta^{2} / 2,
152  *
153  * = (2x + 1)^3 * [\sum_{k} c[k] z^{k+1}]^2 / 2.
154  */
155  double c[] = {
156  +1., /* not used */
157  -(3. * x + 1.) / 3.,
158  +((72. * x + 42.) * x + 7.) / 36.,
159  -(((2700. * x + 2142.) * x + 657.) * x + 73.) / 540.,
160  +((((181440. * x + 177552.) * x + 76356.) * x + 15972.) * x + 1331.) / 12960.
161  };
162  double z = (y - x - 1.) / pow(2. * x + 1., 2.);
163  double z_to_k_p_1 = z; /* = z^{1} */
164  double sum = z_to_k_p_1; /* = c[0] * z^{1} */
165  int k;
166  for(k = 1; k < 5; k++) {
167  z_to_k_p_1 *= z;
168  sum += c[k] * z_to_k_p_1;
169  }
170  return pow(2. * x + 1., 3.) * sum * sum / 2.;
171  } else {
172  double root_1_plus_4xy = sqrt(1. + 4. * x * y);
173  return x + y - root_1_plus_4xy + log((1. + root_1_plus_4xy) / (2. * y));
174  }
175 }
176 
177 
178 /*
179  * equation (56). = \sqrt{2 (phi(xi) - phi(z0))} with the sign of x + 1 -
180  * y. phi(xi) - phi(z0) is given by the expression in equation (84).
181  */
182 
183 
184 static double zeta(double x, double y)
185 {
186  return copysign(sqrt(2. * half_zeta_2(x, y)), x + 1. - y);
187 }
188 
189 
190 /*
191  * equation (98). note: rho, r, and r'sin(theta) are all even functions
192  * of theta. we assume 0 <= theta <= pi
193  */
194 
195 
196 static double theta_over_sin_theta(double theta)
197 {
198  /* Taylor series = 1 + theta^2 / 6 + 7 theta^4 / 360. therefore,
199  * once theta is below 10^-4, the first two terms are sufficient */
200  return theta < 1e-4 ? theta * theta / 6. + 1. : theta / sin(theta);
201 }
202 
203 
204 static double theta_over_sin_theta_primed_sin_theta(double theta)
205 {
206  /* sin x * d/dx (x / sin x)
207  * = (sin x - x cos x) / sin x,
208  * = 1 - x / tan x
209  *
210  * Taylor series is
211  *
212  * theta^2 / 3 + theta^4 / 45 + 2 theta^6 / 945 + ...
213  *
214  * therefore once theta is below 10^-4 the first two terms are
215  * sufficient
216  */
217 
218  return theta < 1e-4 ? (theta * theta / 45. + 1. / 3.) * theta * theta : 1. - theta / tan(theta);
219 }
220 
221 
222 static double rho(double theta_over_sin_theta_, double xi)
223 {
224  return sqrt(theta_over_sin_theta_ * theta_over_sin_theta_ + xi * xi);
225 }
226 
227 
228 static double r(double theta, double y, double xi)
229 {
230  double theta_over_sin_theta_ = theta_over_sin_theta(theta);
231 
232  xi /= theta_over_sin_theta_;
233  return (1. + sqrt(1. + xi * xi)) * theta_over_sin_theta_ / (2. * y);
234 }
235 
236 
237 static double r_primed_sin_theta(double theta, double y, double xi)
238 {
239  double theta_over_sin_theta_ = theta_over_sin_theta(theta);
240 
241  xi /= theta_over_sin_theta_;
242  return (1. + 1. / sqrt(1. + xi * xi)) * theta_over_sin_theta_primed_sin_theta(theta) / (2. * y);
243 }
244 
245 
246 /*
247  * equation (96). f is an even function of theta.
248  *
249  * NOTE: there is a typo in both the final printed version and in the copy
250  * on arXiv, the middle term in the denominator should read "- 2 r(\theta)
251  * \cos(\theta)".
252  *
253  * we rewrite the denominator as (r - cos(theta))^2 + 1 - cos^2(theta) =
254  * ((r - cos(theta))^2 + sin^2(theta)
255  */
256 
257 
258 static double f(double theta, double y, double xi)
259 {
260  double r_ = r(theta, y, xi);
261  double r_minus_cos_theta = r_ - cos(theta);
262  double sin_theta = sin(theta);
263 
264  return (r_primed_sin_theta(theta, y, xi) - r_minus_cos_theta * r_) / (r_minus_cos_theta * r_minus_cos_theta + sin_theta * sin_theta);
265 }
266 
267 
268 /*
269  * equation (97). psi is an even function of theta.
270  */
271 
272 
273 static double psi(double theta, double xi)
274 {
275  double theta_over_sin_theta_ = theta_over_sin_theta(theta);
276  double rho_ = rho(theta_over_sin_theta_, xi);
277  double root_xi2_plus_1 = sqrt(1. + xi * xi);
278 
279  return cos(theta) * rho_ - root_xi2_plus_1 - log((theta_over_sin_theta_ + rho_) / (1. + root_xi2_plus_1));
280 }
281 
282 
283 /*
284  * equation (100), unscaled variables
285  */
286 
287 
288 static void f1_f2(double x, double M, double *f1, double *f2)
289 {
290  double a = x + M;
291  double b = sqrt(4. * x + 2. * M);
292  *f1 = a - b;
293  *f2 = a + b;
294 }
295 
296 
297 /*
298  * ============================================================================
299  *
300  * Implementations by Domain
301  *
302  * ============================================================================
303  */
304 
305 
306 /*
307  * series expansion described in section 3. equation (7).
308  */
309 
310 
311 static double MarcumQ_small_x(double M, double x, double y)
312 {
313  int n = 0;
314  double x_to_n_over_n_factorial = 1.;
315  double sum = 0.;
316 
317  do {
318  double term = x_to_n_over_n_factorial * Q(M + n, y);
319  sum += term;
320  if(term <= 1e-17 * sum)
321  break;
322  n++;
323  x_to_n_over_n_factorial *= x / n;
324  } while(1);
325 
326  return exp(-x) * sum;
327 }
328 
329 
330 /*
331  * asymptotic expansions in section 4.1. equation (37) for y > x, 1 -
332  * equation (39) for x < y, with modifications for x == y.
333  */
334 
335 
336 static double MarcumQ_large_xy(double M, double x, double y, double xi)
337 {
338  double root_y_minus_root_x = sqrt(y) - sqrt(x);
339 
340  /* equation (31) */
341  double sigma = root_y_minus_root_x * root_y_minus_root_x / xi;
342  double rho_ = sqrt(y / x);
343 
344  double rho_to_M_over_root_8_pi = pow(y / x, M / 2.) / sqrt(8. * LAL_PI);
345  double e_to_neg_sigma_xi_over_xi_to_n_minus_half = exp(-root_y_minus_root_x * root_y_minus_root_x) * sqrt(xi);
346 
347  /* Phi_{0}. equation (35). NOTE: that in (35) they show
348  * \sqrt{\sigma \xi} being replaced by \sqrt{y} - \sqrt{x}, but it
349  * should be |\sqrt{y} - \sqrt{x}|. the x ~= y limit is obtained
350  * from
351  *
352  * erfc(a x) / x = 1/x - 2/\sqrt{pi} [ a / (1 0!) - a^3 x^2 / (3 1!) + a^5 x^4 / (5 2!) - a^7 x^6 / (7 3!) + ...]
353  *
354  * the ratio |2nd term| / |0th term| is < 10^-17 when ax < 3 *
355  * 10^-6. using this,
356  *
357  * Phi_{0} = sqrt(pi / sigma) - 2 sqrt(xi)
358  *
359  * to the limits of double precision when |sqrt(y)-sqrt(x)| < 3e-6.
360  * we cheat a bit and call that 1e-5 because our final answer isn't
361  * going to be perfect anyway; and we compute the ratio of
362  * sqrt(pi)/sqrt(sigma), instead of sqrt(pi/sigma), to skirt
363  * numerical overflow until sigma is identically 0.
364  *
365  * the special case for \sigma = 0 is found by substituting
366  * \Phi_{0} from (35) into (36) and determining the desired value
367  * of \Phi_{1}. when this is done we find a net positive power of
368  * \sigma in the term involving \Phi_{0} and so evidently the
369  * correct value for \Phi_{1} will be obtained by setting \Phi_{0}
370  * to any finite value and allowing the factor of \sigma appearing
371  * in (36) to eliminate it naturally when \sigma = 0 (if we let
372  * \Phi_{0} diverge we end up with inf*0 = nan and the code fails).
373  */
374  double Phi = fabs(root_y_minus_root_x) < 1e-5 ? sigma == 0. ? 0. : sqrt(LAL_PI) / sqrt(sigma) - 2. * sqrt(xi) : sqrt(LAL_PI / sigma) * gsl_sf_erfc(fabs(root_y_minus_root_x));
375 
376  /* Psi_{0} from Phi_{0}. this is a special case of equation (38)
377  * that remains valid in the small \rho-1 limit (x ~= y). when n =
378  * 0, from equation (32) we see that A_{0}(\mu) = 1, and
379  * equation (38) can be written
380  *
381  * \Psi_{0} = \rho^{\mu - 1} / \sqrt{8\pi} (\rho-1) * \Phi_{0}.
382  *
383  * meanwhile, from equation (31) we have that
384  *
385  * \sqrt{\sigma} = |\rho - 1| / \sqrt{2 \rho}
386  *
387  * and so equation (35) can be written
388  *
389  * \Phi_{0} = \sqrt{2 \pi \rho} / |\rho - 1| erfc(\sqrt{y} - \sqrt{x}).
390  *
391  * combining these we find
392  *
393  * \Psi_{0} = \rho^{\mu - .5} / 2 erfc(\sqrt{y} - \sqrt{x}) * sgn(\rho - 1).
394  *
395  * NOTE: this result is only used for x < y (\rho > 1). for x > y
396  * (\rho < 1) we compute Q=1-P where from equation (39) -P is the
397  * sum of -\tilde{\Psi}_{n} = \Psi_{n} except -\tilde{\Psi}_{0}
398  * which is given by -1 * equation (40). comparing -1 * (40) to
399  * our result above when \rho < 1, we see that they are nearly
400  * identical except for the sign of the argument of erfc() being
401  * flipped (making it positive, again). therefore, by taking the
402  * absolute value of the erfc() argument we obtain an expression
403  * for the 0th term in the sum valid in both cases. the only
404  * remaining difference is that the sum needs to be initialized to
405  * +1 in the x > y case and 0 in the x < y case.
406  *
407  * if rho is exactly equal to 1, instead of relying on FPUs
408  * agreeing that 1-1=+0 and not -0 we enter the result with the
409  * correct sign by hand. the "correct" sign is decided by the
410  * initialization of the sum. we've chosen to initialize the x==y
411  * case to 0, so the correct sign for Psi_{0} is positive.
412  */
413  double Psi = rho_ == 1. ? .5 : copysign(pow(rho_, M - .5) * gsl_sf_erfc(fabs(root_y_minus_root_x)) / 2., rho_ - 1.);
414 
415  double sum = x > y ? 1. : 0.;
416  int n = 0;
417 
418  do {
419  /* equation (37) */
420  sum += Psi;
421  if(fabs(Psi) <= 1e-17 * fabs(sum))
422  break;
423 
424  /* next n. rho_to_M_over_root_8_pi is providing the factor
425  * of (-1)^n, so we flip its sign here, too. */
426  n++;
427  rho_to_M_over_root_8_pi = -rho_to_M_over_root_8_pi;
428  e_to_neg_sigma_xi_over_xi_to_n_minus_half /= xi;
429 
430  /* Phi_{n} from Phi_{n - 1}. equation (36) */
431  Phi = (e_to_neg_sigma_xi_over_xi_to_n_minus_half - sigma * Phi) / (n - .5);
432 
433  /* Psi_{n} from Phi_{n}. equation (38) */
434  double lnA_n_M_minus_1 = lnA(n, M - 1.);
435  Psi = rho_to_M_over_root_8_pi * exp(lnA_n_M_minus_1) * (1. - exp(lnA(n, M) - lnA_n_M_minus_1) / rho_) * Phi;
436  } while(1);
437 
438  /*
439  * for very small probabilities, improper cancellation can leave us
440  * with a small negative number. nudge back into the allowed range
441  */
442 
443  if(-1e-200 < sum && sum < 0.)
444  sum = 0.;
445 
446  return sum;
447 }
448 
449 
450 /*
451  * recurrence relation in equation (14).
452  */
453 
454 
455 static double MarcumQ_recurrence(double M, double x, double y, double xi)
456 {
457  /* factor in equation (14) we've omitted from our c_mu() */
458  double root_y_over_x = sqrt(y / x);
459 
460  /* where to start evaluating recurrence relation */
461  double mu = sqrt(2. * xi) - 1.;
462  /* make sure it differs from M by an integer */
463  mu = M - ceil(M - mu);
464 
465  double Qmu_minus_1 = XLALMarcumQmodified(mu - 1., x, y);
466  double Qmu = XLALMarcumQmodified(mu, x, y);
467 
468  do {
469  double cmu = root_y_over_x * c_mu(mu, xi);
470 
471  /* compute Q_{\mu + 1} from Q_{\mu} and Q_{\mu - 1} */
472  double Qmu_saved = Qmu;
473  Qmu = (1. + cmu) * Qmu - cmu * Qmu_minus_1;
474  Qmu_minus_1 = Qmu_saved;
475  mu++;
476  } while(mu < M);
477 
478  if(fabs(mu - M) / M > 1e-15)
479  XLAL_ERROR_REAL8(XLAL_ELOSS, "recursion failed to land on correct order: wanted M=%.16g, got M=%.16g", M, mu);
480 
481  return Qmu;
482 }
483 
484 
485 /*
486  * asymptotic expansion in section 4.2
487  */
488 
489 
490 static double MarcumQ_large_M(double M, double x, double y)
491 {
492  /* equation (56) */
493  double zeta_ = zeta(x, y);
494 
495  /* equation (67). if size of Psi array is changed, updated safety
496  * check inside loop below */
497  double e_to_neg_half_M_zeta_2 = exp(-M * half_zeta_2(x, y));
498  double Psi[20];
499  Psi[0] = sqrt(LAL_PI / (2. * M)) * gsl_sf_erfc(-zeta_ * sqrt(M / 2.));
500  Psi[1] = e_to_neg_half_M_zeta_2 / M;
501 
502  double sum = 0.;
503  int k = 1;
504 
505  do {
506  /* equation (71) */
507  double Bk = 0.;
508  int j;
509  for(j = 0; j <= k; j++)
510  Bk += /* FIXME: need f(j, k - j) * */ Psi[j] / pow(M, k - j);
511 
512  sum += Bk;
513  if(Bk <= 1e-17 * sum)
514  break;
515 
516  /* next k */
517  k++;
518 
519  /* equation (68). Psi_{j} from Psi_{j - 2}. note that we
520  * have all j < k, we just need to add the j = k-th term to
521  * the sequence */
522  if(k >= 20)
523  /* FIXME: allow dynamic allocation of Psi */
525  Psi[k] = (k - 1) / M * Psi[k - 2] + pow(-zeta_, k - 1) / M * e_to_neg_half_M_zeta_2;
526  } while(1);
527 
528  /* equation (72) */
529  return gsl_sf_erfc(-zeta_ * sqrt(M / 2.)) / 2. - sqrt(M / LAL_TWOPI) * sum;
530 }
531 
532 
533 /*
534  * quadrature method in section 5. the integrand is an even function of
535  * theta, so we just integrate over [0, +\pi] and double the result.
536  *
537  * NOTE: as theta approaches +\pi, \theta/sin(\theta) diverges in equation
538  * (97) causing the cos(\theta)\rho term to go to -\infty and the argument
539  * of the logarithm to go to +\infy and thus the logarithm to -\infty;
540  * this makes the integrand go to 0 in (95), and so we can back the upper
541  * bound away from +\pi a bit without affecting the result and in doing so
542  * avoid the need to trap overflow errors inside the argument of the
543  * exponential.
544  *
545  * NOTE: the expression in (95) only yields Q if \zeta(x, y) < 0,
546  * otherwise it yields -P.
547  */
548 
549 
551  double M;
552  double y;
553  double xi;
554 };
555 
556 
557 static double integrand(double theta, void *integrand_params)
558 {
559  struct integrand_params params = *(struct integrand_params *) integrand_params;
560 
561  /*
562  * equation (95)
563  */
564 
565  return exp(params.M * psi(theta, params.xi)) * f(theta, params.y, params.xi);
566 }
567 
568 
569 static double MarcumQ_quadrature(double M, double x, double y, double xi)
570 {
571  struct integrand_params params = {
572  .M = M,
573  .y = y,
574  .xi = xi
575  };
576  gsl_function integrand_ = {
577  .function = integrand,
578  .params = &params
579  };
580  gsl_integration_workspace *workspace = gsl_integration_workspace_alloc(20);
581  double integral = 0.;
582  double abserr = 0.;
583 
584  /* "limit" must = size of workspace. upper bound a bit less than
585  * +\pi (see above) */
586  gsl_integration_qag(&integrand_, 0, LAL_PI - 1./512., abserr, 1e-12, 20, GSL_INTEG_GAUSS41, workspace, &integral, &abserr);
587  gsl_integration_workspace_free(workspace);
588 
589  /* instead of dividing by 2*pi and doubling the integral, divide by
590  * pi */
591  integral *= exp(-M * half_zeta_2(x, y)) / LAL_PI;
592 
593  if(x + 1. < y ) /* if zeta < 0 */
594  return integral;
595  return 1 + integral;
596 }
597 
598 
599 /*
600  * ============================================================================
601  *
602  * Exported API
603  *
604  * ============================================================================
605  */
606 
607 
608 /**
609  * The modified form of the Marcum Q function. Used by Gil <i>et al.</i> in
610  *
611  * A. Gil, J. Segura, and N. M. Temme. Algorithm 939: Computation of the
612  * Marcum Q-Function. ACM Transactions on Mathematical Software (TOMS),
613  * Volume 40 Issue 3, April 2014, Article No. 20. arXiv:1311.0681
614  *
615  * The relationship between this function and the standard Marcum Q
616  * function is
617  *
618  * XLALMarcumQmodified(M, x, y) = XLALMarcumQ(M, sqrt(2. * x), sqrt(2. * y)).
619  *
620  * The function is defined for \f$1 \leq M\f$, \f$0 \leq x\f$, \f$0 \leq
621  * y\f$. Additionally, the implementation here becomes inaccurate when
622  * \f$M\f$, \f$x\f$, or \f$y\f$ is \f$\geq 10000\f$.
623  */
624 
625 
626 double XLALMarcumQmodified(double M, double x, double y)
627 {
628  double xi;
629  double f1, f2;
630  double Q_;
631 
632  /*
633  * check input
634  */
635 
636  if(M < 1.)
637  XLAL_ERROR_REAL8(XLAL_EDOM, "require 1 <= M: M=%.16g", M);
638  if(x < 0.)
639  XLAL_ERROR_REAL8(XLAL_EDOM, "0 <= x: x=%.16g", x);
640  if(y < 0.)
641  XLAL_ERROR_REAL8(XLAL_EDOM, "0 <= y: y=%.16g", y);
642 
643  if(M > 10000.)
644  XLAL_ERROR_REAL8(XLAL_ELOSS, "require M <= 10000: M=%.16g", M);
645  if(x > 10000.)
646  XLAL_ERROR_REAL8(XLAL_ELOSS, "require x <= 10000: x=%.16g", x);
647  if(y > 10000.)
648  XLAL_ERROR_REAL8(XLAL_ELOSS, "require y <= 10000: y=%.16g", y);
649 
650  /*
651  * compute some parameters
652  */
653 
654  xi = 2. * sqrt(x * y);
655  f1_f2(x, M, &f1, &f2);
656 
657  /*
658  * selection scheme from section 6
659  */
660 
661  if(x < 30.) {
662  /*
663  * use the series expansion in section 3
664  */
665 
666  Q_ = MarcumQ_small_x(M, x, y);
667  } else if(xi > 30. && M * M < 2. * xi) {
668  /*
669  * use the asymptotic expansion in section 4.1
670  */
671 
672  Q_ = MarcumQ_large_xy(M, x, y, xi);
673  } else if(f1 < y && y < f2 && M < 135.) {
674  /*
675  * use the recurrence relations in (14)
676  */
677 
678  Q_ = MarcumQ_recurrence(M, x, y, xi);
679  } else if(f1 < y && y < f2 && M >= 135.) {
680  /*
681  * use the asymptotic expansion in section 4.2
682  */
683 
684  /* FIXME: missing an implementation of f_{jk} needed by
685  * this code path */
686  XLAL_ERROR_REAL8(XLAL_EDOM, "not implemented: %s(%.16g, %.16g, %.16g)", __func__, M, x, y);
687 
688  Q_ = MarcumQ_large_M(M - 1., x / (M - 1.), y / (M - 1.));
689  } else {
690  /*
691  * use the integral representation in section 5
692  */
693 
694  Q_ = MarcumQ_quadrature(M, x / M, y / M, xi / M);
695  }
696 
697  /*
698  * we're generally only correct to 12 digits, and the error can
699  * result in an impossible probability. nudge back into the
700  * allowed range.
701  */
702 
703  if(1. < Q_ && Q_ < 1. + 1e-12)
704  Q_ = 1.;
705 
706  if(Q_ < 0. || Q_ > 1.)
707  XLAL_ERROR_REAL8(XLAL_ERANGE, "%s(%.17g, %.17g, %.17g) = %.17g", __func__, M, x, y, Q_);
708 
709  return Q_;
710 }
711 
712 
713 /**
714  * The function defined by J.\ Marcum,
715  * \f[
716  * Q_{M}(a, b) = \int_{b}^{\infty} x \left( \frac{x}{a} \right)^{M - 1} \exp \left( -\frac{x^{2} + a^{2}}{2} \right) I_{M - 1}(a x) \,\mathrm{d}x,
717  * \f]
718  * where \f$I_{M - 1}\f$ is the modified Bessel function of order \f$M -
719  * 1\f$.
720  *
721  * The CCDF for the random variable \f$x\f$ distributed according to the
722  * noncentral \f$\chi^{2}\f$ distribution with \f$k\f$ degrees-of-freedom
723  * and noncentrality parameter \f$\lambda\f$ is \f$Q_{k/2}(\sqrt{\lambda},
724  * \sqrt{x})\f$.
725  *
726  * The CCDF for the random variable \f$x\f$ distributed according to the
727  * Rice distribution with noncentrality parameter \f$\nu\f$ and width
728  * \f$\sigma\f$ is \f$Q_{1}(\nu/\sigma, x/\sigma)\f$.
729  *
730  * The probability that a signal that would be seen in a two-phase matched
731  * filter with |SNR| \f$\rho_{0}\f$ is seen to have matched filter |SNR|
732  * \f$\geq \rho\f$ in stationary Gaussian noise is \f$Q_{1}(\rho_{0},
733  * \rho)\f$.
734  *
735  * This function is implemented by computing the modified form used by Gil
736  * <I>et al.</I>,
737  *
738  * XLALMarcumQ(M, a, b) = XLALMarcumQmodified(M, a * a / 2., b * b / 2.).
739  */
740 
741 
742 double XLALMarcumQ(double M, double a, double b)
743 {
744  /*
745  * inverse of equation (5). note: all error checking handled
746  * inside XLALMarcumQmodified().
747  */
748 
749  return XLALMarcumQmodified(M, a * a / 2., b * b / 2.);
750 }
static double theta_over_sin_theta(double theta)
Definition: XLALMarcumQ.c:196
double XLALMarcumQmodified(double M, double x, double y)
The modified form of the Marcum Q function.
Definition: XLALMarcumQ.c:626
static void f1_f2(double x, double M, double *f1, double *f2)
Definition: XLALMarcumQ.c:288
static double MarcumQ_quadrature(double M, double x, double y, double xi)
Definition: XLALMarcumQ.c:569
static double f(double theta, double y, double xi)
Definition: XLALMarcumQ.c:258
static double r(double theta, double y, double xi)
Definition: XLALMarcumQ.c:228
static double integrand(double theta, void *integrand_params)
Definition: XLALMarcumQ.c:557
static double MarcumQ_recurrence(double M, double x, double y, double xi)
Definition: XLALMarcumQ.c:455
static double r_primed_sin_theta(double theta, double y, double xi)
Definition: XLALMarcumQ.c:237
static double MarcumQ_large_M(double M, double x, double y)
Definition: XLALMarcumQ.c:490
static double rho(double theta_over_sin_theta_, double xi)
Definition: XLALMarcumQ.c:222
static double theta_over_sin_theta_primed_sin_theta(double theta)
Definition: XLALMarcumQ.c:204
static double MarcumQ_small_x(double M, double x, double y)
Definition: XLALMarcumQ.c:311
static double c_mu(double mu, double xi)
Definition: XLALMarcumQ.c:117
static double half_zeta_2(double x, double y)
Definition: XLALMarcumQ.c:141
static double psi(double theta, double xi)
Definition: XLALMarcumQ.c:273
static double MarcumQ_large_xy(double M, double x, double y, double xi)
Definition: XLALMarcumQ.c:336
static double Q(double mu, double x)
Definition: XLALMarcumQ.c:99
double XLALMarcumQ(double M, double a, double b)
The function defined by J. Marcum,.
Definition: XLALMarcumQ.c:742
static double zeta(double x, double y)
Definition: XLALMarcumQ.c:184
static double lnA(int n, double mu)
Definition: XLALMarcumQ.c:128
#define LAL_PI
Archimedes's constant, pi.
Definition: LALConstants.h:179
#define LAL_TWOPI
2*pi is circumference of a circle divided by its radius
Definition: LALConstants.h:180
static const INT4 a
Definition: Random.c:79
#define XLAL_ERROR_REAL8(...)
Macro to invoke a failure from a XLAL routine returning a REAL8.
Definition: XLALError.h:752
@ XLAL_ERANGE
Output range error.
Definition: XLALError.h:411
@ XLAL_EMAXITER
Exceeded maximum number of iterations.
Definition: XLALError.h:455
@ XLAL_ELOSS
Loss of accuracy.
Definition: XLALError.h:459
@ XLAL_EDOM
Input domain error.
Definition: XLALError.h:410