[Python & SQL] 몬테 카를로 | 원주율 구하기

안녕하세요.

 

이승혁입니다.

 

오늘은 몬테 카를로 방법(Monte Carlo method)을 사용해 원주율(3.14......)을 구해보도록 하겠습니다.

 

### 몬테 카를로 방법(Monte Carlo method) ? 

 

난수를 생성해 함수의 값을 확률적으로 계산하는 방법

 

계산하려는 값이 닫힌 형식으로 표현되지 않거나 복잡한 경우에 근사적으로 계산할 때 사용됩니다.

 

원주율과 같이 딱 하나의 값으로 수렴하지 않는 경우에 사용합니다.

 

 

## 원 그리기

 

반지름 r, 중심이 (a, b) 인 원의 방정식은 (x - a)2 + (y - b)2 = r2 입니다.

 

반지름 1 , 중심 (0,0) 인 원의 식은 어떻게 될까요 ?

 

x^2 + y^2 = 1

 

파이썬으로 그려서 확인해 보도록 하겠습니다.

 

import matplotlib.pyplot as plt

circle=plt.Circle((0,0),1,fc='w',ec='b') # plt의  Cricle 사용 중심, 반지름, 면적 색, 테두리색
a=plt.axes(xlim=(-1,1),ylim=(-1,1)) # x축 범위 , y축 범위를 a 변수에 할당

a.add_patch(circle) # 만들어진 a 변수(범위)에 circle patch
a.set_aspect('equal') # x ,y 비율은 같게

plt.grid(True, axis='both',color='black', alpha=1, linestyle='-') # 격자 표시
plt.show()    

 

 

 

이렇게 중심이 0,0 이고 반지름이 1인 원을 그려 보았습니다.

 

우리가 사용할 부분은 2사분면 입니다.

 

# 직사각형의 꼭지점 좌표

 

1. 0,0

2. 1,0

3. 0,1

4. 1.1

 

그 중 원의 테두리에 걸치는 점은 2번과 3번입니다.

 

x^2 + y^2 = 1 , 원의 방정식에 맞는 모습이죠?

 

 

### 원주율 (3.14...) 도출하기

 

그럼 이제 3.14를 구하기 위해서는 어떻게 해야 할까요 ?

 

x 좌표 y좌표 값을 각각 제곱해 더합니다 ( x^2 + y^2 )

 

그 값이 1이면 원의 테두리 입니다 ( x^2 + y^2 = 1 )

 

그럼 그보다 같거나 작으면 ? ( x^2 + y^2 <= 1 )

 

원의 안쪽 입니다. 

 

이때의 x , y 좌표의 개수를 확인해 보겠습니다

 

### Python 원주율 구하기

 

1. 100개의 난수

 

import numpy as np
import matplotlib.pyplot as plt

N=100
x=np.random.rand(N)
y=np.random.rand(N)

var_x=[i for i,j in zip(x,y) if i**2+j**2<1]
var_y=[j for i,j in zip(x,y) if i**2+j**2<1]
    
print('%s 개의 점 중 %s개'%(N,len(var_x)))

c=plt.Circle((0,0),1,fc='w',ec='b')
a=plt.axes(xlim=(-1,1),ylim=(-1,1))
a.add_patch(c)
a.set_aspect('equal')
plt.grid(True, axis='both',color='black', alpha=1, linestyle='-')
plt.plot(var_x,var_y)
plt.show()    

 

보시는 바와 같이 2사분면의 안쪽에 찍혔고, 100개 난수 중 79개가 해당됩니다.

난수의 개수를 늘려서 개수를 확인하겠습니다.

 

2. 난수 100000

 

 

 

10000개 중 7810개로 안쪽을 빼곡하게 채웠습니다.

 

그럼 이때의 넓이는 어떻게 될까요 ?

 

1 x 1 x 3.14 / 4 = 0.785 입니다.

 

안쪽에 찍힌 점은 7810개 입니다.

 

뭔가 감이 오시나요 ?

 

7810개에 해당하는 개수는 부채꼴(1/4)입니다.

 

따라서 전체 개수 중 해당하는 비율에 4를 곱하게 되면 원주율이 나오게 됩니다.

 

import numpy as np
import matplotlib.pyplot as plt

N=10000
x=np.random.rand(N)
y=np.random.rand(N)

var_x=[i for i,j in zip(x,y) if i**2+j**2<1]
var_y=[j for i,j in zip(x,y) if i**2+j**2<1]
    
print('%s 개의 점 중 %s개'%(N,len(var_x)))
print('%s / %s * 4 = %s'%(len(var_x), N ,len(var_x) / N * 4))

c=plt.Circle((0,0),1,fc='w',ec='b')
a=plt.axes(xlim=(-1,1),ylim=(-1,1))
a.add_patch(c)
a.set_aspect('equal')
plt.grid(True, axis='both',color='black', alpha=1, linestyle='-')
plt.rcParams["figure.figsize"] = (10,10)
plt.plot(var_x,var_y)

plt.show()    

 

 

### SQL 원주율 구하기

 

SQL을 통해서도 원주율을 구해보았습니다.

0과 1사이의 난수를 10만개 생성해 제곱한 결과를 저장합니다.

제곱해서 더한 결과가 1보다 작은 개수를 카운트합니다.

 

전체 10만개 중 개수의 비율에 곱하기 4를 하여 원주율을 구합니다.

select (count(*)/100000)*4 pi
  from(
      select  (power(dbms_random.value(0,1),2) + power(dbms_random.value(0,1),2)) as "x^2+y^2"
       from dual
       connect by level <= 100000
     )
  where "x^2+y^2"<=1;

 

 

 

오늘은 몬테카를로 방법을 이용해서 원주율을 구해 보았습니다.

 

또 다시 프로그래밍에 있어 수학의 필요성을 느끼게 됩니다.

 

원주율을 구하고 시각화 하는 과정을 진행하며

파이썬으로 원 그리기, 난수 생성하기, 그래프 겹치기, 그리드 설정하기 등 

추가적인 공부가 더 되어서 좋은 실습이었습니다.

 

많은 도움이 되었길 바라면서 마치겠습니다.

 

읽어주셔서 감사합니다!