裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

π を明示的には使わない,ビュフォンの針のシミュレーション(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))

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ややこしいプログラム -- R の eval() と quote()

2020年11月23日 | ブログラミング

R の eval() は,Python の eval() とは全く異なるものである。ということを書きたい。

以下の例示プログラムを使って説明しよう。

 1:  func = function(x = NA, y = NA) {
 2:    z = function(f0) {
 3:      return(f0(10))
 4:    }
 5:    f.body = quote(2 * x + 3 * y + 5)
 6:    if (is.na(x)) {
 7:      x = z(function(x) eval(f.body))
 8:      print(x)
 9:    } else if (is.na(y)) {
10:      y = z(function(y) eval(f.body))
11:      print(y)
12:    }
13:  }

このプログラムは

func(x = 1)
func(y = 2)

のように利用される。

func(x = 1) の場合,y = NA なので,9 行目の条件判定が成立し,10 行目の y = z(function(y) eval(f.body)) が実行される。
eval(f.body) は,5 行目で定義されている f.body の右辺の quote のなかみと置き換えられる。つまり,y = z(function(y) 2 * x + 3 * y + 5) と書かれるのと同じ効果を生む。
z は 2〜4 行目で定義されている関数で,引数は一つだけで,その引数は別の関数である。
10 行目は,無名関数 function(y) 2 * x + 3 * y + 5 を引数として z を呼ぶことになる。この無名関数の変数は y であり,x は呼ばれたときに x が持っている定数である。
関数 z は 仮引数の関数 f0 の変数に 10 を代入した結果を返す(実際の z はもっと複雑なことをするのだろうが)。
f0(y) = 2 * x + 3 * y + 5 で,f0(10) を求めるということである。x の値は 1 なので,f0(10) = 2 * 1 + 3*10 + 5 = 37 が表示される。

func(y = 1) の場合,x = NA なので,6 行目の条件判定が成立し,7 行目の x = z(function(x) eval(f.body)) が実行される。
これは前述の場合と同じように,x = z(function(x) 2 * x + 3 * y + 5) と書かれるのと同じで,無名関数 function(x) 2 * x + 3 * y + 5 を引数として z を呼ぶ。この無名関数の変数は x であり,y は定数である。
関数 z は f0(x) = 2 * x + 3 * y + 5 で,y の値は 2 であり,f0(10) = 2 * 10 + 3 * 2 + 5 = 31 が返される。

つまり,z に渡される関数の本体はいずれの場合にも 2 * x + 3 * y + 5 であるが,本体中に現れる x, y のどちらが変数であるか,function(x) なのか function(y) なのかが異なる。関数本体が 5 行目で表されるような簡単なものならば,7, 10 行目に直接書き込んでも良いのだが,長くて複雑な場合や,関数本体に複数の変数を含むような場合にはそれぞれに対応して 8,10 行目以外にも複数箇所で同じ関数本体を重複して書かねばならないことになる。それを避けるための手段が eval() と quote() を使うということである。

さて,長くなったが,これを Python で書くと

 1:  import numpy as np
 2:  from numpy import nan as NA
 3:  
 4:  def func(x = NA, y = NA):
 5:      def z(f0):
 6:          return f0(10)
 7:      if np.isnan(x):
 8:          def f_x(x):
 9:              return 2 * x + 3 * y + 5
10:          x = z(f_x)
11:          print(x)
12:      elif np.isnan(y):
13:          def f_y(y):
14:              return 2 * x + 3 * y + 5
15:          y = z(f_y)
16:          print(y)

となり,

func(x = 1)
func(y = 2)

で呼び出すと,R のプログラムと同じ働きをする。
R のような eval() と quote() の対がないので,関数本体が同じ(9, 14 行目)だが,どれが変数かを示すところが違う(8,13 行目)関数定義を2箇所に書かねばならない。
関数が短い場合は,以下のように lambda 記法で書けば行数が少なくて済むこともあるが,同じようなことを複数回書かねばならないことに違いはない。

    if np.isnan(x):
        x = z(lambda x: 2 * x + 3 * y + 5)
        print(x)
    elif np.isnan(y):
        y = z(lambda y: 2 * x + 3 * y + 5)
        print(y)

ところで,Python にも eval() はある。文字列で表されている式を評価するものである。
body = '2 * x + 3 * y + 5'
で,x = 10; y = 2 のとき,eval(body)  は 2 * x + 3 * y + 5 に x, y を代入して,式を評価し,結果 37 を返す。
x = 1; y = 10 のとき,eval(body)  は 2 * x + 3 * y + 5 に x, y を代入して,式を評価し,結果 31 を返す。
なんだ,同じようなことをやってるじゃないか?とお思いかもしれないが,以下のようなプログラムは動かない。

 1:  import numpy as np
 2:  from numpy import nan as NA
 3:  
 4:  def func(x = NA, y = NA):
 5:      def z(f0):
 6:          return f0(10)
 7:      body = "2 * x + 3 * y + 5"
 8:      if np.isnan(x):
 9:          x = z(lambda x: eval(body))
10:          print(x)
11:      elif np.isnan(y):
12:          y = z(lambda y: eval(body))
13:          print(y)

2 箇所の eval(body) は,x, y いずれかが nan なので,eval(body) は nan になる。

というように,R の eval() は Python の eval とは全くの別物ということである。

Python に R の eval() のようなものはないのだろうか??

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村