b_chi {DPQ} | R Documentation |
Compute E[\chi_\nu] / \sqrt{\nu}
useful for t- and chi-Distributions
Description
b_\chi(\nu) := E[\chi(\nu)] / \sqrt{\nu} = \frac{\sqrt{2/\nu}\Gamma((\nu+1)/2)}{\Gamma(\nu/2)},
where \chi(\nu)
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\bigl(\frac{\sqrt{2/\nu}\Gamma((\nu+1)/2)}{\Gamma(\nu/2)}\bigr),
is even easier to compute via lgamma
and log
,
and I have used Maple to derive an asymptotic expansion in
\frac{1}{\nu}
as well.
Note that lb_\chi(\nu)
also appears in the formula
for the t-density (dt
) and distribution (tail) functions.
Usage
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 ?
Arguments
nu |
non-negative numeric vector of degrees of freedom. |
one.minus |
logical indicating if |
c1 , c2 |
boundaries for different approximation intervals used:
FIXME: |
order |
the polynomial order in The default, |
Details
One can see that b_chi()
has the properties of a CDF of a
continuous positive random variable: It grows monotonely from
b_\chi(0) = 0
to (asymptotically) one. Specifically,
for large nu
, b_chi(nu) = b_chiAsymp(nu)
and
1 - b_\chi(\nu) \sim \frac{1}{4\nu}.
More accurately, derived from Abramowitz and Stegun, 6.1.47 (p.257) for a= 1/2, b=0,
\Gamma(z + 1/2) / \Gamma(z) \sim \sqrt(z)*(1 - 1/(8z) + 1/(128 z^2) + O(1/z^3)),
and applied for b_\chi(\nu)
with z = \nu/2
, we get
b_\chi(\nu) \sim 1 - (1/(4\nu) * (1 - 1/(8\nu)) + O(\nu^{-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 := \frac{1}{4\nu}
,
b_\chi(\nu) = c_\chi(r) = 1 - r + \frac{1}{2}r^2 +
\frac{5}{2}r^3 - \frac{21}{8}r^4 - \frac{399}{8}r^5 + O(r^6).
Value
a numeric vector of the same length as nu
.
Author(s)
Martin Maechler
References
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.
See Also
The t-distribution (base R) page pt
;
our pntJW39()
.
Examples
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)