[확률과 통계적 추론] 2-3.1 포아송분포 접근 with R

 

저번 포스팅(https://moogie.tistory.com/67)에서 포아송 분포를 유도하였습니다.

 

포아송 분포를 유도할때 전제 조건으로는 전체구간을 1로 설정하여 구간을 n등분하였고, 각 하위구간이 서로 독립이며 확률이 $\frac{\lambda}{n}$인 베르누이 시행을 따른다고 가정했습니다.

이럴 때 n이 무한대로 발산한다면 각 하위구간에서 성공 횟수를 합한 값들의 분포가 포아송분포를 따른다고 볼 수 있습니다.

이런 방식이 통해 실제 포아송 분포와 유사한지 확인하기 위해 시뮬레이션을 돌려보았습니다.

# Approximately Poisson Distribution
# 개념 
# 1 unit의 구간을 n등분 후 각 독립적인 subinterval이 확률이 lambda/n인 bernoulli trial이라고 생각
# n이 무한에 가깝다면 possion(lambda)와 유사한 분포를 가짐

bernoulli_to_possion <- function(times, n, lambda){
  res <- vector(mode="integer", length=as.integer(times))
  for(i in 1:times){
    sum = 0
    for(j in 1:n){
     sum = sum + sample(c(0,1), size=1, prob=c(1-lambda/n, lambda/n))
    }
    res[i] <- sum
  }
  return(res)
}

 

n=1000인 시뮬레이션으로 얻은 결과와 $\lambda=5$인 포아송분포에서 얻은 결과를 아래와 같이 막대그림으로 표현하였습니다.

library(tidyverse)
a <- bernoulli_to_possion(times=1000, n=1000, lambda=5)
tibble(num=a) %>% ggplot(mapping=aes(x=num)) + geom_bar() +
  xlab(NULL) + ylab("Count")
  
real <- rpois(n=1000, lambda=5)
tibble(num=real) %>% ggplot(mapping=aes(x=num)) + geom_bar() +
  xlab(NULL) + ylab("Count") + ggtitle("Real Poisson")

• 실제 포아송 분포와 유사한 결과를 보이는 것을 알 수 있습니다.

 

또한,  포아송분포는 평균과 분산이 같다는 특징이 있는데 n과 times, 그리고 포아송 모수를 다르게 해서 평균과 분산이 비슷한지 확인해 봤습니다.

res <- crossing(times=c(25,250), n=c(25, 100, 1000), lambda=c(3,5,10)) %>% 
  mutate(res = pmap(.l=list(times=times, n=n, lambda=lambda), .f=bernoulli_to_possion))

res %>% mutate(mean = map_dbl(.x=res, .f=mean),
               variance = map_dbl(.x=res, .f=var))

• crossing 함수를 사용해서 times, n, lambda에 해당하는 각 값들의 Cartesian Product 하였으며  res변수에 각 인스턴스에 해당되는 포아송분포를 따르는 난수를 list형태로 저장하였습니다.

난수의 수(times)와 구간수(n)에 상관없이 평균은 대체로 잘 추정하나 n이 작을수록 분산을 과소추정하는 경향이 있네요. 다만 실제 전제조건인 n이 무한대로 갈수록 포아송분포에 근접한다는 조건과 같이 n이 커질수록 평균과 분산이 유사해지는 것을 확인할 수 있습니다.