The asymptotic and numerical inversion of the Marcum $Q-$function

The generalized Marcum functions appear in problems of technical and scientific areas such as, for example, radar detection and communications. In mathematical statistics and probability theory these functions are called the noncentral gamma or the noncentral chi-squared cumulative distribution functions. In this paper we describe a new asymptotic method for inverting the generalized Marcum $Q-$function and for the complementary Marcum $P-$function. Also, we show how monotonicity and convexity properties of these functions can be used to find initial values for reliable Newton or secant methods to invert the function. We present details of numerical computations that show the reliability of the asymptotic approximations.


Introduction
We define the generalized Marcum Q-function by using the integral representation where μ > 0 and I μ (z) is the modified Bessel function. We also use the complementary function and the complementary relation reads There are other notations for the generalized Marcum function in the literature. Among them, probably the most popular is the following: where we have added a tilde in the definition to distinguish it from the definition we are using (1). For μ = 1 the definitions coincide with the original definition of the Marcum Q-function [1]. The relation with the notation we use is simple: and similarly for the P function. The generalized Marcum Q-function is an important function used in many applications in science and engineering, and notably in radar detection and communications. These functions also occur in statistics and probability theory, where they are called noncentral chi-squared or noncentral gamma cumulative distributions. Noncentral distributions play an important role in statistics because they arise in the power analysis of statistical tests. The central gamma cumulative distribution is in fact the incomplete gamma function, and the relation of the Marcum functions with the incomplete gamma functions is given in Section 2. For references on these application areas we refer to our recent publication [2], in which we have described reliable numerical algorithms for computing P μ (x, y) and Q μ (x, y) for a wide range of the positive real parameters μ, x, y (μ ≥ 1).
In this paper, we describe a new method for inverting the generalized Marcum functions for large values of the parameter μ. For this we use an asymptotic representation, and essential steps in the inversion algorithm are based on the inversion of the complementary error function. We have used the same approach in our recent paper [3] for numerical and asymptotic inversion algorithms for the incomplete gamma ratios. Also, we show how monotonicity and convexity properties of these functions can be used to find initial values for reliable Newton or secant methods to invert the functions. A combination of both asymptotic and numerical methods gives an efficient method of inversion for positive real values of the variables (with μ ≥ 1).

Properties of Marcum functions
Considering the Maclaurin series for the modified Bessel function and integrating term by term in (1) we obtain the series expansions P μ (x, y) = e −x ∞ n=0 x n n! P μ+n (y), x n n! Q μ+n (y).
These expansions are in terms of the incomplete gamma function ratios defined by where for μ > 0 the standard incomplete gamma functions are defined by We have the complementary relation P μ (x) + Q μ (x) = 1. Note that for the incomplete gamma function ratios P μ (y) and Q μ (y) appearing in (6) algorithms are given in [3]. The series expansion for P μ (x, y) and Q μ (x, y) given in (6) provide the standard definition for the noncentral gamma distribution in statistics: if Y is a random variable with a noncentral gamma distribution with parameter a > 0 and noncentral parameter λ ≥ 0, then the lower and upper tail probabilities for the distribution function of Y are given by [4] Here, p(n|λ) = e −λ λ n n! is the point probability of a Poisson distribution, and F(y|a), F c (y|a) are the lower and upper tail probabilities, respectively, of a central gamma distribution function with parameter a; that is, F(y|a), F c (y|a) are the incomplete gamma function ratios P a (y) and Q a (y), respectively. In statistical terminology, it is said that the noncentral gamma distribution is a mixture of central gamma distributions with Poisson weights. From the noncentral gamma, the noncentral chi-squared distribution function is derived very easily: if the random variable Y has a noncentral gamma distribution with parameters a and noncentrality parameter λ, then X = 2Y has a noncentral chi-squared distribution with parameter n = 2a and with noncentrality parameter 2λ. Particular values can be obtained from (1) and (6): and similar complementary relations for P μ (x, y). As will follow from the relations given later (and in less detail from the relations in (10)), the transition in the (x, y) quarter plane from small values of Q μ (x, y) to values close to unity occurs for large values of μ, x, y across the line y = x + μ. Above this line in the (x, y) quarter plane Q μ (x, y) is smaller than P μ (x, y) and below this line the complementary function P μ (x, y) is the smaller one. For more details on this transition line we refer to Example 1 in Section 5.2.
As we see next, this transition line is close to the inflection points of the graphs of Marcum functions both as a function of x and y.

