BASICでモンテカルロ法

質問

モンテカルロ法で円周率の推定値を計算することを最近習ったのですが、定積分でもそれが可能なのを知り、どうやってプログラムを組めばいいのか分からず、困っています。
例えば、定積分∫[0→1]x^2dx=1/3～0.333（[0→1]というのは、積分範囲です。）をモンテカルロ法で計算すると、どういうプログラムを組めばいいのでしょうか？
わかる範囲で書いてみたのですが…積分の範囲をどうやってプログラミングすればいいのか、いまいち分かりませんでした。
教えていただけると、助かります。よろしくお願いします。

RANDOMIZE

INPUT n
SET WINDOW -0.1,1.1, -0.1,1.1
DRAW GRID
SET POINT STYLE 1
LET sumin=0

FOR i=0 TO n
LET x=RND
LET y=RND
SET POINT STYLE 2
IF y<x*x THEN
SET POINT COLOR 4
LET sumin=sumin+1
END IF
! PRINT USING "(%.####, %.####)": x,y
PLOT POINTS: x,y
NEXT i

PRINT 1*1*sumin/n
END

回答No.3

cametan_42 · ベストアンサー

ええと、専門用語では「提案分布」と言うんですが、それが積分範囲を「近似しつつ」かつ「計算が簡単な」分布を目的の積分範囲の外側になきゃいけません。
僕はBASICは良く知らないのですが、つまり擬似コード的にはこう言う事です。

1.積分したい関数(例では0≦x≦1の範囲でのx^2)を設定する。
2.上の「被積分関数」を上手く近似しそうな適当な分布を選ぶ。

とは言っても、実際の「確率分布」ではなくって、この場合「0≦x≦1の範囲」ってのをヒントとして、0≦x≦1、0≦y≦1と言う「正方形」を考えます。これが1の「被積分関数」を上手く囲んでいる事が分かるでしょう。

3.(x, y)の乱数を生成して「提案分布」内に乱数をプロットしていく。
4.3.の乱数がy≦x^2を満たしていたら「当たり」、そうじゃなかったら「ハズれ」として、数をカウントしていく。
5.当たりの数/(当たりの数＋ハズレの数)×「提案分布の面積」が目的の被積分関数の面積、になる

こう言うプログラムがオーソドックスな「モンテカルロ法での定積分」の手法となるんではないでしょうか?

お礼

ありがとうございました。参考にさせていただきました。

回答No.4

cametan_42

僕はBASICを全く知らないので、代わりにSchemeと言う言語で問題のモンテカルロ法のソースを書いてみました。
何かの参考にして頂ければ幸いです。

(define (montecarlo n)
(let loop ((i 0) (j 0));カウンターi、成功数jの初期値
(if (= i n) ;終了条件
(exact->inexact (/ j i)) ;解＝成功数/カウント
(let ((x (random)) (y (random)));x,yを設定
(loop (+ i 1) (+ j (if (<= y (* x x));繰り返し計算
1
0)))))))

これをScheme処理系(例えばフリーウェアのDrScheme)で走らせると以下のような結果になりました。

> (montecarlo 1000000)
0.333365

と言うわけで、目的は達成できたとは思います。

お礼

ありがとうございました。参考にさせて頂きました。

回答No.2

asuncion

0～1の範囲の乱数を2つ求めて、それらをx, yとします。
y <= x*x を満たすとき、個数（例示のコードではsumin）を
１つ増やします。

乱数の組（x, y）を求める回数を十分大きく取れば、
sumin / 回数
の値は0.33333... に近づいていくことでしょう。

お礼

ありがとうございました。参考にさせて頂きました。

回答No.1

Hayashi_Trek

提示しているプログラムは、モンテカルロ法による円周率の計算のようですね。
「定積分」ということは面積です。
半径１の円を考えます。（式１：x^2 + y^2 = 1）
計算を簡単にするために、x,yが共に正の部分（第一象限）のみに注目すると
1/4円の面積は、π/4となります。（理論値）
1/4円の面積を求めるために、ｘ軸を10等分して
円に収まる棒グラフのような四角形を考えます。
各棒グラフの高さは、式１より
y = √（1 - x^2）
各棒グラフの面積を合計すれば、円の面積の近似値が得られます。

あとはｘ軸の分割数を大きくして、より詳しく計算するだけです。
各棒グラフを単純な四角形から台形にすればさらに正確な近似値が得られます。

お礼

ありがとうございました。参考にさせて頂きました。