본문 바로가기

코드이야기/R

[R프로그래밍] 몬테카를로 시뮬레이션으로 원주율(π) 구하기


몬테카를로 시뮬레이션이란 난수를 이용한 시뮬레이션을 여러 번 반복하여 문제의 근사해를 구하는 방법을 말합니다. 수소폭탄의 개발자인 스타니스와프 울람이라고 하는 폴란드계 미국인 과학자에 의해 모나코의 유명한 도박도시 몬테카를로의 이름을 본따 만들어진 것입니다. 수학이나 물리학에 주로 사용되며 엔리코 페르미가 중성자의 특성을 연구하기 위해 이 방법을 사용하기도 하였습니다.



지난 포스팅에서 시뮬레이션의 기본인 난수생성방법을 알아보았습니다. 여기서는 시뮬레이션이 정의와 프로세스, 그리고 R을 이용해 몬테카를로 시뮬레이션을 R로 구현하는 방법에 대해 이야기 하겠습니다.

시뮬레이션의 정의와 프로세스

위키피디아에는 '시뮬레이션은 실제로 실행하기 어려운 실험을 간단히 행하는 모의실험을 뜻한다'라고 정의되어 있습니다. 여기서 핵심어는 '모의'입니다. 한자어로는 '본뜰 모(模)'에 '비길 의(擬)'를 사용합니다. 실제의 것을 흉내내어 그대로 해본다는 말입니다.


실제로 해보면 되는 것이지 왜 굳이 모의실험을 해볼까요? 세상은 직접 해보기에는 너무 복잡합니다. 그 실험이 대상이 사람의 생명을 다루는 일이라면 모의실험 자체가 불가능합니다. 모형을 적절히 설정하고, 컴퓨터의 힘을 빌리면 이런 실험도 만족할 만한 결과를 얻어낼 수 있습니다. 시뮬레이션의 목적을 설정하고 설계하는 과정이 실제 시뮬레이션보다 중요합니다. 시뮬레이션은 컴퓨터의 힘을 빌릴 수 있지만 시뮬레이션 설계는 우리의 머리로 해야하기 때문입니다.


이런 시뮬레이션을 하기위해서는 일반적으로 아래와 같은 과정을 따릅니다.

1. 목적정의 : 왜 시뮬레이션을 하려고 하는지 그것을 통해 해결하려는 문제를 정의합니다.

2. 상황정의 : 목적에 따라 시뮬레이션 해야할 상황을 정의합니다.

3. 흐름도 작성 : 시뮬레이션이 흐름을 개략적으로 그려봅니다.

4. 알고리즘 결정 : 흐름도를 컴퓨터로 처리할 수 있도록 모델을 만듭니다.

5. 시뮬레이션 수행 : 결과를 얻기 위해 알고리즘을 코딩하고 수행합니다.

6. 결과검토 : 필요한 통계처리를 수행하고 결과를 분석해 결론을 도출합니다.

7. 문서화 : 도출한 결론을 그래프와 같은 시각화 도구들을 포함해 문서화 합니다.

R을 이용해 몬테카를로 시뮬레이션 하기

이제 몬테카를로 시뮬레이션을 실제로 수행해 보도록 하겠습니다. 오늘 해볼 작업은 몬테카를로 시뮬레이션을 통해 원주율(π) 구하기 입니다. 우리는 π가 3.141592... 로 나가는 무리수라는 사실을 알고 있습니다. 그런데 이 값은 어떻게 구하는 지는 크게 생각해 보지 않았습니다. 다만, 반지름이 r인 원의 둘레는 2πr이고 면적은 πr2이라는 사실만 기억하고 있습니다.


π는 원주와 지름의 비로 정의합니다. 이 π는 지름의 크기와 관계없이 항상 일정한 특성을 가지고 있고, 그 비례 상수를 그리스어 περίμετρο(perimeter)이 머릿글자를 따 그리스 문자 π로 나타내는 것입니다.


π를 R로 구하기 위해 시뮬레이션 과정을 정의해 보겠습니다.

1. (1, 1), (1, -1), (-1, 1), (-1, -1)을 정점으로 하는 평면을 생각합니다.

2. 이 평면에 (x좌표, y좌표)를 나타내는 두개의 난수를 임의로 생성합니다.

3. 이 점이 원점으로부터 거리가 1이내인지 판단하고 맞으면 카운트하고, 아니면 버립니다.

4. 수많은 점에 대해서 3번의 작업 수행 후 전체 점에서 1이내에 있는 점의 비율을 구합니다.

5. 이 비율이 면적이 4인 정사각형에 대한 원의 면적과 같으므로 이비율을 4배 합니다.

6. 5에서 구한 숫자가 원주율에 근사하게 됩니다.

이제 코드로 구현해 보겠습니다.

> pi <- function(n) {
+     counter <- 0
+     for (i in 1:n) {
+         coordinates <- runif(2, min=-1, max=1)
+         if (sqrt(coordinates[1]^2 + coordinates[2]^2) < 1) {
+             counter <- counter + 1
+         }
+     }
+     return(4 * counter / n)
+ }
> 
> pi(1)
[1] 4
> pi(10)
[1] 2.4
> pi(100)
[1] 3.24
> pi(1000)
[1] 3.128
> pi(10000)
[1] 3.12
> pi(100000)
[1] 3.13352
> pi(1000000)
[1] 3.14234
> pi(10000000)
[1] 3.141888

몬테카를로 시뮬레이션이 금방 이해하기 쉽지는 않습니다. 또한, 시뮬레이션을 위한 모형에 대한 알고리즘을 생각하고 코드로 구현하는 것 역시 간단한 작업은 아닙니다. 하지만, 공부해 두면 충분히 효과를 주는 도구임은 분명합니다. 사진: NicoCanali