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"))

Arguments

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).

Value

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

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.

Yacoub, M.D., 2007, The kappa-mu distribution and the eta-mu distribution: IEEE Antennas and Propagation Magazine, v. 49, no. 1, pp. 68--81.

Examples

Run this code

# NOT RUN {
x <- seq(0,3, by=0.5)
para <- vec2par(c(0.69, 0.625), type="kmu")
cdfkmu(x, para, marcumQ=TRUE, marcumQmethod="chisq")
cdfkmu(x, para, marcumQ=TRUE, marcumQmethod="delta")
cdfkmu(x, para, marcumQ=FALSE) # about 3 times slower
# }
# NOT RUN {
para <- vec2par(c(0.69, 0.625), type="kmu")
quahi <- supdist(para, delexp=.1)$support[2]
cdfkmu(quahi, para, quahi=quahi)
# }
# NOT RUN {
delx <- 0.01
x <- seq(0,3, by=delx)

plot(c(0,3), c(0,1), xlab="RHO", ylab="cdfkmu(RHO)", type="n")
para <- list(para=c(0, 0.75), type="kmu")
cdf <- cdfkmu(x, para)
lines(x, cdf, col=2, lwd=4)
para <- list(para=c(1, 0.5625), type="kmu")
cdf <- cdfkmu(x, para)
lines(x, cdf, col=3, lwd=4)

kappas <- c(0.00000001, 0.69, 1.37,  2.41, 4.45, 10.48, 28.49)
mus    <- c(0.75, 0.625,  0.5,  0.375, 0.25,  0.125, 0.05)
for(i in 1:length(kappas)) {
   kappa <- kappas[i]
   mu    <- mus[i]
   para <- list(para=c(kappa, mu), type="kmu")
   cdf <- cdfkmu(x, para)
   lines(x, cdf, col=i)
}
# }
# NOT RUN {
delx <- 0.005
x <- seq(0,3, by=delx)
nx <- 20*log10(x)
plot(c(-30,10), 10^c(-4,0), log="y", xaxs="i", yaxs="i",
     xlab="RHO", ylab="cdfkmu(RHO)", type="n")
m <- 1.25
mus <- c(0.25, 0.50, 0.75, 1, 1.25, 0)
for(mu in mus) {
   col <- 1
   kappa <- m/mu - 1 + sqrt((m/mu)*((m/mu)-1))
   para <- vec2par(c(kappa, mu), type="kmu")
   if(! is.finite(kappa)) {
      para <- vec2par(c(Inf,m), type="kmu")
      col <- 2
   }
   lines(nx, cdfkmu(x, para), col=col)
}
mtext("Yacoub (2007, figure 4)")
# }
# NOT RUN {
# The Marcum Q use for the CDF avoid numerical integration of pdfkmu(), but
# below is an example for which there is some failure that remains to be found.
para <- vec2par(c(10, 23), type="kmu")
# The following are reliable but slower as they avoid the Marcum Q function
# and use traditional numerical integration of the PDF function.
A <- cdfkmu(c(0.10, 0.35, 0.9, 1, 1.16), para, marcumQ=FALSE)
# Continuing, the first value in c() has an erroneous value for the next call.
B <- cdfkmu(c(0.10, 0.35, 0.9, 1, 1.16), para, marcumQ=TRUE)
# But this distribution is tightly peaks and well away from the origin, so in
# order to snap the erroneous value to zero, we need a successful median
# computation.  We can try again using the qualo argument to pass through to
# quakmu() like the following:
C <- cdfkmu(c(0.10, 0.35, 0.9, 1, 1.16), para, marcumQ=TRUE, qualo=0.4)
# The existance of the median for the last one also triggers a truncation of
# the CDF to 0 when negative solution results for the 0.35, although the
# negative is about -1E-14.
# }
# NOT RUN {
# Does the discipline of the signal litature just "know" about the apparent
# upper support of the Kappa-Mu being quite near or even at pi?
"simKMU" <- function() {
   km <- 10^runif(2, min=-3, max=3)
   f <- cdfkmu(pi, vec2par(km, type="kmu"))
   return(c(km, f))
}
EndStudy <- sapply(1:1000, function(i) { simKMU() } )
boxplot(EndStudy[3,])
# }

Run the code above in your browser using DataLab