IDL로 Monte Carlo 방법을 테스트해봅시다

몬테카를로 방법(Monte Carlo method)은 난수를 이용하여 확률적으로 값을 계산해내는 알고리즘을 뜻하는 용어로서 수학, 물리학, 공학 등 많은 분야에서 널리 사용되어오고 있는 기법입니다. 이 기법에 대한 가장 일반적인 예제라면 아무래도 원(Circle)의 면적 또는 원주율을 계산하는 예제일 것 같습니다. IDL에서도 이러한 예제를 직접 구현하여 테스트해보는 것이 가능합니다. 그래서 나름 재미있는 소재가 될 것 같아서 그 과정을 한번 다뤄보기로 하겠습니다. 먼저 원의 1/4만큼에 해당되는 원호 궤적을 그리기 위하여 다음과 같이 코드를 작성해보았습니다.

angles = [0:90:0.1]

HELP, angles

xs = COS(angles*!dtor)

ys = SIN(angles*!dtor)

win = WINDOW(DIMENSIONS=[600, 600], /NO_TOOLBAR)

pl = PLOT(xs, ys, THICK=2, ASPECT_RATIO=1, MARGIN=0.08, $

  FONT_SIZE=11, /CURRENT)

여기서 angles는 0~90까지 0.1 간격의 각도(Degree) 값들에 해당됩니다. 이렇게 생성된 총 901개의 각도값들에 대하여 COS을 취하여 산출한 X좌표 값들을 xs 배열로 그리고 SIN을 취하여 산출한 Y좌표 값들을 ys 배열로 정의하였습니다. 이 xs, ys를 PLOT 함수에 투입하여 원의 1/4 궤적을 표출한 모습은 다음 그림과 같습니다.

이 그림에서 볼 수 있듯이 xs, ys 각각은 0~1의 범위의 값을 가집니다. 즉 반경이 1인 원의 궤적의 1/4에 해당되는 부분입니다. 그러면 그림상에서 원의 안쪽에 해당되는 부분의 면적은 원의 총 면적의 1/4이므로 pi(원주율)의 1/4인 0.7854 정도가 될 것입니다. 그런데 이러한 면적의 값을 Monte Carlo 방법으로 한번 계산해보고자 합니다. 이를 위해서는 위의 그림에서 보이는 XY 공간 내에 임의의 점들을 여러개 뿌려야 합니다. 만약 임의의 위치를 가질 수 있는 100개의 점들을 뿌린다고 가정하면, 확률적으로 원의 안쪽에 위치할 확률이 바깥쪽에 위치할 확률보다 클 것입니다. 이 확률을 계산하면 원의 1/4 면적에 근접하는 값을 얻게 될 것입니다. 물론 뿌리는 점의 갯수가 많을수록 점점 더 실제에 근접한 결과를 얻게 될 것이라 기대할 수 있습니다.

먼저 임의의 위치를 갖는 100개의 점들을 생성하기 위하여 다음과 같이 RANDOMU 함수를 사용하였습니다. RANDOMU 함수는 제가 예전에 이 블로그에서 다룬 적이 있는데(관련게시물 참조), 균일분포(Uniform Distribution)를 갖는 난수들을 생성할 수 있으며, 이렇게 생성된 난수들은 기본적으로 0~1의 범위의 값을 갖게 됩니다. 따라서 xp, yp는 100개의 임의의 점들 각각의 X, Y 좌표값들로 구성된 배열이 됩니다. 당연히 xp와 yp는 0~1의 범위의 값을 임의로 가집니다.

n = 100

xp = RANDOMU(seed, n)

yp = RANDOMU(seed, n)

그런데 이제 임의의 점 각각에 대하여 그 점이 원 궤적의 안쪽에 있느냐 바깥쪽에 있느냐를 판단해야 합니다. 이를 위해서는 원점으로부터의거리에 대하여 반경 1.0을 기준으로 판별하면 됩니다. 즉 각 점의 X, Y 좌표에 대하여 원점으로부터의 거리가 1보다 작거나 같으면 원 궤적 안쪽에 있는 것으로 간주하고, 1보다 크면 바깥쪽에 있는 것으로 간주해 봅시다. 이러한 작업은 다음과 같이 WHERE 함수를 사용하면 간단하게 해결이 됩니다.

ww = WHERE(SQRT(xp^2+yp^2) LE 1, count, COMPLEMENT=wn)

이렇게 하면 n개의 점들 중에서 원점으로부터의 거리가 1보다 작거나 같은 경우에 해당되는 인덱스들이 ww가 되고 그렇지 않은 나머지 경우에 해당되는 인덱스들은 wn이 됩니다. 그리고 count 변수의 값은 거리가 1보다 작거나 같은 경우의 총 갯수가 됩니다. 따라서 count의 값을 총 갯수인 n으로 나눠주면 그 비율이 됩니다. 다음과 같이 이 비율값을 계산하여 실제면적 값과 비교해보면 됩니다.

ratio = FLOAT(count)/n

PRINT, ratio, !pi/4

