今、指数関数の積分をモンテカルロ法を用いて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.1
- 回答日時:
MatLabもモンテカルロ法も使ったことないけど、たぶん
>X2 = a + rand()*(1-0);
が
>X2 = 0 + rand()*(1-0);
の書き間違いで、たまたまa = 0だから、うまくいくんじゃないかと思うんだけど…
#Wikipediaのモンテカルロ法の説明見る限り、最大値が予め分かってないといけないんだよね?多分。
No.2
- 回答日時:
ん~何かMatlabっぽくない感じのプログラムですね。
そもそもこれってループ使わなくても書けるような。まあこれは本筋じゃないのでおいといて。
とりあえずわからないのは、
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);
と書きました。
すみませんがよろしくお願いいたします。
No.3
- 回答日時:
「モンテカルロ法で積分」といってもいくつかの戦略が考えられます. プログラムを示すだけでなく, あなたのとった戦略を言葉で説明して
ください. もちろん, その中で X1 や X2 がどのような意味を持つのかも明らかにしてください.この回答への補足
回答、ありがとうございます。
exp(x)を区間aからbまでの積分をします。
X1は区間ab内に一様分布に従う乱数を発生させ、それをexp()に代入します。
もうひとつ乱数を発生させ、それとexp(b)とをかけたものが先ほど求めたexp(x1)よりも小さければその点をカウントをします。
そして、最後に、全部の総面積を求める為、
区間abの長さを底辺、高さをexp(b)、その中に入る点の数の割合を(カウントされた点の和/全体の試行数N)とし、その積を面積としました。
説明が下手で申し訳ありませんがよろしくお願いいたします。
No.4
- 回答日時:
「モンテカルロ積分」をどこまで理解できているのか不安.
それは「戦略」じゃなくてせいぜい「戦術」だと思う... 「戦略」としては
対象領域 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を0からexp(b) までの一様乱数として、下の部分もそれに合わせて組み直し、
どうにか出来ました。
ありがとうございます。
No.5ベストアンサー
- 回答日時:
A No.2のKulesです。
Tacosanさんの回答で思い出しました。
そうか、長方形内部の点をランダムにうつから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列の一様分布擬似乱数が得られます。
・if文は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文で1つ1つの要素を調べなくても、
一度に調べることができます。
以上より、例えば今回であれば
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(,)で発生させるのとでは、何が異なるのでしょうか?
すみませんがよろしくお願いいたします。
すみません、別問題なんですが、
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(試行回数)
で積分値を得る、とも考えました。
ご教示いただけますでしょうか。よろしくお願いいたします。
No.6
- 回答日時:
別法としては
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))
と言うように比較して計算できないかと考えたのですが、
間違っていますでしょうか。
よろしくお願いいたします。
No.7
- 回答日時:
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まで積分ってことでいいですかね?
正直重積分まで来ると私の想像力では太刀打ちできないですね…
役立たずで申し訳ないです。
参考になれば幸いです。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
関連するカテゴリからQ&Aを探す
おすすめ情報
- ・漫画をレンタルでお得に読める!
- ・人生のプチ美学を教えてください!!
- ・10秒目をつむったら…
- ・あなたの習慣について教えてください!!
- ・牛、豚、鶏、どれか一つ食べられなくなるとしたら?
- ・【大喜利】【投稿~9/18】 おとぎ話『桃太郎』の知られざるエピソード
- ・街中で見かけて「グッときた人」の思い出
- ・「一気に最後まで読んだ」本、教えて下さい!
- ・幼稚園時代「何組」でしたか?
- ・激凹みから立ち直る方法
- ・1つだけ過去を変えられるとしたら?
- ・【あるあるbot連動企画】あるあるbotに投稿したけど採用されなかったあるある募集
- ・【あるあるbot連動企画】フォロワー20万人のアカウントであなたのあるあるを披露してみませんか?
- ・映画のエンドロール観る派?観ない派?
- ・海外旅行から帰ってきたら、まず何を食べる?
- ・誕生日にもらった意外なもの
- ・天使と悪魔選手権
- ・ちょっと先の未来クイズ第2問
- ・【大喜利】【投稿~9/7】 ロボットの住む世界で流行ってる罰ゲームとは?
- ・推しミネラルウォーターはありますか?
- ・都道府県穴埋めゲーム
- ・この人頭いいなと思ったエピソード
- ・準・究極の選択
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
配列数式の解除
-
2つ以上の変数を比較して最大数...
-
VBA 1次元配列を2次元に追加する
-
for each の現在の配列ポインタ...
-
【MATLAB】任意の多次元配列か...
-
VBAで近似曲線の係数取得
-
subの配列引数をoptionalで使う...
-
ExcelVBA Rangeの扱いについて
-
配列に同じ値を入れる方法
-
2次元動的配列の第一引数のみを...
-
AES暗号にて、AES_set_encrypt_...
-
linest関数に配列を渡す
-
CGIでカスタム配列でソート
-
配列を任意の数値で埋める方法
-
Array配列の末尾に追加したい。
-
配列数式って何ですか??
-
ファイルの列挙
-
Excel VBA配列をFunctionに渡す
-
[Excel2000_VBA] 型が一致しま...
-
VB6 配列を初期化したい
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
配列数式の解除
-
2つ以上の変数を比較して最大数...
-
VBA 1次元配列を2次元に追加する
-
VB6 配列を初期化したい
-
特定のセル範囲で4文字以上入力...
-
ListViewで、非表示列って作れ...
-
《エクセル2000》A列・B列の商...
-
配列変数の添字が範囲外ですと...
-
Excel-VBAの配列「Public Const...
-
subの配列引数をoptionalで使う...
-
for each の現在の配列ポインタ...
-
配列を任意の数値で埋める方法
-
Dim は何の略ですか?
-
VBのFunctionで、配列を引数...
-
配列内の内容を全て表示する方法
-
2次元動的配列の第一引数のみを...
-
Excel VBA配列をFunctionに渡す
-
VBA Match関数の限界
-
Array配列の末尾に追加したい。
-
AES暗号にて、AES_set_encrypt_...
おすすめ情報