モンテカルロ法を用いた積分計算のプログラミング

今、指数関数の積分をモンテカルロ法を用いてMatlab上で計算しようとしています。

途中までプログラムを組みましたが、間違っている個所があるようなのです。
rand()の使い方が間違っているのでしょうか？

また、rand()は一様分布に従う疑似乱数を返すという解釈で合ってますでしょうか？
プログラム上にあるX2がとる値の大きさと、randのとる値の大きさは？と聞かれたのですが、
どのように答えていいか分りません。。。

下にプログラムを載せました。
よろしくお願いいたします。

N = 10000; %the number of calculation
a =0; %begin of interval
b =1; %end of interval
c = exp(b); %constant of max
f =0;

for n = 1:N;
X1 = a +rand()*(b-a);
%Generate values from the uniform distribution on the interval [a, b]
X2 = a + rand()*(1-0);
if (c * X2 ) < (exp(X1))
f = f +1; %count the number of point in the area
end

end
y = f*c*(b-a) / N % f/N:ratio of points, c*(b-a):area of integral
Int = exp(b) - exp(a)
err = abs(Int - y)

No.5 ベストアンサー 回答者： Kules 回答日時： 2011/05/05 13:57

A No.2のＫｕｌｅｓです。

Ｔａｃｏｓａｎさんの回答で思い出しました。
そうか、長方形内部の点をランダムにうつからX1,X2は独立した乱数でないといけないですね。
で、今回は(a,0),(b,0),(a,exp(a)),(b,exp(b))を頂点に持つ長方形内にランダムにN点(X1,X2)を打って、
その中でX2<exp(X1)になってるものの数を数える、ということでしょうか。

とりあえずMatlab初心者ということなのでアドバイスのようなものを少し
（これも質問に対する回答からは少しずれていると思いますが）

・配列で何とかできることは配列の形でしてしまう
C言語などではスカラー量しか扱えない関数も、Matlabではベクトルの形で入力することができます。
randでいうと、
http://www.mathworks.co.jp/help/ja_JP/techdoc/re …
を見てもわかるとおり、rand(m,n)でm行n列の擬似乱数値の行列を作れます。
よって、今回はforループでN回randを呼び出さなくても、
X1=a+rand(N,1)*(b-a);
でN行1列の一様分布擬似乱数が得られます。

・ｉｆ文はlogicalを使って書き変えられないか考える
例えば
A=rand(5,1);
B=A<0.5;
とすると、例えばBは
[0;1;0;1;1]のように0or1で表される配列になります。
(B(k)=1(A(k)<0.5),0（else）というようになります)
ということは、これを使えば
sum(B)とすると、「A<0.5となる要素の数」を表すことができます。
（sum(B)は、Bの合計を出します）
このようにすることで、わざわざif文で１つ１つの要素を調べなくても、
一度に調べることができます。


以上より、例えば今回であれば
X1=a+rand(N,1)*(b-a);
X2=0+rand(N,1)*(c-0);
f=sum(X2<exp(X1));
とすれば、forループ部分で求めるべきfは求まります。

まあよっぽど非力なPCを使ってない限り双方の書き方に時間的な差は出ないと思いますし、
Nが極端に大きくなると配列を用意するメモリが足りなくなるので万能ではありませんが。


参考になれば幸いです。

この回答への補足

ご回答、ありがとうございます。

アドバイスの方も、今後の参考にさせていただきたいです。

もしよろしければお教えいただきたいのですが、
X1=a+rand()*(b-a)
で一様分布を発生させるのと、
rand(,)で発生させるのとでは、何が異なるのでしょうか？

すみませんがよろしくお願いいたします。

補足日時：2011/05/18 21:37

この回答へのお礼

すみません、別問題なんですが、

I = int_0^1 int_0^1-x exp(-x)*exp(-y)

の二重積分を考えております。
この場合ですと、前回と同様に、

X = a + rand()*(1-0);
Y = a + rand()*(1-X-a);
と一様分布を発生させて、

X2 = 0 + rand()*(exp(-1));
if ( X2 ) < (exp(-Y))

と言うように比較して計算できないかと考えたのですが、
こちらでも大丈夫でしょうか？
p=rand(1,2)とし、
P(1)+p(2) <= 1
のとき、
f = exp(-p(1))*exp(-p(-2))
を計算し、
S = f/n(試行回数)
で積分値を得る、とも考えました。

ご教示いただけますでしょうか。よろしくお願いいたします。

お礼日時：2011/05/18 21:51

No.7 回答者： Kules 回答日時： 2011/05/18 23:28

A No.2,5のKulesです。

まず、
>X1=a+rand()*(b-a)
>で一様分布を発生させるのと、
>rand(,)で発生させるのとでは、何が異なるのでしょうか？

について。
結論から言うと、出てくる結果は同じになります。
http://www.mathworks.co.jp/help/ja_JP/techdoc/re …
で書かれているようなものを使って乱数ストリームを制御し、同じ乱数が出力されるようにすれば、
for n=1:N;
X1=a+rand()*(b-a);
end;
でX1に出てきたものを記録したものと、
X1=a+rand(N,1)*(b-a);
として出力されるものは全く同じになります。
要はN個の乱数を作るのに
「マジメに乱数を1個ずつN回作る」か「N個一気に作って、まとめて計算するか」の違いなので、
そこまでNが大きくなければ好みの問題です。
私は配列で一気に作って一気に計算するのが好きなので後者を選択します。

