phi <- 10 n <- 100000 p <- 8 J <- rpois(n, phi) U <- rchisq(n, p+2*J) mean(U); var(U) p + 2 * phi 2 * p + 8 * phi X2.nc <- rchisq(n, p, 2*phi) mean(X2.nc); var(X2.nc) # cbind(J,U) win.graph(height=5.5, width=7.0) par(mfrow=c(2,1)) ### Histogram of X^2_p(2*phi) (Only those < 60) hist(U[U<=60],breaks=seq(0,60,1),main="Direct Generation of NC Chi-Square", freq=F) lines(seq(0,60,1),dchisq(seq(0,60,1),p,2*phi)) hist(X2.nc[X2.nc<=60],breaks=seq(0,60,1),main="Sampling of NC Chi-Square", freq=F) lines(seq(0,60,1),dchisq(seq(0,60,1),p,2*phi))