최근에 새삼스럽게 '몬테카를로 시뮬레이션'을 다시 파고 있다. 추론과 예측에 효과적인 방법임을 절감했기 때문이다. 복습 삼아 몬테카를로 시뮬레이션의 유명한 예제 중에 '원주율 구하기'를 찾아 봤다가, 재미 삼아 엑셀에서 VBA 없이 구현해 보았다.
난수로 점을 찍어 정사각형 넓이와 해당 정사각형의 한 변의 길이와 동일한 길이의 지름인 원의 넓이의 비로서 원주율 구하기
출처: 몬테카를로 방법으로 원주율 구하기 (http://www.playsw.or.kr/repo/cast/109)
먼저 정말 친절하게 설명한 몬테카를로 방법으로 원주율 구하기(http://www.playsw.or.kr/repo/cast/109) 문서를 꼭 참조하길 바란다. 진국은 이미 저 문서에 있다. 내가 하려는 건, 프로그래밍과 무관하게 살아 온 사람이 엑셀로만 몬테카를로 시뮬레이션을 경험해 보기이다.
1) 우선 정사각형 안에 충분한 개수의 점을 찍어야 한다.
아래 수식을 셀에 복사하고 백만 줄 전체에 붙여 넣어도 된다. ([Ctrl]+[↓]를 알아서 잘 쓰시길.)
X좌표 수식 |
Y좌표 수식 |
=RAND()*(IF(MOD(RANDBETWEEN(1, 2), 2)=0,1,-1)) |
좌동 |
위 수식은 (1,1), (1,-1), (-1, 1), (-1,-1)을 꼭지점으로 하는 정사각형 내에 점을 찍는 취지로서, 이런 식으로 보이게 된다.
X좌표 예시 |
Y좌표 예시 |
0.066007122 |
-0.268583866 |
2) 점과 원 중심 간의 거리를 구한다.
수식 |
예시 |
=SQRT((A1-0)^2+(B1-0)^2) | 0.276575909 |
3) 점이 원 밖인지 원 안인지 구분한다. (지름이 2인 원)
수식 |
예시 |
=IF(C1>1, "밖", "원") | 원 |
아래와 같은 셀을 복사하여 대량으로 붙여 넣으면 된다. PC 성능에 따라 백만 줄 정도 만들 때에는 수십 초가 걸리기도 한다.
◀ | A | B | C | D |
1 | 0.152992 | 0.618691 | 0.637327 | 원 |
2 | -0.33291 | -0.68673 | 0.763168 | 원 |
3 | -0.54012 | -0.8925 | 1.043209 | 밖 |
4) 충분히 복사한 원 안팎에 있는 난수 점들의 비율을 계산하여 대략적인 원주율을 구한다.
시험 삼아 천 개, 십만 개, 백만 개의 점들을 두고 원주율을 계산해 보았다.
회수 |
수식 |
결과 |
오차 |
1,000회 |
=COUNTIF($D1:$D1000, "원")/COUNTA($D1:$D1000)*4 | 3.16 | -0.018407346 |
십만회 |
=COUNTIF($D1:$D100000, "원")/COUNTA($D1:$D100000)*4 | 3.13516 | 0.006432654 |
백만회 |
=COUNTIF($D:$D, "원")/COUNTA($D:$D)*4 | 3.142059 | -0.000466673 |
(=PI()-결과값) |
위 결과는 아름답게도 백만 회의 난수 생성을 거친 결과가 실제 π와 제일 유사하지만, 몇 차례 돌려 본 결과로는 십만 개의 점을 가지고 원주율을 구했을 때의 결과가 π에 더 가까웠던 적도 많았다. 그렇다면 십만 정도로 해도 현실적으로는 효과가 있다고 판단할 만하며, 미심쩍으면 이삽십만 개만 만들어도 충분하다고 본다.
예제 파일은 아래 링크에 두었다.
링크 https://1drv.ms/f/s!AAk_f8uIr69wgpVy 폴더에서 GettingPiWithExcel.xlsx 파일을 열거나 다운로드 하세요.
이렇게 엑셀의 난수 만들기(RAND, RANDBETWEEN) 기능이 꽤 강력하니, 중간과정의 수식을 좀 더하면 프로그래밍 없이도 몬테카를로 시뮬레이션이 가능하다. 막연한 느낌에 휘둘리는 게 싫었다면, 보다 일관적인 예측과 추정으로 이끌어 주는 몬테카를로 시뮬레이션의 세계를 접해 보길 바란다. 재미 들리면 VBA부터 시작하여 여러 프로그래밍 언어로 개발을 시작하게 될지도 모른다. ☺
'BI' 카테고리의 다른 글
4할이 꿈의 타율인 이유는? (0) | 2017.09.18 |
---|---|
분석가에게 멍석 깔아주기 (0) | 2017.01.19 |
SAS에서 이전 연월 문자열 만들기 (0) | 2016.12.01 |