また、
>I = int_0^1 int_0^1-x exp(-x)*exp(-y)
についてですが、
exp(-x)*exp(-y)をxは0から1まで、yは0から1-xまで積分ってことでいいですかね？
正直重積分まで来ると私の想像力では太刀打ちできないですね…
役立たずで申し訳ないです。

参考になれば幸いです。

この回答へのお礼

ご回答、ありがとうございました！
無事、プログラムを組み終える事が出来ました。

Kulesさんのご親切に感謝いたします。

お礼日時：2011/06/02 04:06

No.6 回答者： Tacosan 回答日時： 2011/05/07 21:52

別法としては

average(exp(a+rand(N, 1)*(b-a)))*(b-a)
のような感じでもできそう.

この回答への補足

何度もご親切にありがとうございます。
大変、嬉しく思います。

新しく問題を立てるのが適切かと思いましたが、こちらに失礼します。

今、別問題で、

I = int_0^1 int_0^1-x exp(-x)*exp(-y)

の二重積分を考えております。
この場合ですと、前回と同様に、

X = a + rand()*(1-0);
Y = a + rand()*(1-X-a);
と一様分布を発生させて、

X2 = 0 + rand()*(exp(-1));
if ( X2 ) < (exp(-Y))

と言うように比較して計算できないかと考えたのですが、
間違っていますでしょうか。

よろしくお願いいたします。

補足日時：2011/05/18 21:47

No.4 回答者： Tacosan 回答日時： 2011/05/05 12:35

「モンテカルロ積分」をどこまで理解できているのか不安.

それは「戦略」じゃなくてせいぜい「戦術」だと思う... 「戦略」としては
対象領域 D を含んだ長方形領域 R の内部に (独立かつ) 一様ランダムに点を取ると「(D の内部にある点の数)/(生成した点の数)」という比が「(D の面積)/(R の面積)」という比に一致する
ということだね. だから (X1, X2) が考えている長方形領域 R の内部にあるように X1, X2 を生成するのが正しい.

で, 結論としては #1 で正解... というか, 上の戦略を忠実に記述するなら (今の場合 R の縦方向は 0 から exp(b) なので) X2 は「0 から exp(b) までの一様乱数」とすべき. 「0 から 1 までの一様乱数」に c = exp(b) を掛けてもいいけど, ちょっと回りくどい.

この回答への補足

ご回答、ありがとうございます。

Monte Carlo積分に着いて今、色々と文献を参考にしながら勉強している途中でして、
お察しの通り、まだ理解できていないです。

X2を０からexp(b) までの一様乱数として、下の部分もそれに合わせて組み直し、
どうにか出来ました。
ありがとうございます。

補足日時：2011/05/18 21:08

No.3 回答者： Tacosan 回答日時： 2011/05/05 01:38

「モンテカルロ法で積分」といってもいくつかの戦略が考えられます. プログラムを示すだけでなく, あなたのとった戦略を言葉で説明して

ください. もちろん, その中で X1 や X2 がどのような意味を持つのかも明らかにしてください.

この回答への補足

回答、ありがとうございます。

exp(x)を区間aからbまでの積分をします。
X1は区間ab内に一様分布に従う乱数を発生させ、それをexp()に代入します。

もうひとつ乱数を発生させ、それとexp(b)とをかけたものが先ほど求めたexp(x1)よりも小さければその点をカウントをします。

そして、最後に、全部の総面積を求める為、
区間ａｂの長さを底辺、高さをexp(b)、その中に入る点の数の割合を(カウントされた点の和/全体の試行数N)とし、その積を面積としました。

説明が下手で申し訳ありませんがよろしくお願いいたします。

補足日時：2011/05/05 02:04

No.2 回答者： Kules 回答日時： 2011/05/05 00:27

ん～何かＭａｔlabっぽくない感じのプログラムですね。

そもそもこれってループ使わなくても書けるような。
まあこれは本筋じゃないのでおいといて。

とりあえずわからないのは、
X1 = a +rand()*(b-a);
は閉区間[a,b]内の点をランダムに選んでるとして、
X2は何をしているんですかね？
if文にある (c * X2 ) < (exp(X1))
を見るとy=exp(x)上の点（X1,exp(X1)）と
y=c*x上の点(X2,c*X2)の大小関係を見てるみたいなんですが、
それだとX1とX2に同じ乱数の値を使わないといけない気がするんですが。

X2の意味と、
X2 = a + rand()*(1-0);
の1,0辺りの意味合いを補足いただければもう少し有益な回答ができるかも知れません。
（私の不勉強で申し訳ないですが）

以上、参考になれば幸いです。

この回答への補足

ご回答、ありがとうございます。

プログラミングの経験も浅く、Matlabをこの春から勉強しだしました。

X2はX1よりも上にあるか下にあるか、という比較の為に使いました。

(c * X2 ) < (exp(X1)) から、
Cはexp(b)であるため、C*X2を使うと、
X2は区間[0,1]内にないとX1を越えてしまうのではないかと思い、
X2 = a + rand()*(1-0);
と書きました。

すみませんがよろしくお願いいたします。

補足日時：2011/05/05 00:49

No.1 回答者： himajin100000 回答日時： 2011/05/04 22:49

MatLabもモンテカルロ法も使ったことないけど、たぶん

>X2 = a + rand()*(1-0);
が
>X2 = 0 + rand()*(1-0);
の書き間違いで、たまたまa = 0だから、うまくいくんじゃないかと思うんだけど…

#Wikipediaのモンテカルロ法の説明見る限り、最大値が予め分かってないといけないんだよね？多分。