R 기하분포

기하분포

본 내용은 기하분포 내용을 R로 시뮬레이션하고, 그 결과를 이론값과 비교해본다.

예를 들어, 승부차기를 한다. 이때 골을 넣을 확률이 0.2라고 하자. 그럼, 실패할 확률은 0.8이 된다.

그럼, 골 넣을 때까지 차다고 한다면 평균 몇 번을 차면 되는가? 그냥 상식적으로 보면 5번 정도 차면 1번 들어가지 않을까?

자, 이걸 일반화하면 이렇게 된다.

• 골 넣은 확률 = 성공확률 = p = 0.2
• 노골 = 실패확률 = 1 - p = q = 0.8
• 시도 횟수 = r
• 1번에 성공할 확률 = P(X = 1) = 0.2
• 2번에 성공할 확률 = P(X = 2) = 실패확률 * 성공확률 = 0.8 * 0.2
• 3번에 성공할 확률 = P(X = 3) = 실패확률 * 실패확률 = 0.8 * 0.8 * 0.2
• r번에 성공할 확률 = P(X = r) = 실패확률^(r-1) * 성공확률 = 0.8^(r-1) * 0.2 = q^(r-1) * p
• 성공할 때까지 차는 평균 횟수 = 기대치 = 합(확률변수값 * 확률) = E(x) = sum (r * P(X = r))

골 넣을 확률이 0.2 이니까 아래와 같은 숫자 집합에서 1을 뽑을 확률과 같다.

goal <- c(1, 0, 0, 0, 0)

• 여기서 1번만 뽑을 경우를 1000번 수행해보고 1이 나온 경우가 몇번인지 세보자.
• 그리고 2번만 뽑을 경우를 1000번 수행해보고 마지막에 1이 한번만 나온 경우가 몇번인지 세보자.
• 마찬가지로 30번을 뽑을 경우를 1000번 수행해보고 마지막에 1이 한번만 나온 경우가 몇번인지 세보자.

library(ggplot2)

set.seed(1234)

goal <- c(1, 0, 0, 0, 0)
getGoal <- function (n, r)
{
  # r 횟수째에 성공한 건수
  count <- 0
  
    for(i in 1:n)
    {
      # r 횟수 만큼 복원 추출한다.
      a <- sample(goal, r, replace=T)
      
      # r번째에 성공했다면
      # 만약 합이 1이고 맨 끝에 숫자가 1이라면
      if(sum(a) == 1 && a[r] == 1)
      {
        count <- count + 1
      }
    }
    
  return (count)
}

# 100번 시도해서 1번째에 성공한 경우
a <- getGoal(1000, 1)
a

## [1] 190

success_count <- c()
# 1번째에 성공하는 경우, 2번째에 성공하는 경우
# , 3번째에 성공하는 경우... 30번째에 성공하는 경우
for(i in 1:30)
{
  # 100번 시도 했을 때
  a <- getGoal(1000, i)
  success_count <- c(success_count, a)
}

success_count

##  [1] 205 163 117 103  94  85  50  35  40  32  25  13  12  12   6   4   4
## [18]   5   1   4   2   4   1   0   2   1   0   0   0   0

위의 내용을 확률 분포로 살펴보자. 확률분포로 그릴때는 상대도수로 그리는게 좋다.

total_success_count = sum(success_count)
success_count_rel = success_count / total_success_count
plot(success_count_rel)

1000번을 차면 첫번째에 골을 넣을 가능성(확률)은 약 200번이 된다.

이와 같은 분포를 기하분포라고 한다.

기하분포의 이론적 확률은 다음과 같이 정의된다.

확률 P(X = r) = q^(r - 1) * p

이렇게 정의된 내용을 위의 모의실험된 확률분포와 겹쳐서 그려본다.

# 이론적 확률변수의 값을 구한다.
p <- 0.2    # 성공확률
q <- 1 - p  # 실패확률

r <- seq(1:30) # 30번 시도
P_r <- q^(r-1) * p # r번째에 성공할 확률값

# 실험값과 이론값을 데이터프레임으로 구성한다.
res_simulation <- data.frame(x = r, y = success_count_rel, type='simulation')
res_theory <- data.frame(x = r, y = P_r, type="theory")
res <- rbind(res_simulation, res_theory)

# 그래프로 확인한다.
g <- ggplot(res, aes(x = x, y = y, color=type)) + geom_line()
g

기하확률분포는 r = 1 일 경우 가장 높은 확률값을 가진다. 시도 횟수가 많아질 수록 확률값은 기하급수적으로 작게 된다.

이론적인 기하분포의 기대치는 E(X) = 1 / p가 되고 분산은 Var(X) = q / p^2이 된다.