LAL 7.6.1.1-4859cae
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
99static 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
117static 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
128static 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
141static 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
184static 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
196static 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
204static 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
222static 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
228static 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
237static 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
258static 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
273static 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
288static 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
311static 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
336static 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
455static 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
490static 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
557static 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
569static 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
626double 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
742double 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