Comparing simulated to exact distributions

Question

I am looking to simulate a stochastic process using Monte Carlo methods. The stochastic RV is known to have a non-central $\chi^2$ distribution, so I draw pseudo-random numbers using R's rchisq function.

When I compare the probability density of the exact distribution to empirically measured probabilities from the simulated distribution, the two distributions correspond to each other very well for "benign" parameters. However, for extreme parameters, I find there is a major discrepancy between the two distributions, especially near 0 and on the right tail. This issue remains even when pushing the Monte Carlo parameters (time steps, number of scenarios) to very high levels.

It is unclear to me whether the issue lies with the random number generator or with the iterative procedure used for the simulation. Has someone come across such issues before, for any stochastic process whose distribution is known? Is there a way to use importance/rejection sampling or any other techniques to help get a better match for all parameters? Any references which deal with this issue would be much appreciated.

Answer 1

After running the simulation myself, I find no indication that R rchisq function is not an accurate random generator. Here is an illustration for 3 degrees of freedom and a range of non-centrality parameters:

Produced by the following R code:

 ncp=c(1,5,10,50,100,500)
    for (i in 1:6){
    x=rchisq(1e6,df=3,ncp=ncp[i])
   hist(x,prob=TRUE,nclass=1234)
    curve(dchisq(x,df=3,ncp=ncp[i]),add=TRUE)}

The warning in the R documentation:

The code for non-zero ‘ncp’ is principally intended to be used for moderate values of ‘ncp’: it will not be highly accurate, especially in the tails, for large values.

must correspond to much higher values of 'ncp' as I did not spot any discrepancy even for ncp=1e12.

As a side remark, an reasonable alternative to using rchisq with non-zero ncp is to represent the $\chi^2_d(\lambda)$ variable as the sum of a $\chi^2_d$ variable and of a $\mathcal{N}(\sqrt{\lambda},1)$ variable.