R23. Poisson Ditribution in R

Given a rate of lambda for an event in some interval, the Poisson distribution finds k such events, in the same interval. We know that the expected value of k has to be lambda.

For a particular k, we have the formula lambda^k/factorial(k) multiplied by a normalization constant. We can rewrite lambda^k/factorial(k) as:

(lambda/k)*(lambda/(k-1))*(lambda/(k-2))*...*(lambda/2)*(lambda) for k not 0, and 1 for k equal to 0

The function get_num finds this value multiplied by the normalization constant, which is exp(-lambda). With the normalization, the total probability for all k is 1.

Besides calculating, we can use the dpois function.

# ex23.R

N <- 20

get_num <- function(lambda,kv) {
  coeff <- exp(-lambda)
  if (kv == 0) return(coeff) 
  cv <- lambda/1:kv
  coeff*prod(cv)
}

poisson <- function(lambda, calc) {
  k <- 0:N
  dis.calc <- mapply(get_num,lambda,k)
  dis.R <- dpois(k,lambda)
  if (calc) {
    r <- dis.calc
  } else {
    r <- dis.R
  }
  r
}

c0 <- 0.8
pr <- 6
plot(0:N,poisson(1, FALSE), col = rgb(0,1-c0,c0),
     type = 'o', pch = pr, lty = pr,
     xlab = 'k', ylab = 'p',
     main = "Poisson Distribution (lambda = 1,3,5,7,9)")
for (i in seq(3,9,2)) {
  c0 <- c0-.1
  pr <- (i-1)/2 + 1
  points(0:N,poisson(i, TRUE), 
         col = rgb(0,1-c0,c0), type = 'o', pch = pr, lty = pr)
}
err <- poisson(9,TRUE)-poisson(9,FALSE)
cat('err for lambda of 9 =', sd(err))
# err for lambda of 9 = 1.755738e-17

Output: