Comparing simulated to exact distributions

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.

• Could you post your code and output for readers to check the issue? Thank you. Sep 27 '15 at 11:34

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)


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.

• Thanks for your answer. It appears that R's chi-squared RNG is indeed highly error prone for extreme parameters. This link has a table comparing the errors for R & C++ Boost RNGs: boost.org/doc/libs/1_39_0/libs/math/doc/sf_and_dist/html/… Sep 27 '15 at 11:50
• Interesting like but what do you call extreme? I went up to $\lambda=10^6$ and could not spot a difference between histogram and true pdf. Sep 27 '15 at 12:02
• Since this is a simulated stochastic process, the error accumulates over several time steps. As a result, there's a substantial difference between the histogram & true pdf at termination. Sep 27 '15 at 12:16
• Have you tried the alternative simulation of a central $\chi^2_{p-1}$ plus an independent squared normal? Sep 27 '15 at 14:37
• That is not possible for my particular study, as there would be 2 sources of randomness. Sep 27 '15 at 15:26