Monotonicity and convexity properties
For the inversion process it is important to describe the monotonicity and convexity with respect to x or y, with μ fixed.
First we mention the recurrence relations and the related three-term recurrence which is also satisfied by P μ (x, y).
Taking the derivative with respect to y in (1) and using (11) we have and similarly By using the relations in (11) it follows that and we see that Q μ (x, y) (P μ (x, y)) is an increasing (decreasing) function of x and a decreasing (increasing) function of y. With respect to μ, Q μ (x, y) is increasing and P μ (x, y) is decreasing.
Regarding the convexity properties, using (12), (13), and (15) we obtain with c μ (x, y) as defined in (12). It is easy to prove, using similar ideas as in [5], that c ν (x, y) is decreasing as a function of x and increasing as a function of y. On the other hand, we have c ν (0, y) = y/ν. Therefore, if μ + 1 > y (both fixed), we have c μ+1 (x, y) < 1 for all x > 0, because c μ+1 (x, y) decreases as a function of x. Therefore For y > μ + 1 the convexity may change as a function of x (only once, because c μ+1 is monotonic). The inflection points with respect to x and y can be estimated quite sharply, as shown in [6] (by estimating the values for which c μ+1 (x, y) = 1). It can be proved that if μ > 0 then and if μ ≥ 3 2 then As we discuss later, the monotonicity and convexity properties are useful for finding initial values for the Newton or the secant methods, and they provide a reliable numerical inversion method for the Marcum functions.

Inversion of Marcum's Q-function
In statistics, the inversion of the Marcum functions refers to two different kinds of problems: (a) Computation of the quantiles of the distribution; that is, compute the p-quantile y for which F(y|a, λ) = p. (b) Computation of the noncentrality parameter of the distribution given the lower or upper tail probability; that is, for given y, a, and p obtain λ such that F(y|a, λ) = p or F c (y|a, λ) = p.
Newton algorithms for solving both problems have been proposed in the literature; see for example, [7] for problem (a) and [8]. No special study when the parameter a is large is discussed in any of these references.
To discuss the inversion of the Marcum Q-function, we follow the process described by Helstrom [9], where the inversions are linked to a specific problem in radiometry. In this reference, the inversion of the Q-function is performed in two steps; for the interpretation of these steps with respect to applications in radiometry we refer to Helstrom's paper.
In the two steps described by Helstrom we need two given numbers q 0 , q 1 , satisfying 0 < q 0 ≤ q 1 < 1. For the asymptotic inversion we assume that μ is a large parameter. The two steps are: (1) Find y from the equation and denote this value by y 0 . Recall, see (10), that Q μ (0, y) = Q μ (y) (the normalized incomplete gamma function). (2) Find x from the equation and denote this value by x 1 . The value y 0 is obtained in Step 1.

For
Step 1 we refer to [10]; see also [11], chapter 10] and [12]. For Step 2 we use an asymptotic representation of Q μ (x, y), and for these details we refer to Section 5.1. This inversion process is with respect to x with y = y 0 a fixed value obtained in Step 1.
We also consider the inversion with respect to y with fixed x; see Section 5.2.
In statistics, the inversion of Q μ (x, y) with respect to x corresponds to the problem of inverting the distribution function with respect to the noncentrality parameter given the upper tail probability. On the other hand, the inversion of P μ (x, y) with respect to y with fixed x corresponds to the problem of computing the p-quantiles of the distribution function.
Before this, we discuss the convergence of iterative methods (Newton, secant). The asymptotic methods that will be described in Section 5 will then be used to boost the convergence of iterative methods, or, when μ is large enough, they can be used as the sole method of inversion.

