(** 인터넷 익스플로러에서 보면 파이썬 코드 줄이 제대로 적용되어 보이지 않을 수 있습니다.**)
몬테카를로 방법은 반복해서 무작위 값을 만들어 원하는 값을 계산하는 유형의 알고리즘들을 가리킨다. 스타니스와프 울람이 모나코의 유명한 도박의 모시 몬테카를로의 이름을 따서 지었다. 이번에는 몬테카를로 알고리즘 두 가지를 파이썬으로 작성 및 실행하여 각각의 경우에 원하는 수치를 계산할 수 있음을 확인한다.
먼저 파스칼의 문제를 몬테카를로 방법으로 해결해보자. 파스칼의 문제는 두 개의 주사위를 24번 반복해서 던지는 게임에 관한 것이다. 동시에 두 개의 주사위가 모두 6이 나오면 이 게임에서 이긴다. 파스칼의 문제는 이 게임에서 이길 확률을 구하는 것이다.
확률 과목을 수강한 학생이라면 다음과 같이 진행할 수 있다.
- 두 개의 주사위가 모두 6이 나올 확률 : 1/6 x 1/6 = 1/36
- 반대로 두 개의 6이 나오지 않을 확률 : 1 - 1/36 = 35/36
- 24번 반복해서 두 개의 6이 나오지 않을 확률 : (35/36)x ... x(35/36) = (35/36)^24
- 24번 반복했을 때 두 개의 6이 나올 확률 : 1 - (35/36)^24
파스칼의 문제에서 이길 확률은 1-(35/36)^24이다. 파이썬으로 이 값을 계산해보면 그 확률은 다음과 같다.
>>> 1 - (35/36.0)**24
0.4914038761309034
위와 같이 가능한 경우의 수를 따져 확률을 구하는 방법을 모르는 시절에는 이 것을 어떻게 계산했을까? 처음 파스칼의 문제를 따져보던 시기는 확률에 대한 수학이 없었다. 이렇게 경우의 수를 따져 확률을 계산하는 대신 몬테카를로 방법으로 이 확률(의 근사값)을 계산할 수 있다.
기본적인 아이디어는 주사위 2개를 던져 동시에 2개의 6이 나올 때까지 24번 던지는 게임을 반복해서 실제로 확률을 시뮬레이션해보는 것이다.
import random
def rollDie():
return random.choice([1,2,3,4,5,6])
def checkPascal(numTrials):
numWins = 0.0
for i in range(numTrials):
for j in range(24):
d1 = rollDie()
d2 = rollDie()
if d1 == 6 and d2 == 6:
numWins = numWins + 1
break;
print('24번 중 2개의 6이 나와 이기는 확률', numWins/numTrials)
위 파이썬 프로그램에서 rollDie()는 주사위를 던져 1부터 6 사이의 숫자를 구하고, checkPascal(numTrials) 함수는 numTrials만큼 파스칼 문제의 게임을 반복해서 게임에 승리할 확률을 시뮬레이션으로 구한다. 이 함수들을 자세히 이해하기 앞서 일단 실행하고 파스칼 문제의 확률을 구해보자.
>>> checkPascal(1000000)
24번 중 2개의 6이 나와 이기는 확률 0.492315
>>> checkPascal(10)
24번 중 2개의 6이 나와 이기는 확률 0.4
>>> checkPascal(5)
24번 중 2개의 6이 나와 이기는 확률 0.6
>>>
백 만번 반복해서 checkPascal 함수를 실행하면 0.492315 확률을 얻는데 이 수치는 경우의 수를 따져 계산한 확률 0.4914038761309034과 상당히 비슷한 것을 알 수 있다.
checkPascal() 함수는 승리 횟수를 저장할 변수 numWins를 초기화하며 시작한다. 두 개의 내포된 반복문에서 바깥쪽 for 반복문은 지정한 반복 횟수만큼 게임을 되풀이하고, 안쪽 for 반복문은 24번 두 개의 주사위를 던지는 과정을 시뮬레이션한다. d1과 d2는 두 개의 주사위를 던져 얻은 숫자이고, 둘 다 6이면 승리 횟수 numWins을 1증가시키고 break문으로 안쪽 for 반복문을 마친다. 즉, 24번 중 아직 남은 경우가 있더라도 더 진행하지 않는다.
checkPascal() 함수에서 반복 횟수를 10이나 5로 크게 줄이면 시뮬레이션 결과로 이길 확률을 각각 0.4와 0.6으로 얻게 된다. 이 확률 값은 수식으로 계산한 확률과 다소 차이가 있다. 특히 5회 반복하는 경우 파스칼 문제에서 승리할 확률이 질 확률보다 높다는 편향된 결과를 보여준다.
다음 문제는 몬테카를로 방법으로 원주율 파이(pi)를 계산하는 것이다. 중학교 수학에서 원주율이 3.141592... 이라는 것을 배워서 이미 알고 있긴하다.
원주율을 몬테칼르로 방법으로 계산하는 방법의 기본적인 아이디어는 한 변의 길이가 2인 정사각형과 그 내부를 채우는 원을 그린다. 원의 지름이 2가 되기 때문에 반지름은 1이 된다. 이 정사각형 안에 무작위로 점을 찍어서 원 내부에 포함된 점의 수(M)와 전체 점의 수(N)의 비율이 원의 면적과 정사각형의 면접의 비율과 일치할 것이다.
- M / N = pi x r^2 / (2r x 2r) = pi / 4
- 따라서 pi = 4 x M / N이다.
즉, 몬테카를로 방법으로 M과 N의 비율을 구하면 pi를 계산할 수 있다.
def throwNeedles(numNeedles):
inCircle = 0
for Needles in range(1, numNeedles + 1):
x = random.random()
y = random.random()
if ( (x*x + y*y)**0.5 <= 1.0 ): # x,y가 원 안에 있다면
inCircle = inCircle + 1
return 4 * (inCircle / float(numNeedles))throwNeedles(numNeedles) 함수는 지정된 갯수 numNeedles의 바늘을 정사각형 안에 던져 찍힌 위치가 원 내부에 있는 경우를 따져서 전체 바늘 던진 횟수에 대한 비율(4배)을 계산하여 원주율의 근사값을 구한다. 일단 실행해서 원주율을 구해보자. 천 만번 반복해서 시뮬레이션하면 우리가 아는 원주율과 상당히 가까운 값을 얻을 수 있음을 알 수 있다.
>>> throwNeedles(1000000)
3.139532
>>> throwNeedles(10000000)
3.1414944
>>>
throwNeedles() 함수는 먼저 원 내부에 포함된 좌표의 수를 저장할 inCircle을 0으로 초기화한다. 지정한 숫자 numNeedles 만큼 다음 과정을 반복한다. 던진 바늘이 찍힌 위치 x와 y를 무작위로 구한고, x, y 좌표가 원 내부에 있으면 (또는 이 좌표에서 원점까지의 거리가 반지름의 길이 이하이면) 원 내부에 포함된 좌표의 수를 1 증가시킨다. 이 반복된 과정이 끝나면 전체 좌표들의 수 numNeedles 대비 원 내부에 포함된 좌표들의 수 inCircle을 구하고 앞서 설명한 아이디어대로 이 비율에 4를 곱하여 원주율을 구해 리턴한다.
사실 우리가 원주율을 알고 있기 때문에 몬테카를로 시뮬레이션 결과가 원하는 수치를 낸다는 것을 이해한다. 하지만 만일 우리라 원주율을 모른다면 앞에서 얻은 3.139532와 3.1414944 중 어느 것을 더 선호해야하는지 판단하기 어렵다. 무작정 많이 반복하면 더 좋은 근사값을 얻을 수도 있지만 그 만큼 계산 시간이 오래걸린다.
즉, 몬테카를로 시뮬레이션에서 몇 번 반복하면 충분히 납득할만한 수치를 얻을 수 있는지 궁금하다. 한 가지 요령은 원하는 정밀도(precision)을 지정하고 시뮬레이션에서 얻은 수치들의 표준 편차가 그 정밀도 이하가 될 때까지 반복하는 것이다.
다음 파이썬 프로그램은 이 전략에 따라서 동적으로 반복 횟수를 결정하여 원주율 근사값을 구한다. stdDev(X)는 리스트 X의 원소들에 대한 표준 편차를 구하는 함수이다. estPi(precision, numTrials)는주어진 정밀도 precision에 도달할 때 까지(예: 표준편차가 정밀도/2보다 작을 때까지) numTrials만큼 반복해서 원주율 근사값들을 각각 구한다. 만일 앞서 시뮬레이션에서 얻은 원주율들의 표준편차가 precision에 도달하지 못하면 1000으로 시작한 바늘의 갯수를 2배로 늘려 다시 numTrials 만큼의 원주율 근사값들을 구하도록 시뮬레이션을 다시 진행한다.
def stdDev(X):
mean = float(sum(X))/len(X)
tot = 0.0
for x in X:
tot += (x - mean)**2
return (tot/len(X))**0.5
def getEst(numNeedles, numTrials):
estimates = []
for t in range(numTrials):
piGuess = throwNeedles(numNeedles)
estimates.append(piGuess)
sDev = stdDev(estimates) # 표준편차(standard deviation, stdDev)
curEst = sum(estimates)/len(estimates)
print('Est. = ' + str(round(curEst, 5)) +
', Std. dev. = ' + str(round(sDev, 5)) +
', Needles = ' + str(numNeedles))
return (curEst, sDev)
def estPi(precision, numTrials):
numNeedles = 1000
sDev = precision # 표준편차(standard deviation, sDev)
while sDev >= precision/2.0:
curEst, sDev = getEst(numNeedles, numTrials)
numNeedles = numNeedles * 2
return curEst
estPi() 함수를 실행해서 원하는 정밀도의 원주율을 구해보자. 처음에는 정밀도 0.01, 100개의 원주율 근사값들을 구해서 최종 원주율 근사값을 결정한다. 3.141650... 이라는 좋은 결과를 얻었고, 맨 마지막에 표준 편차 값이 0.00445로 0.01/2 보다 작아서 시뮬레이션을 정지하였음을 알 수 있다. 이 때 던진 바늘의 수는 128000개다.
두번째에는 0.001로 정밀도를 더 높이고, 구할 원주율 근사값들의 수는 100개로 동일하다. 최종 결과는 3.141601...로 더 좋은 결과이다. 이때 사용한 바늘의 수는 16,384,000개다. 지정한 정밀도 0.001/2보다 작은 표준 편차 0.00037을 갖는 100개의 원주율 근사값들을 시뮬레이션에서 얻었기 때문에 종료하였다.
>>> estPi(0.01, 100)
Est. = 3.14408, Std. dev. = 0.0572, Needles = 1000
Est. = 3.13976, Std. dev. = 0.03623, Needles = 2000
Est. = 3.13968, Std. dev. = 0.02542, Needles = 4000
Est. = 3.13962, Std. dev. = 0.01816, Needles = 8000
Est. = 3.14117, Std. dev. = 0.01302, Needles = 16000
Est. = 3.14202, Std. dev. = 0.00951, Needles = 32000
Est. = 3.14179, Std. dev. = 0.00688, Needles = 64000
Est. = 3.14165, Std. dev. = 0.00445, Needles = 128000
3.1416506249999996
>>> estPi(0.001, 100)
Est. = 3.15016, Std. dev. = 0.05635, Needles = 1000
Est. = 3.1435, Std. dev. = 0.03882, Needles = 2000
Est. = 3.13984, Std. dev. = 0.02466, Needles = 4000
Est. = 3.14126, Std. dev. = 0.01774, Needles = 8000
Est. = 3.14261, Std. dev. = 0.01428, Needles = 16000
Est. = 3.14215, Std. dev. = 0.00829, Needles = 32000
Est. = 3.14227, Std. dev. = 0.00657, Needles = 64000
Est. = 3.14189, Std. dev. = 0.00426, Needles = 128000
Est. = 3.14139, Std. dev. = 0.00336, Needles = 256000
Est. = 3.14148, Std. dev. = 0.00253, Needles = 512000
Est. = 3.14164, Std. dev. = 0.0018, Needles = 1024000
Est. = 3.14163, Std. dev. = 0.00128, Needles = 2048000
Est. = 3.14165, Std. dev. = 0.00079, Needles = 4096000
Est. = 3.14168, Std. dev. = 0.00057, Needles = 8192000
Est. = 3.1416, Std. dev. = 0.00037, Needles = 16384000
3.1416019580078114