매번 난수를 발생시켜서 결과를 산출하기 때문에 이 결과는 실행할 때마다 달라지게 됩니다. 제가 5회에 걸쳐서 얻은 결과는 다음과 같습니다. 여러분이 해보시면 ratio 변수의 값은 또 다르게 나올 것입니다.

     0.750000     0.785398

     0.810000     0.785398

     0.810000     0.785398

     0.750000     0.785398

     0.740000     0.785398

그런데 이 값들을 보면 ratio의 값이 실제 면적값인 0.7854에 어느 정도는 근접한 것 같긴 합니다. 사실 n이 100이면 갯수가 아직은 작긴 합니다. 일단 이 시점에서 점들의 위치를 그림상에 중첩하여 표출해 봅시다. 이러한 작업은 다음과 같이 SCATTERPLOT 함수를 사용하면 되는데, 앞서 WHERE 함수로 얻은 결과도 함께 이용해야 원 궤적의 안쪽과 바깥쪽에 속하는 점들을 구분하여 표출할 수 있습니다.

IF count NE 0 THEN BEGIN

  sc_inside = SCATTERPLOT(xp[ww], yp[ww], SYMBOL='circle', /SYM_FILLED, $

    SYM_COLOR='green', SYM_SIZE=0.5, CLIP=0, /OVERPLOT)

  sc_outside = SCATTERPLOT(xp[wn], yp[wn], SYMBOL='circle', /SYM_FILLED, $

    SYM_COLOR='crimson', SYM_SIZE=0.5, CLIP=0, /OVERPLOT)

ENDIF

이렇게 표출을 하면 다음과 같은 그림을 얻게 됩니다. 아무래도 n이 100개이다보니 상당히 썰렁합니다.

그러면 이번에는 n을 100개대신 1000개로 10배 늘려봅시다. 앞의 내용에서 변수 n의 값만 바꿔주면 됩니다. 그러면 출력되는 결과들은 다음과 같습니다. 역시 제가 5회 실행하여 얻은 결과들입니다.

     0.776000     0.785398

     0.778000     0.785398

     0.775000     0.785398

     0.781000     0.785398

     0.789000     0.785398

그리고 n이 1000개일 경우의 그림은 대략 다음과 같습니다.

아까 100개일 때보다는 그래도 좀 북적거리는 느낌이 듭니다. 어쨌든 n을 더 크게 설정하면 좀 더 실제에 근접한 결과를 얻게 될 것입니다. 만약 n을 10000개로 하면 결과값들 및 그림은 다음과 같습니다.

     0.786800     0.785398

     0.781900     0.785398

     0.788400     0.785398

     0.782000     0.785398

     0.785700     0.785398

이 정도면 꽤 볼만한 그림이 나온 것 같습니다. 그리고 이번에는 더 늘려서 n을 100000개로 설정하면 그 결과는 다음과 같습니다.

     0.784690     0.785398

     0.787390     0.785398

     0.789250     0.785398

     0.787890     0.785398

     0.785600     0.785398

이와 같이 이제는 전반적으로 많은 점들로 매우 빼곡하게 채워진 것을 볼 수 있습니다. 어차피 n의 크기는 여러분의 컴퓨터가 감당할 수 있을만큼 크게 해도 됩니다. 그래서 제가 n을 1백만으로 해봤더니 결과값이 0.785558로 산출되었습니다. 0.785398에 더 근접한 결과가 나온 셈입니다. 이와 같이 원의 면적을 Monte Carlo 방법을 사용하여 수치계산적인 방식으로 계산하는 것이 가능합니다. 다음 그림은 n의 값이 1백만일 경우입니다. 점들의 갯수가 워낙 많아서 이제는 빈 틈이 거의 없이 완전히 채워진 것을 볼 수 있습니다.

그런데 이와 유사한 방식으로 면적 대신 원주율의 값을 계산하기도 합니다. 왜냐하면 지금 우리가 보고 있는 것과 같은 반경이 1인 원의 1/4 영역의 면적은 원주율의 1/4이 될 것이기 때문에, 앞서 우리가 계산했던 ratio의 값에 4를 곱하면 원주율에 근접하는 값을 얻을 수 있기 때문입니다. 따라서 다음과 같이 관련 값들을 출력하는 내용을 추가하여 실행하면, 계산된 원주율의 값이 실제 값인 3.141592에 어느 정도 근접하는가를 확인할 수 있습니다. 여기서는 pi_approx 변수의 값을 확인하면 됩니다.

pi_approx = ratio*4

PRINT, n, count, ratio, pi_approx

만약 n이 10000개일 경우에는 결과값들은 다음과 같습니다. 그림은 생략합니다.

   10000        7765     0.776500      3.10600

   10000        7806     0.780600      3.12240

   10000        7799     0.779900      3.11960

   10000        7796     0.779600      3.11840

   10000        7906     0.790600      3.16240

그런데 만약 n이 1백만개가 되면 결과는 대략 다음과 같습니다.

     1000000      785348     0.785348      3.14139

     1000000      785045     0.785045      3.14018

     1000000      785754     0.785754      3.14302

확실히 실제 원주율의 값에 점점 근접하는 것을 확인할 수 있습니다.

따라서 이와 같은 방식으로 결과들을 다양하게 산출해보면 Monte Carlo 방법에 대하여 좀 더 직관적으로 이해하는 것이 가능하지 않을까 생각해 봅니다.