Generate positive real numbers from rpois()

Question

I am trying to create a Poisson simulation using rpois(). I have a distribution of two decimal place interest rates and I want to find out if these have Poisson rather than a normal distribution.

The rpois() function returns positive integers. I want it to return two decimal place positive numbers instead. I have tried the following

set.seed(123)
trialA <- rpois(1000, 13.67) # generate 1000 numbers
mean(trialA)
13.22 # Great! Close enough to 13.67
var(trialA)
13.24 # terrific! mean and variance should be the same
head(trialA, 4) 
6 7 8 14 # Oh no!! I want numbers with two decimals places...??????

# Here is my solution...but it has a problem

# 1) Scale the initial distribution by multiplying lambda by 100

trialB <- rpois(1000, 13.67 * 100)

# 2) Then, divide the result by 100 so I get a fractional component
trialB <- trialB / 100
head(trialB, 4) # check results
13.56 13.62 13.26 13.44 # terrific !

# check summary results
mean(trialB)
13.67059  # as expected..great!
var(trialB)
0.153057 # oh no!! I want it to be close to: mean(trialB) = 13.67059

How can I use rpois() to generate positive two decimal place numbers that have a Poisson distribution.

I know that Poisson distributions are used for counts and that counts are positive integers but I also believe that Poisson distributions can be used to model rates. And these rates could be just positive integers divided by a scalar.

Answer 1

If you scale a Poisson distribution to change its mean, the result is no longer Poisson, and the mean and variance are no longer equal -- if you scale the mean by a factor s, then the variance changes by a factor s^2.

You probably want to use the Gamma distribution. The mean of the Gamma is shape*scale and the variance is shape*scale^2, so you have to use scale=1 to get real, positive numbers with equal mean and variance:

set.seed(1001)
r <- rgamma(1e5,shape=13.67,scale=1)
mean(r) ## 13.67375
var(r)  ## 13.6694

You can round to two decimal places without changing the mean and variance very much:

r2 <- round(r,2)
mean(r2) ## 13.67376
var(r2)  ## 13.66938

Compare with a Poisson distribution:

par(las=1,bty="l")
curve(dgamma(x,shape=13.67,scale=1),from=0,to=30,
      ylab="Probability or prob. density")
points(0:30,dpois(0:30,13.67),type="h",lwd=2)