This function computes the Barnes Extended Hypergeometric function, which in lmomco is useful in applications involving expectations of order statistics for the Generalized Exponential Poisson (GEP) distribution (see lmomgep). The function is
$$
F_{p,q}(\bm{\mathrm{n}},\bm{\mathrm{d}}; \lambda) = \sum_{k=0}^\infty \frac{\lambda^k}{\Gamma(k+1)}\frac{\Pi_{i=1}^{p} \Gamma(n_i + k)\Gamma^{-1}{(n_i)}}{\Pi_{i=1}^{q} \Gamma(d_i + k)\Gamma^{-1}{(d_i)}}\mbox{,}
$$
where \(\bm{\mathrm{n}} = [n_1, n_2, \ldots, n_p]\) for \(p\) operands and \(\bm{\mathrm{d}} = [d_1, d_2, \ldots, d_q]\) for \(q\) operands, and \(\lambda > 0\) is a parameter.
BEhypergeo(p,q, N,D, lambda, eps=1E-12, maxit=500)An R
list is returned.
The value for the function.
The number of iterations \(j\).
The error of convergence.
An integer value.
An integer value.
A scalar or vector associated with the \(p\) summation (see Details).
A scalar or vector associated with the \(q\) summation (see Details).
A real value \(\lambda > 0\).
The relative convergence error on which to break an infinite loop.
The maximum number of interations before a mandatory break on the loop, and a warning will be issued.
W.H. Asquith
For the GEP both \(\bm{\mathrm{n}}\) and \(\bm{\mathrm{d}}\) are vectors of the same value, such as \(\bm{\mathrm{n}} = [1, \ldots, 1]\) and \(\bm{\mathrm{d}} = [2, \ldots, 2]\). This implementation is built around this need by the GEP and if the length of either vector is not equal to the operand then the first value of the vector is repeated the operand times. For example for \(\bm{\mathrm{n}}\), if n = 1, then n = rep(n[1], length(p)) and so on for \(\bm{\mathrm{d}}\). Given that n and d are vectorized for the GEP, then a shorthand is used for the GEP mathematics shown herein:
$$F^{12}_{22}(h(j+1)) \equiv F_{2,2}([1,\ldots,1], [2,\ldots,2]; h(j+1))\mbox{,}$$
for the \(h\) parameter of the distribution.
Lastly, for lmomco and the GEP the arguments only involve \(p = q = 2\) and \(N = 1\), \(D = 2\), so the function is uniquely a function of the \(h\) parameter of the distribution:
H <- 10^seq(-10,10, by=0.01)
F22 <- sapply(1:length(H), function(i) BEhypergeo(2,2,1,1, H[i])$value
plot(log10(H), log10(F22), type="l")
For this example, the solution increasingly wobbles towards large \(h\), which is further explored by
plot(log10(H[1:(length(H)-1)]), diff(log10(F22)), type="l", xlim=c(0,7))
plot(log10(H[H > 75 & H < 140]), c(NA,diff(log10(F22[H > 75 & H < 140]))),
type="b"); lines(c(2.11,2.11), c(0,10))
It can be provisionally concluded that the solution to \(F^{12}_{22}(\cdot)\) begins to be suddenly questionable because of numerical difficulties beyond \(\log(h) = 2.11\). Therefore, it is given that \(h < 128\) might be an operational numerical upper limit.
Kus, C., 2007, A new lifetime distribution: Computational Statistics and Data Analysis, v. 51, pp. 4497--4509.
lmomgep