Inversion by iterative methods
We discuss in detail the inversion with respect to x, that is, for Step 2. For the inversion with respect to y similar results can be obtained.
First we observe that values q 1 ≤ q 0 in Step 2 do not make sense. For an explanation, consider the quadrant (x ≥ 0, y ≥ 0). When we have found a value y 0 from Step 1 in the inversion process, satisfying q 0 = Q μ (0, y 0 ) = Q μ (y 0 ), the inversion of Q μ (x, y 0 ) = q 1 can be interpreted as the inversion with respect to x in the (x ≥ 0, y ≥ 0) quadrant along the horizontal line y = y 0 . When we follow Q μ (x, y 0 ) along this line, starting at x = 0 with value q 0 , observing that this function is increasing with increasing x, it can never be equal to q 1 when q 1 < q 0 . When q 1 = q 0 the solution of the inversion is x = 0.
Because, with y 0 fixed the function Q μ (x, y 0 ) is increasing as a function of x, it changes from concave to convex. This is a favorable situation to ensure convergence of the Newton method in the second step of inversion, that is, for the computation of the value of x (called converges monotonically to x 1 for any starting valuex 0 ∈ [x + , Therefore, withx 0 equal to x + and x − we have guaranteed convergence in the cases x 1 > x + and x 1 < x − , respectively. On the other hand, in the interval (x − , x + ) the second derivative is small and the graph is close to a straight line, in which case the Newton method is safe. For this reason, any starting value in [x − , x + ] produces convergence. If y − μ − 1 < 0, we can take as starting valuex 0 = 0.
For the Newton method, it is necessary to compute the derivative of f (x), which can be written in terms of the Bessel function (15). An alternative is the secant method, which has a slightly smaller convergence rate but avoids the use of the derivative. Convergence also holds for the secant method, for similar graphical reasons as for Newton method, and the choice of initial values x − and x + is a natural and safe selection; if these values are negative, we take small positive values of x as starting values.
Depending on the value of q 1 , it may be more interesting to invert using the P-function instead of the Q-function. If q 1 is close to 1, it is more interesting to consider the value p 1 = 1 − q 1 and invert the equation P μ (x, y 0 ) = p 1 . Inversion for P can also be carried out reliably with the Newton or the secant method and the same starting values used for Q can be used for P (monotonicity an convexity properties for P immediately follow from those of Q considering (3)).
Inversion with respect to y can be done in a similar way. Q(x, y) is decreasing as a function of y and with the convexity properties of Equation (18) we conclude that choosing y − = x + μ − 3 2 and y + = x + μ − 1 is a safe selection ensuring convergence (and if these values are negative, taking instead small starting values for y).
Inversion using the secant method has been exhaustively tested and convergence holds in parameter ranges where our algorithm for the Marcum-Q function works [2]; more that 10 8 test points have been considered.

Asymptotic inversion
For the asymptotic inversion methods we use the representation Here, erfc z is the complementary error function defined by and R μ (ζ ) can be written in the form The first coefficient is The quantity ζ is given by with sign(ζ ) = sign(y − x − 1).
For details on the representation in (22) we refer to Appendix B. For the inversion process we consider the two cases: inversion with respect to x and to y.

REMARK 1. For the complementary function we have
where R μ (ζ ) and S μ (ζ ) are the same, with expansion as in (24).

REMARK 2.
For the incomplete gamma functions we have similar representations (for convenience some of the notation is as in (22)) where The function R μ (η) has the asymptotic representation The expansion is valid uniformly with respect to y ≥ 0. The coefficient C 0 (η) is regular at the transition point y = 1, and the same for all higher coefficients. For more details we refer to [13, Chapter 11].

Asymptotic inversion with respect to x
This corresponds to Step 2 described in Section 3. Throughout this case we take y = y 0 , the value obtained in Step 1, that satisfies q 0 = Q μ (0, y 0 ) = Q μ (y 0 ), and the present inversion problem is We use the method described for the incomplete gamma functions; see [3], [10], [12], and [11,Chapter 10]. We use representation (22) and start with solving the equation We call the solution ζ 0 , and considering q 1 as a function of ζ 0 , we have Considering q 1 in (31) as a function of ζ , we have (see (15)) Upon dividing, and replacing x, y with μx, μy, we obtain where Next we use an asymptotic representation of the modified Bessel function that is valid for large values of μ, uniformly with respect to ξ ≥ 0. We have 1 where η = 1 + ξ 2 + ln ξ and the coefficients U k ( p) are polynomials in p. The first two are Using the first part of (37) in (35) we obtain where which is the same as the relation in (26), and in which As in the asymptotic inversion of the incomplete gamma functions, where we obtained a similar equation as in (40), we solve this differential equation by substituting an expansion of the form where we have ζ 0 computed from Equation (32). After substituting this expansion in (40) and considering the first order of approximation for large μ, we obtain for ζ 1 the relation where f 0 (ζ ) is the part of f (ζ ) defined in (42) by putting T μ (ξ ) = 1. That is, which can be written in the form Observe that this quantity is related to the coefficient d 0 (ζ ) defined in (25). We have f 0 (ζ 0 ) = 1 − ζ 0 d 0 (ζ 0 ). In Equation (47) we may write x also with subscript 0, because it is related to ζ 0 via Equation (26) with ζ replaced by ζ 0 . We have to invert this equation to find x when ζ 0 is available. This can be done by using standard equation solvers, such as Newton's method; see Appendix A for more details.
When we have the value x corresponding to ζ 0 , we can compute f 0 (ζ 0 ) by using (47), and then ζ 1 from (45). This gives the the second-order approximation ζ ∼ ζ 0 + ζ 1 /μ. When we have higher order coefficients ζ n in the expansion of ζ , we use this ζ to find x from (26), and so on. These higher coefficients can be obtained by expanding f (ζ ) in negative powers of μ (after substitution of (44), and also by expanding the exponential function in (40)). The comparison of the coefficients of equal powers of μ gives the relations for the ζ j .
To explain a few steps, we write with f 0 (ζ ) defined in (46). Next we write Then we find where we have used

Asymptotic inversion with respect to y
In this case the inversion problem is defined by with x a given fixed value. Observe that this time q may be any value in the interval (0, 1), because for any positive x, we have have Q μ (x, 0) = 1 and Q μ (x, y) is monotonically decreasing to 0 as y → ∞; see (10) and (15). We proceed as in the previous case, computing ζ 0 from equation and by using (15). We obtain the relation (cf. (35)) We replace the Bessel function by using I μ−1 (μξ ) = I μ (μξ ) + (1/ξ )I μ (μξ ) and we use the expansion for the derivative where η and p are the same as in (38) and the coefficients V k ( p) are polynomials in p. The first two are This gives the analog of (40) in the form where in which When we take T μ (ξ ) = 1 and W μ (ξ ) = 1 we obtain g 0 (ζ 0 ), which turns out to be the same as f (ζ 0 ) given in (47). In this way we obtain for ζ the approximation ζ ∼ ζ 0 + ζ 1 /μ, where ζ 1 is given in (45). EXAMPLE 1. As an application we use q = 1 2 . Then the equation gives ζ 0 = 0 and we have where d 0 (x) is shown in the expansion (A.9), and using this in expansion (A.5) we find This is for the scaled variables in Q μ (μx, μy). For the real life variables x, y in Q μ (x, y) we have: Q μ (x, y) = 1 2 when This gives a description of the transition in the quadrant (x ≥ 0, y ≥ 0) from small values to values near unity of the Marcum functions.

REMARK 3. For applications in mathematical statistics it is of interest to consider the inversion for the P-function in the form
By using q = 1 − p we can use the inversion of the Q-function, but when p is very small the evaluation q = 1 − p does not make sense. We can repeat the analysis for the P-function, and the only change we have to make is to change the sign of ζ 0 , and to assume as earlier sign(ζ 0 ) = sign(y − x − 1). This follows from the representation in Remark 2, see (27).
In Appendix A we give more details on the inversion process.  * For the meaning of the relative errors δ 0 and δ 1 we refer to the text of Section 6.

Numerical examples
Again we consider two cases: inversion with respect to x and to y. We describe the inversion using the asymptotic methods. As we will see, we obtain good accuracy even for relatively small values of μ. If more accuracy is needed for small μ, it is always possible to do inversion using directly the secant method for the Marcum function, which can be computed using the algorithm in [2], as done in Section 4; the performance of the secant method is improved by using starting values provided by the asymptotic methods when these are accurate.
For the inversion with respect to x described in Section 5.1, we compare our method with the numerical results shown in [9]. In that paper the secant method has been used for finding the value of x when y and μ are given. Rather large values of μ and rather small values of q 0 and 1 − q 1 are considered. From the point of view of the application in radiometry considered in [9], this is due to the fact that q 0 represents a false-alarm probability (with typical low values) and q 1 represents a probability of detection (with typical large values). We use the same values of that paper q 0 = {10 −6 , 10 −8 }, q 1 = {0.9, 0.999}, and one other set of q 0 and q 1 .
In Table 1, we give the relative errors of the computations. For Step 1 of the inversion process the value y 0 satisfying Q μ (y 0 ) = q 0 is computed, and with this y 0 the value q 0 = Q μ (y 0 ) by using the algorithms for the incomplete gamma function ratios described in [3]. The displayed relative error is δ 0 = |q 0 / q 0 − 1|.
For the inversion in Step 2 we have used only the approximations obtained from the asymptotic inversion process described in Section 5.1, and not the algorithms for the Marcum function Q μ (x, y) themselves. The relative errors in Table 1 are much smaller than those given in [9], which proves the accuracy of our asymptotic inversion method in Step 2. The results for Step 1 are also better, but for that case we have used not only asymptotic methods, but highly accurate numerical algorithms from [3].
For applying Newton's method in Step 2 we have to use different starting choices.
(1) When y < 1 (see Section 1.1.1 and Figure 1) we fit a polynomial with the values f (0), f (0) and f (0). (2) When y ≥ 1 and q 1 < 1 2 (see Section 1.1.2 and Figure 2 [left and middle]) we need to compute the x 1 left of y − 1, and start at x = 0, without further fitting.
(3) When y ≥ 1 and q 1 > 1 2 (see Figure 2 [left and right]) we need the x 1 on the right of y − 1. Observe that f (x) > g(x) = x + y − √ 1 + 4x y if x > y − 1. The function g(x) has zeros at x = y − 1 and x = y + 1 (and between these values it is negative). A suitable starting value is the solution of the equation g(x) = 1 2 ζ 2 , which root is located on the right of x = y + 1.
With the values of ζ 0 , x = x 1 and y = y 0 we can compute f 0 (ζ 0 ) given in (47), and then ζ 1 of (45), either using their explicit expressions or the expansions in (A.7). For the numerical examples shown in Table 1 we have used these series expansions with terms up to k = 10 when needed.
Then we can compute ζ ∼ ζ 0 + ζ 1 /μ (see (44)), and using this value of ζ in (A.1) yields a new value x. In the described inversion process and in the computations, the x and y are scaled variables. The real life x and y are replaced with μx and μy, and the scaled values solve the equation Q μ (μx, μy) = q 1 .
For the inversion with respect to y we apply the method described in Section 5.2, which is straightforward. We show in Table 2 the results of numerical computations for several values of μ and in each case we take x = μ. We take in the approximation ζ ∼ ζ 0 + ζ 1 /μ in the columns indicated by ζ 0 , δ 0 only the term ζ 0 , and in the columns indicated by ζ 1 , δ 1 , the approximation including ζ 1 . The latter case has a better performance, as expected.
The values δ 0 , δ 1 are the corresponding relative errors |q/ q − 1|, where q is the given value and q = Q μ (x, y), with the given μ and x, and the y obtained by inversion.

Appendix A: More details on the Inversion Process
We discuss details on the inversion of Equation (26) with sign(ζ ) = sign(y − x − 1) with respect to x (with y fixed) or to y (with x fixed). The starting value ζ = ζ 0 follows from solving one of the equations (see For the inversion with respect to y we use It is also convenient to have the following expansions where f (ζ 0 ) and ζ 1 are defined in (45) and (46). The first coefficients are A.1. Details on the inversion with respect to x. We give a few details of this inversion, because different cases have to be considered. Let us denote the right-hand side of (A.1) by f (x), that is, Then and both f (x) and f (x) vanish at x = y − 1, with f (y − 1) = 1/(2y − 1), this value being well defined because x = y − 1 and y = 1 2 cannot happen, because then x = − 1 2 (a special case of the vanishing of √ 1 + 4x y ).
A.1.1. The case y < 1. For this case we refer to Figure 1. When y < 1, f (x) is monotonically increasing, starting with f (0) = y − 1 − ln(y), and the inversion of 1 2 ζ 2 0 = f (x) with respect to x can be done straightforwardly. The only point to verify is whether indeed 1 2 ζ 2 0 > f (0), otherwise there is no real positive root of the equation 1 2 ζ 2 0 = f (x). To verify this point, observe that the present ζ 0 follows from the first equation in (A.2), and because we assume that y < 1, the values q 0 , q 1 should satisfy q 1 > q 0 > 1 2 ; see Figure 1 (right). This means that the inversion of Q μ (μx, μy) = q 1 happens in the quadrant (x ≥ 0, y ≥ 0) below the line y = x + 1 (scaled variables) on the horizontal line y = y 0 , on which no transition point can be found, on which Q μ (μx, μy) is increasing, and on which ζ of representation (22) is increasing in absolute value (with ζ < 0). The starting value of this ζ (at x = 0) is the η in the representation of the incomplete gamma function in (28), which satisfies 1 2 η 2 = f (0). Hence, the ζ = ζ 0 in the relation (A.1) corresponding to the x and y values satisfying Q μ (μx, μy) = q 1 is such that 1 2 ζ 2 0 > f (x), as shown in Figure 1 (left). On the right we see the graph of Q μ (μx, μy) with no transition point (inflection point).
A.1.2. The case y > 1. When y > 1 the graph of f (x) is as in Figure 2 (left), where we see a minimum at x = y − 1. There are two possible subcases.
A.2. Details on the inversion with respect to y. In this case we solve the second equation in (A.2) with respect to y, with given x. For small values of ζ and ζ 0 we use the expansion given in (A.2) and for other values a Newton process. The function on the right-hand side of (A.1) (let us denote it by g(y)) has the derivative g (y) = y − 2x y − 1 + (y − 1) √ 1 + 4x y y 1 + √ 1 + 4x y . (A.12) It vanishes with g at y = x + 1 and g is a convex function of y. For small values of ζ we can use the series in (A.5) and for larger values Newton's method.

Appendix B: Asymptotic Representation of the Function Q µ (x, y)
The asymptotic representation Q μ (μx, μy) in Section 5 can be derived from the contour integral For P μ (μx, μy) a similar representation is valid when we take c > ρ. A slightly different form was derived in [14, section 4], and for details we refer to this paper. In (B.1) we can take c = t 0 , where t 0 the positive saddle point of φ(t), which is given by t 0 = (1 + 1 + ξ 2 )/ξ . The saddle point coalesces with the pole at ρ when y = x + 1, and we assume that 0 < t 0 < ρ. This implies y > x + 1, the domain in the (x, y)-plane where Q μ (μx, μy) ≤ P μ (μx, μy) (approximately, but true for large μ, x, and y).