b_chi {DPQ} | R Documentation |
b_chi(nu) := E[chi(nu)] / sqrt(nu) = sqrt(2/nu) Gamma((nu+1)/2) / Gamma(nu/2),
where χ(ν) denotes a chi-distributed random variable, i.e.,
the square of a chi-squared variable, and Gamma(z) is
the Gamma function, gamma()
in R.
This is a relatively important auxiliary function when computing with
non-central t distribution functions and approximations, specifically see
Johnson et al.(1994), p.520, after (31.26a), e.g., our pntJW39()
.
Its logarithm,
lb_chi(nu) := log(sqrt(2/nu) Gamma((nu+1)/2) / Gamma(nu/2)),
is even easier to compute via lgamma
and log
,
and I have used Maple to derive an asymptotic expansion in
\frac{1}{ν} as well.
Note that lb_chi(nu) also appears in the formula
for the t-density (dt
) and distribution (tail) functions.
b_chi (nu, one.minus = FALSE, c1 = 341, c2 = 1000) b_chiAsymp(nu, order = 2, one.minus = FALSE) #lb_chi (nu, ......) # not yet lb_chiAsymp(nu, order) c_dt(nu) # warning("FIXME: current c_dt() is poor -- base it on lb_chi(nu) !") c_dtAsymp(nu) # deprecated in favour of lb_chi(nu) c_pt(nu) # warning("use better c_dt()") %---> FIXME deprecate even stronger ?
nu |
non-negative numeric vector of degrees of freedom. |
one.minus |
logical indicating if 1 - b() should be returned instead of b(). |
c1, c2 |
boundaries for different approximation intervals used:
FIXME: |
order |
the polynomial order in 1/nu of the asymptotic expansion of b_χ(ν) for nu -> Inf. The default, |
One can see that b_chi()
has the properties of a CDF of a
continuous positive random variable: It grows monotonely from
b_χ(0) = 0 to (asymptotically) one. Specifically,
for large nu
, b_chi(nu) = b_chiAsymp(nu)
and
1 - b_chi(nu) ~= 1/(4 nu).
More accurately, derived from Abramowitz and Stegun, 6.1.47 (p.257) for a= 1/2, b=0,
Γ(z + 1/2) / Gamma(z) ~= sqrt(z)*(1 - 1/(8z) + 1/(128 z^2) + O(1/z^3)),
and applied for b_χ(ν) with z = ν/2, we get
b_χ(ν) \sim 1 - (1/(4ν) * (1 - 1/(8ν)) + O(ν^{-3})),
which has been implemented in b_chiAsymp(*, order=2)
in 1999.
Even more accurately, Martin Maechler, used Maple to derive an asymptotic expansion up to order 15, here reported up to order 5, namely with r := 1/(4nu),
b_chi(nu) = c_chi(r) = 1 - r + 1/2 r^2 + 5/2 r^3 - 21/8 r^4 - 399/8 r^5 + O(r^6).
a numeric vector of the same length as nu
.
Martin Maechler
Johnson, Kotz, Balakrishnan (1995)
Continuous Univariate Distributions,
Vol 2, 2nd Edition; Wiley.
Formula on page 520, after (31.26a)
Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover. https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain.
The t-distribution (base R) page pt
;
our pntJW39()
.
curve(b_chi, 0, 20); abline(h=0:1, v=0, lty=3) r <- curve(b_chi, 1e-10, 1e5, log="x") with(r, lines(x, b_chi(x, one.minus=TRUE), col = 2)) ## Zoom in to c1-region rc1 <- curve(b_chi, 340.5, 341.5, n=1001)# nothing to see e <- 1e-3; curve(b_chi, 341-e, 341+e, n=1001) # nothing e <- 1e-5; curve(b_chi, 341-e, 341+e, n=1001) # see noise, but no jump e <- 1e-7; curve(b_chi, 341-e, 341+e, n=1001) # see float "granularity"+"jump" ## Zoom in to c2-region rc2 <- curve(b_chi, 999.5, 1001.5, n=1001) # nothing visible e <- 1e-3; curve(b_chi, 1000-e, 1000+e, n=1001) # clear small jump c2 <- 1500 e <- 1e-3; curve(b_chi(x,c2=c2), c2-e, c2+e, n=1001)# still ## - - - - c2 <- 3000 e <- 1e-3; curve(b_chi(x,c2=c2), c2-e, c2+e, n=1001)# ok asymp clearly better!! curve(b_chiAsymp, add=TRUE, col=adjustcolor("red", 1/3), lwd=3) if(requireNamespace("Rmpfr")) { xm <- Rmpfr::seqMpfr(c2-e, c2+e, length.out=1000) } ## - - - - c2 <- 4000 e <- 1e-3; curve(b_chi(x,c2=c2), c2-e, c2+e, n=1001)# ok asymp clearly better!! curve(b_chiAsymp, add=TRUE, col=adjustcolor("red", 1/3), lwd=3) grCol <- adjustcolor("forest green", 1/2) curve(b_chi, 1/2, 1e11, log="x") curve(b_chiAsymp, add = TRUE, col = grCol, lwd = 3) ## 1-b(nu) ~= 1/(4 nu) a power function <==> linear in log-log scale: curve(b_chi(x, one.minus=TRUE), 1/2, 1e11, log="xy") curve(b_chiAsymp(x, one.minus=TRUE), add = TRUE, col = grCol, lwd = 3)