モンテカルロ法で円周率を求める

import random
from re import I
N = 1000000
p_in = 0
p_out = 0
for _ in range(N):
    rx = random.random()
    ry = random.random()
    if rx**2+ry**2 <=1:
        # inside
        p_in += 1
    else:
        p_out += 1

all = p_in+p_out
target_area = p_in/all
circumference = target_area*4

print(circumference)

実行結果
3.142424

違う実装(座標値を配列に追加しているので遅い)

import random
import matplotlib.pyplot as plt
import numpy as np

N = 100000
p_in_x = np.array([])
p_in_y = np.array([])
p_out_x = np.array([])
p_out_y = np.array([])
for _ in range(N):
    rx = random.random()
    ry = random.random()
    if rx**2+ry**2 <=1:
        # inside
        p_in_x = np.append(p_in_x,rx)
        p_in_y = np.append(p_in_y,ry)
    else:
        p_out_x = np.append(p_out_x,rx)
        p_out_y = np.append(p_out_y,ry)

all = len(p_in_x)+len(p_out_y)
target_area = len(p_in_x)/all
circumference = target_area*4

print(circumference)
plt.plot(p_in_x,p_in_y,'bo')
plt.plot(p_out_x,p_out_y,'ro')
plt.axis('equal')
plt.show()

実行結果
3.14348
f:id:seinzumtode:20220117140842p:plain