裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

π を明示的には使わない,ビュフォンの針のシミュレーション(Python)

2020年11月23日 | Python

ビュフォンの針のシミュレーションプログラムはたくさん書かれているが,プログラム内で定数 π を使っているものがほとんどである。

π を使わないシミュレーションプログラムを書いたので記録に残しておく。

方針は,原点を左下の頂点とする一辺  1 の正方形内のランダムな点の座標(x2,y2) のうち,原点からの距離が 1 以下の点について sin(theta) は三角関数の定義により y2/sqrt(x2**2 + y2**2) であることを使う。

以下のプログラムで DEBUG = True として実行すると,theta の分布は一様分布になっていることが確認できる。

 

なお,プログラム中の REMARK と書かれた if 文を外すと,半径が 1 より大きい点も採用されてしまうので,theta は一様分布にはならない。

import random
import math

DEBUG = False
if DEBUG:
    import numpy as np
    import matplotlib.pyplot as plt

n = 100000
d = 10
l = d / 2
i = 0
cross = 0
if DEBUG: theta = np.zeros(n)
while i < n:
  y = random.uniform(0, d / 2) # 針の中心から最も近い平行線までの距離
  # 0 <= theta <= π/2 の一様分布を求める
  x2 = random.uniform(0, 1) # 一辺 1 の正方形内の座標 (x2, y2)
  y2 = random.uniform(0, 1)
  r2 = x2**2 + y2**2 # 原点から (x2, y2) までの平方距離
  if r2 < 1: ## REMARK 円の内部の点に限定する
    r = math.sqrt(r2) # 原点から (x2, y2) までの距離
    sine_theta = y2 / r # sine の定義により
    if DEBUG: theta[i] = math.asin(sine_theta)
    i += 1
    if y <= l / 2 * sine_theta: # 平行線と交わる
      cross += 1

if DEBUG: plt.hist(theta, bins=np.arange(0, np.pi/2, np.pi/60))
print(2 * l * n / (cross * d))

 

コメント   この記事についてブログを書く
« ややこしいプログラム -- R ... | トップ | Python で直線回帰 »
最新の画像もっと見る

コメントを投稿

ブログ作成者から承認されるまでコメントは反映されません。

Python」カテゴリの最新記事