[R] unexpected results involving the skellam distribution

[Ernest Adrogué]

Hi all,
I am new to R. I made this script that plots the
deviation of a skellam-distributed random variable
with respect to the skellam distribution.

I would expect to get random errors, but the plot
sistematically shows a non-random pattern (first a
peak and then a low). I don't know how to explain
this. Can somebody enlighten me??


require(skellam)

n <- 5000
lam1 <- 5.45
lam2 <- 2.78
p1 <- rpois(n, lam1)
p2 <- rpois(n, lam2)

rho <- cor(p1,p2)

z <- hist(p1-p2, plot=FALSE)

mu1 <- lam1 - (rho*sqrt(lam1*lam2))
mu2 <- lam2 - (rho*sqrt(lam1*lam2))

x <- z$breaks[1:length(z$breaks)-1]
y <- dskellam(x, mu1, mu2)

plot(x,z$density-y)


-- 
Cheers.
Ernest