몬테카를로 시뮬레이션 03

Matt Lee·2021년 10월 26일
0
post-thumbnail

이번에는 조건부 확률에 대한 몬테카를로 시뮬레이션을 수행 해 보겠습니다. 시뮬레이션을 위한 상황은 다음과 같습니다.

공평한 A,B 주사위 두 개를 던졌을 때 두 개의 주사위의 합이 7이 나왔다면 A 주사위의 눈이 2일 확률에 대한 몬테카를로 시뮬레이션을 수행 하겠습니다.

NOTE: 위에 경우에 대한 수학적 확률은 Conditional Probability로 아래와 같이 계산이 가능합니다.

P(A=2A+B=7)=P(A=2  A+B=7)P(A+B=7)=P(A=2  2+B=7)P(A+B=7=P(A=2  B=5)P(A+B=7)=P({2,5})P({(1,6),(2,5),(3,4),(4,3),(5,2),(6,1)})=136636=16\begin{aligned} P(A=2|A+B=7) &= \frac{P(A=2 \; \cap A+B=7)}{P(A+B=7)} \\ &= \frac{P(A=2 \; \cap 2+B=7)}{P(A+B=7} \\ &= \frac{P(A=2 \; \cap B=5)}{P(A+B=7)} \\ &= \frac{P(\{2,5\})}{P(\{ (1,6),(2,5),(3,4),(4,3),(5,2),(6,1) \})} \\ &= \frac{\frac{1}{36}}{\frac{6}{36}} \\ &= \frac{1}{6} \end{aligned}

몬테카를로 시뮬레이션을 통해 수확적 확률로 근사 되는지 테스트를 해 보겠습니다.

 

코드 설명
n:= 백만번 시뮬레이션을 수행합니다.
ctr:= 초기값 0을 할당합니다. ctr을 통해 sum이 7이 나온 경우에만 카운트를 갱신 하여 백만번의 카운트를 제어합니다.
trial:= 시행에 대한 시뮬레이션을 정의 합니다. sample 함수를 통해 주사위 두개를 던지는 시행을 정의합니다.
success:= 주사위의 합이 7이 나온 경우 중에 trial의 첫 번째 원소가 2인 경우에만 성공이라고 판단합니다.
simlist[ctr]:= index 기준을 ctr로 하여 success의 값을 원소로 받습니다. while문 실행 후에 simlist는 원소의 갯수가 백만인 벡터가 됩니다.

set.seed(7)
n <- 1000000
ctr <- 0
simlist <- numeric(n)
while (ctr < n) {
    trial <- sample(1:6, 2, replace = TRUE)
    if (sum(trial) == 7) {
        success <- if (trial[1] == 2) 1 else 0
        ctr <- ctr + 1
        simlist[ctr] <- success
    }
}
mean(simlist)
## [1] 0.166645

 
결과해석
mean(simlist)의 결과가 0.166645로 수학적 확률인 16\frac{1}{6}로 잘 근사되었음을 확인 할 수 있습니다. 이론적으로는 시행을 무한번으로 한다면 정확히 16\frac{1}{6}로 수렴 됩니다. 여기서 확인 할 수 있는 재밌는 점은 P(A=2A+B=7)=16=P(A=2)P(A=2|A+B=7) = \frac{1}{6} = P(A=2)로 주사위의 합이 7이라는 정보가 주사위 A가 2인 확률에 아무런 영향을 주지 않았다는 점입니다. 즉, 주사위 A가 2인 사건과 주사위 A와 주사위 B의 합이 7이라는 사건은 서로 독립사건이란 의미입니다.

NOTE:

몬테카를로 시뮬레이션의 3원칙: 시행에 대한 시뮬레이션, 시행의 성공에 대한 판단, 복제 부분을 명확하게 확인하기 위해 코드는 반복문 형식으로 작성 되었습니다. 동일한 시뮬레이션에 대해 vectorized 형식의 코드도 한번 시도해 보세요.

Reference

본 포스팅은 Probability: With Applications and R 2nd Edition by Amy S. Wagaman (Author), Robert P. Dobrow ISBN-13: 978-1119692386 을 참조 하였습니다.

Book link => click here

profile
미국에 서식 중인 응용 수학과 대학원생, 아직은 잉여지만 그래도 행복 :)

0개의 댓글