Random number generator for non-central chi-squared with non-integer dimension

Question

Does someone know of a random number generation algorithm for a non-central chi-squared distribution with a non-integer dimension?

PS By algorithm, I am interested in the detailed procedure or methodology underlying the algorithm.

Answer 1

I am assuming you know how to generate random draws from a central chi-squared distribution, or from its equivalent gamma version; see below for the details. I also suggest possible readings on algorithms for the generation of random variates from a gamma distribution, in case this premise does not hold; see also the references provided in the comments of your post.

Solution 1. As per comments, the function rchisq of R does implement what you are looking for. From the help page of this function:

Usage

rchisq(n, df, ncp = 0)

Arguments

n: number of observations. If length(n) > 1, the length is taken to be the number required.

df: degrees of freedom (non-negative, but can be non-integer).

ncp: non-centrality parameter (non-negative).

Here is a simple R code for generating $10^4$ samples from the chi-square with $\nu = 2.56$ degrees of freedom and non-centrality parameter ncp=6.52.

hist(rchisq(10^4, df=2.56, ncp=6.52))

Solution 2. Expanding over whuber's comment, you can generate numbers from this distribution by noting the fact that its density is an infinite mixture of central chi-squared distributions with Poisson weights. That is, a non-central chi-squared random variable with non-centrality parameter $\lambda$, i.e. $\chi_\nu^2(\lambda)$ has density function

$$ f(x;\nu, \lambda) = \sum_{r=0}^{\infty} \frac{e^{-\lambda/2} (\lambda/2)^r}{r!}f_{\chi_{\nu+2r}^2}(x) = E\left(f_{\chi_{\nu+2R}^2}\right), $$ where $R\sim \text{Poisson}(\lambda/2)$ and $\chi_\nu^2$ denotes the central chi-squared distribution.

To generate a random variate $x$ from the $\chi_\nu^2(\lambda)$ distribution then do:

• draw $r$ from $\text{Poisson}(\lambda/2)$
• draw $x$ from $\chi^2_{\nu +2r}$

Here is a simple R code for all this.

# implement a simple function 
my_rchisq <- function(df, ncp) {
  if(ncp==0)
    stop("Please use rchisq.")
  n = rpois(1, lambda = ncp/2)
  out = rchisq(1, df= df + 2*n, ncp = 0)

  return(out)
}

# call the function N times
gen2 <- sapply(1:N, function(x) my_rchisq(dof,ncp))

# draw a histogram and compare it with true pdf
hist(gen2,breaks = 30, probability = TRUE)
plot(function(x) dchisq(x, df=dof, ncp=ncp), 
n=200, 
add=TRUE, 
lwd=2, xlim = c(0,50))

Comments The code presumes you know how to generate random variates from $\chi_{\nu}^2$, for any real $\nu>0$ (again, see rchisq in R). However, since this distribution is a special case of a gamma distribution, i.e. if $Y \sim \chi_{nu}^2$ then $Y$ has a gamma distribution with shape parameter $\nu/2$ and scale parameter equal to 2, it boils down to generating random draws from the gamma distribution (rgamma in R).

Since you are interested in $0<\text{ncp}<1$, this means that you have to draw from the gamma distribution with a shape parameter between zero and one. It turns out that the generation of random draws, in this case, may be troublesome due to close-to-zero values. The documentation of rgamma reads

Note that for smallish values of shape (and moderate scale) a large parts of the mass of the Gamma distribution is on values of xx so near zero that they will be represented as zero in computer arithmetic. So rgamma may well return values which will be represented as zero. (This will also happen for very large values of scale since the actual generation is done for scale = 1.)

Now, if you also wish to know how to generate random draws from a gamma distribution, I suggest looking at the references on the help page of rgamma, which are specific to this issue. In particular, for the problem of shape between 0 and 1 check

Ahrens, J. H. and Dieter, U. (1974). Computer methods for sampling from gamma, beta, Poisson and binomial distributions. Computing, 12, 223–246.