I want to find the support of a non-central chi-squared distribution ($99.9 \%$ of the energy). For example, If I have a Gaussian distribution with parameters $\mu$ and $\sigma$, I know $99.9 \%$ of the energy is in $[\mu - 3\sigma, \mu + 3\sigma]$. How can I find a similar support for non-central chi-squared distribution with degrees of freedom $\nu$ and non-centrality parameter $\lambda$?
$\begingroup$
$\endgroup$
2
-
2$\begingroup$ By "energy" do you mean "probability"? Assuming that's the case, could you tell us how you computed the answer for the Gaussian distribution--and then explain why the same method would not work for a non-central $\chi^2$ distribution? $\endgroup$– whuber ♦Commented Sep 7, 2015 at 16:17
-
$\begingroup$ For Gaussian distribution, I did not compute the answer but it is known (see tables in en.wikipedia.org/wiki/Normal_distribution for F(μ + nσ) − F(μ − nσ) for different n). $\endgroup$– bassirCommented Sep 7, 2015 at 21:14
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
You do it the same way you could do it for a normal distribution, by using the cdf (cumulative distribution function) $F$ or better, its inverse, the quantile function $Q=F^{-1}$. There is only one complication, the normal distribution is unimodal and symmetric, so it is quite clear that the shortest interval containing a given probability $p$ is the central interval (around the symmetry point) given by $(Q(p/2), Q(1-p/2))$ (assuming $p>1/2$.)
The non-central chi-squared is not symmetric, so the equal-tail interval given above will not be the shortest one (the shortest one will be the HPD interval Proof that the HPD region is the smallest. As an example for the noncentral chisquare, in R, the equal-tail interval is
qchisq(c(0.005, 0.995), df=10, ncp=1)
[1] 2.380353 27.588429
while the HPD interval will be more work, but there is a dedicated R package:
HDInterval::hdi(function(p)qchisq(p, df=10, ncp=1), 0.99)
lower upper
1.660555 25.793967
attr(,"credMass")
[1] 0.99
which must be using numerical optimization.
answered Sep 16, 2019 at 9:14
$\begingroup$
$\endgroup$
HDRs for standard probability distributions can also be computed using the stat.extendpackage. This gives a user friendly output for the region in set form. Here is an example of computing the region for a particular choice of parameters.
#Compute and print HDR for the chi-squared distribution
DF <- 10
NCP <- 1
HDR <- stat.extend::HDR.chisq(0.99, df = DF, ncp = NCP)
HDR
Highest Density Region (HDR)
99.00% HDR for chi-squared distribution with 10 degrees-of-freedom and
non-centrality parameter = 1
Computed using nlm optimisation with 8 iterations (code = 1)
[1.66055476349383, 25.7939667370053]
It can easily be established that the endpoints of the interval have the same density value (to within a small tolerance from the optimisation algorithm) and that the interval has the required coverage level.
#Check endpoints of HDR
dchisq(min(HDR), df = DF, ncp = NCP)
[1] 0.002842709
dchisq(max(HDR), df = DF, ncp = NCP)
[1] 0.002842709
#Check coverage probability
pchisq(max(HDR), df = DF, ncp = NCP) - pchisq(min(HDR), df = DF, ncp = NCP)
[1] 0.99
