cdfkmu: Cumulative Distribution Function of the Kappa-Mu Distribution

Description

This function computes the cumulative probability or nonexceedance probability of the Kappa-Mu ($\kappa:\mu$) distribution given parameters ($\kappa$ and $\mu$) computed by parkmu. The cumulative distribution function is complex and numerical integration of the probability density function pdfkmu is used. Alternatively, the cumulative distribution function may be defined in terms of the Marcum Q function $$ F(x) = 1 - Q_\nu\biggl(\sqrt{2\kappa\mu},\, x\sqrt{2(1+\kappa)\mu}\biggr)\mbox{,} $$ where $F(x)$ is the nonexceedance probability for quantile $x$ and $Q_v(a,b)$ is the Marcum Q function defined by $$ Q_\nu(a,b) = \frac{1}{\alpha^{\nu-1}}\int_b^\infty t^\nu \, \exp(-(t^2 + a^2)/2) \, I_{\nu-1}(at)\; \mathrm{d}t\mbox{,} $$ which can be numerically difficult to work with and particularly so with real number values for $\nu$. $I_\nu(a)$ is the ``$\nu$th-order modified Bessel function of the first kind.''

Following an apparent breakthrough(?) by Shi (2012), $\nu$ can be written as $\nu = n + \Delta$ where $n$ is an integer and $0 < \Delta \le 1$. The author of lmomco refers to this alternative formulation as the “delta nu method”. The Marcum Q function for $\nu > 0$ ($n = 1,2,3, \cdots)$ is $$ Q_\nu(a,b) = Q_\Delta(a,b) + \exp(-(a^2 + b^2)/2) \, \sum_{i=0}^{n-1}\biggl(\frac{b}{a}\biggr)^{i+\Delta} \, I_{i+\Delta}(ab)\mbox{,} $$ and the function for $\nu \le 0$ ($n=-1,-2,-3,\cdots$) is $$ Q_\nu(a,b) = Q_\Delta(a,b) - \mathrm{exp}(-(a^2 + b^2)/2) \times \sum_{i=n}^{-1}\biggl(\frac{b}{a}\biggr)^{i+\Delta} \mathrm{I}_{i+\Delta}(ab)\mbox{,} $$ and the function for $\nu = 0$ is $$ Q_\nu(a,b) = Q_\Delta(a,b) + \mathrm{exp}(-(a^2 + b^2)/2)\mbox{.} $$

Shi (2012) concludes that the “merit” of these two expressions is that the evaulation of the Marcum Q function is reduced to the numerical evaluation of $Q_\Delta(a,b)$. This difference can result in measurably faster computation times (confirmed by limited analysis by the author of lmomco) and possibly better numerical performance.

Shi (2012) uses notation and text that implies evaluation of the far-right additive term (the summation) for $n=0$ as part of the condition $\nu > 0$. To clarify, Shi (2012) implies for $\nu > 0; n = 0$ (but $n=0$ occurs also for $-1 < \nu <= 0$) the following computation $$ Q_\nu(a,b) = Q_\Delta(a,b) + \mathrm{exp}(-(a^2 + b^2)/2) \times \biggl[\biggl(\frac{b}{a}\biggr)^{\Delta} \mathrm{I}_{\Delta}(ab) + \biggl(\frac{b}{a}\biggr)^{\Delta-1} \mathrm{I}_{\Delta-1}(ab)\biggr] $$ This result produces incompatible cumulative distribution functions of the distribution using $Q_\nu(a,b)$ for $-1 < \nu < 1$. Therefore, the author of lmomco concludes that Shi (2012) is in error (or your author misinterprets the summation notation) and that the specific condition for $\nu = 0$ shown above and lacking $\sum$ is correct; there are three individual and separate conditions to support the Marcum Q function using the “delta nu method”: $\nu \le -1$, $-1 < \nu < 1$, and $\nu \ge -1$.

Usage

cdfkmu(x, para, paracheck=TRUE, getmed=TRUE, qualo=NA, quahi=NA,
                marcumQ=TRUE, marcumQmethod=c("chisq", "delta", "integral"))

Value

Nonexceedance probability ($F$) for $x$.

Arguments

x: A real value vector.
para: The parameters from parkmu or vec2par.
paracheck: A logical controlling whether the parameters and checked for validity.
getmed: Numerical problems rolling onto the distribution from the right can result in erroneous $F$ being integrated of pdfkmu. This option is used to interrupt recurrsion, but if TRUE, then the median will be computed and for those $x$ values less than the median and $F$ initially computing as greater than 50 percent, are reset to 0. Users are unlikely to need this option changed. But the hack can be turned off by setting getmed=FALSE as the user level.
qualo: A lower limit of the range of $x$ to look for a uniroot of $F(x) = 0.5$ to estimate the median quantile that is used to mitigate for erroneous numerical results. This argument is passed along to quakmu but also used as a truncation point for which $F=1$ is returned if $x <$ qualo. Lastly, see the last example below.
quahi: An upper limit of the range of $x$ to look for a uniroot of $F(x) = 0.5$ to estimate the median quantile that is used to mitigate for erroneous numerical results. This argument is passed along to quakmu but also used as a truncation point for which $F=1$ is returned if $x >$ quahi. Lastly, see the last example below.
marcumQ: A logical controlling whether the Marcum Q function is used instead of numerical integration of pdfkmu.
marcumQmethod: Which method for Marcum Q computation is to be used (see source code).

Author

W.H. Asquith

References

Shi, Q., 2012, Semi-infinite Gauss-Hermite quadrature based approximations to the generalized Marcum and Nuttall Q-functions and further applications: First IEEE International Conference on Communications in China---Communications Theory and Security (CTS), pp. 268--273, ISBN 978--1--4673--2815--9,12.

Weinberg, G.V., 2006, Poisson representation and Monte Carlo estimation of generalized Marcum Q-function: IEEE Transactions on Aerospace and Electronic Systems, v. 42, no. 4, pp. 1520--1531.