에 관한 몬테카를로 시뮬레이션을 구현 하려면 다음의 3가지 원칙을 따라야 합니다.
시행에 대한 시뮬레이션: 컴퓨터의 random numbers를 사용해서 model이나 translate를 구현합니다. 컴퓨터의 random numbers를 사용해서 한번 실험을 수행 한 것을 시행이라고 정의합니다.
시행의 성공에 대한 판단: 시행의 결과를 가지고 사건 가 발생 했는지 안했는지의 여부를 판단합니다. 사건 가 발생했다고 판단 되는 경우에 이것을 성공이라고 정의 합니다.
복제: 위의 1,2번 두 단계를 가능한한 많이 반복합니다. 성공적 시행의 비율이 의 시뮬레이션 된 추정값입니다.
NOTE: 앞뒤가 나올 가능성이 동일한 동전에 대해서 3번 동전을 던져서 전부 앞면이 나올 확률 수학적 확률은
코드 설명
시행에 대한 정의 부분입니다. trial := random sample {0,1}, 3번 연속 복원 추출
trial <- sample(0:1, 3, replace = TRUE)
코드 설명
성공에 대한 정의 부분입니다. trial의 합이 3이면 1 아니면 0
if (sum(trial) == 3) 1 else 0
코드 설명
복제에 대한 정의 부분입니다. for문을 통해서 동전을 3번 던지는 실험을 백만번 수행 합니다. 수행 중에 합이 3인 경우만 simlist 벡터에 값을1로 할당합니다. 그 후 mean 함수를 통해서 동전이 3번 나온 경우(즉, 성공한 경우)에 대한 proportion을 계산 합니다. 계산 결과는 수행을 많이 하면 할 수록(n 값을 크게 잡을 수록) 수학적 확률인 에 수렴합니다. 수렴하는 수리적 이유는 the law of large number(큰수의 법칙)로 설명 할 수 있습니다. 이 내용은 추후에 다루 겠습니다.
n := 시행 횟수 정의
simlist := 성공 및 실패에 대한 정보를 원소로 저장할 벡터 정의
mean(simlist) := 앞면이 3번 나옹 경우의 proportion 계산
set.seed(7)
n <- 1000000
simlist <- numeric(n)
for (i in 1:n) {
trial <- sample(0:1, 3, replace = TRUE)
success <- if (sum(trial) == 3) 1 else 0
simlist[i] <- success
}
mean(simlist)
결과해석
mean(simlist)의 결과가 0.125775로 수학적 확률인 0.125로 소수점 이하 세 째 자리 까지 잘 근사되었음을 확인 할 수 있습니다. 이론적으로는 시행을 무한번으로 한다면 정확히 0.125로 수렴 됩니다.
위의 코드는 각 단계별로 시행 -> 성공 -> 복제에 대한 부분을 파악 하기 쉽게 for문을 사용 했는데요.
R의 vectorized 연산(선형대수학의 element-wise operation과 유사: 예를 들어 dot product, scalar vector multiplication etc.)을 이용 하기 위해 다음과 같이 수정할 수 있습니다.
simdivis <- function() {
trial <- sample(0:1, 3, replace = TRUE)
if (sum(trial) == 3) 1 else 0
}
mean(replicate(1000000, simdivis()))
본 포스팅은 Probability: With Applications and R 2nd Edition by Amy S. Wagaman (Author), Robert P. Dobrow ISBN-13: 978-1119692386 을 참조 하였습니다.
Book link => click here