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$?

  • 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
    Sep 7 '15 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$
    – bassir
    Sep 7 '15 at 21:14

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 
[1] 0.99

which must be using numerical optimization.


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)

        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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.