numpy를 활용해 원주율을 구하는 Mote Carlo Simulatino을 시각화 합니다.
numpy로 원주율값 구하기
이미지 출처 : https://byun1114.tistory.com/1071
-
충분히 많은 점을 랜덤으로 찍을때, 네모의 1/4 넓이 안에 들어갈 확률은 4분원의 크기를 네모의 크기로 나눈 값과 같다.
-
이를 활용해 원주율값을 구할 수 있다.
1
import numpy as np
-
원의 넓이는 1 * 1 이므로 1
-
4분원의 넓이는 원의 넓이인 𝝅 * R의 제곱을 4로 나눈 값이다.
-
즉 𝝅/4 가 된다. (현재 네모에서 R은 1)
-
원의 넓이 : 네모의 넓이
pi / 4 : 1 = (4분원 안에 생성된 점의 개수) : 전체 (점을 찍는)시도 횟수
pi = 4 * (4분원 안에 생성된 점의 개수) / 전체 시도 횟수
랜덤으로 점을 찍는 총 횟수 표현
1
2
3
4
5
6
7
8
# 전체 시도 횟수를 천만번 (10의 7제곱)이라 가정
total_trial = int(1e7)
# 점을 찍은 횟수
# 0 ~ 1 사이의 값을 가져오기 위해 rand사용
points = np.random.rand(total_trial,2)
points
array([[0.49354421, 0.75833051], [0.86829397, 0.92375882], [0.72372899, 0.3364979 ], ..., [0.27542163, 0.9882378 ], [0.06996705, 0.31490832], [0.55447043, 0.12300586]])
1
2
3
# 행은 1000만개, 열은 2개가 있는 행렬 생성
# 열은 x,y 좌표, 10000만개는 각각의 점을 의미
points.shape
(10000000, 2)
4분원의 넓이 구하기
-
위의 그림에서 파란선 안쪽의 범위
-
x^2 + y^2이 >= 1이면 원안의 좌표이다.
1
2
3
4
# 각 점의 제곱한 값
points ** 2
# 여기서 제곱한 값에서 x + y를 계산해야 한다. -> np.sum활용
array([[0.24358589, 0.57506517], [0.75393442, 0.85333036], [0.52378365, 0.11323084], ..., [0.07585707, 0.97661396], [0.00489539, 0.09916725], [0.30743746, 0.01513044]])
1
2
# 행렬에서 '열' 을 따라 계산 (axis=1)
np.sum(points ** 2, axis=1)
array([0.81865105, 1.60726477, 0.63701449, ..., 1.05247103, 0.10406264, 0.3225679 ])
위의 값이 1보다 작으면 원 안에 들어가는 것
1
2
# x+y가 1보다 작은지 필터링
np.sum(points ** 2, axis=1) < 1
array([ True, False, True, ..., False, True, True])
-
np.sum(points ** 2, axis=1) 의 결과는 벡터
-
~ <1 의 결과는 scalar 이다.
-
따라서 numpy broadcasting이 적용되면서 boolean이 된다.
- 여기서 최종적으로 알고자 하는 것은 1이상인 수의 ‘개수’ 이다. -> np.sum() 활용 (True는 1인 것을 활용)
1
2
# 천만개의 점중에서 4분원 안에 들어간 점의 개수
np.sum(np.sum(points ** 2, axis=1) < 1)
7851866
pi = 4 * (4분원 안에 생성된 점의 개수) / 전체 시도 횟수
1
4 * np.sum(np.sum(points ** 2, axis=1) < 1) / total_trial
3.1407464
완벽하게 pi의 값이 정확하게 나오지는 않지만 근사치의 값은 구할 수 있다.
댓글남기기