• ベストアンサー

凖モンテカルロ法 or モンテカルロ + MT

3次元に広がる物理量を面積分する方法について検討中です。 凖モンテカルロという方法でLDSという数列を用いて計算をすると、積分計算の 誤差の軽減が速まるらしいということがわかりました。原理はまだ勉強中ですが、 対象の空間を区切って、それぞれの細空間に乱数のポイントを配置することが いい方向に働くように見られます。 これに対して、通常のモンテカルロ法でメルセンヌツイスター(MT)を用いた場合との 違いがいまいちわかりません。MTで一様乱数ができるのであれば、凖モンテカルロ のようにする必要はないのでしょうか?

noname#40366
noname#40366

質問者が選んだベストアンサー

  • ベストアンサー
  • goma_2000
  • ベストアンサー率48% (62/129)
回答No.1

乱数を用いた計算はどちらも同じですが、乱数の収束性が違います。 たとえば、一様といったときに、200点くらいでは一様になりませんよね。ではどのくらいの点数があれば一様と見做せるのか？が準乱数とMTなどの乱数の違いです。これが誤差の軽減に効くということです。 ちなみに、準乱数は乱数と名が付いていますが擬似乱数とは似て非なるものです。超一様分布列なんて言い方をする人も居ます。 以下のサイトも参考になるかと。 http://www.math.tohoku.ac.jp/akama/2006/lds.html

参考URL：
http://www.math.tohoku.ac.jp/akama/2006/lds.html
noname#40366
質問者

お礼

回答ありがとうございます。 やはり収束性の点からは凖モンテカルロ法がいいのですね。 教えていただきましたページをもとにサンプルプログラムを 作成しました。うまく動いているようです。

  1. カテゴリ
  2. 学問・教育
  3. 数学・算数

関連するQ&A

  • モンテカルロ法

    モンテカルロ法を使って計算を行う際に、 ある方法を用いて乱数列を発生させ、それを用いて行いました。 この乱数列は50の区分に分けた後、 カイ２乗分布による危険率１％の頻度検定 に合格したものを用いたのですが、 これではモンテカルロ法による計算を行うには 正しいとはいえないとのことでした。 これはなぜなのでしょうか。 確かに、50区分の頻度検定だけでは、 極端な話、周期50の数列がちょうど50区分に分かれて発生していて 結果として一様乱数となっているような場合も考えられます。 それでは、周期が計算に必要な乱数の数に対して十分に大きいとき、 モンテカルロ法として正しいと言えるようになるのでしょうか？ 解説をよろしくお願いします。

  • 乱数　メルセンヌツイスターについて

    メルセンヌツイスター法（MT法）について少し調べることになったのですが、メルセンヌツイスターのサイトにC言語のソースがあったので動かしてみました。 確かに乱数が発生するのですが、ソースをみてもどういう動きで乱数を生成してるのかがわかりません。 メルセンヌツイスターはどうやって乱数をつくりだしているのでしょうか？ わかる方がいらっしゃいましたら、教えてください。

  • 乱数生成、メルセンヌツイスターの使い方

    http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html 乱数生成にメルセンヌツイスターというものを使おうとしたのですが、 色々試してみてもいまいち使い方が分かりませんでした。 例えば０から１００までの間の乱数を得たい場合などはどのように 書けばよいのでしょうか・・？

  • メルセンヌツイスターを使った2次元乱数

    Mersenne Twisterを使って2次元の乱数を下記のように 生成しています。 1. 乱数を取得. x座標の値とする。 2. 1)で用いた乱数生成を利用して乱数を取得. y座標の値とする。 こうした作成したx,y座標のデータを見ますと、一様性が あまりないように見えます。 これは、2次元の乱数の扱いが間違っているのでしょうか? あるいは、周期が非常に長い乱数でも、2次元的に一様性を 保つためには、凖乱数を使うのがいいのでしょうか。

  • メルセンヌツイスターによる乱数の使い方

    メルセンヌツイスターによる乱数を生成したいのですが、下記コードにすると、mt_rand関数を呼ぶたびにシードが初期化され他乱数が生成されます。 やりたいことは、main関数で一度シード１０１で初期化した後は、mt_rand関数内では、初期化することなく継続した乱数を生成したいのです。 そのためa_mt_rand関数のようにすると、mtが定義されていないとなるエラーとなります。 namespaceの問題と思うのですが、色々試してみましたが、できませんでした。 どのようにすればできるでしょうか。 vc++11、windows11　使用 参考サイト C++ 乱数ライブラリ std::random の使い方 リンクはうまく貼れませんでした。 #pragma hdrstop #include <iostream> #include <stdio.h> #include <random> using namespace std; void mt_rand(void); int main(int argc, char *argv[]) { int ptr; std::mt19937 mt(101); // メルセンヌツイスターの32ビット版、引数は初期シード std::uniform_int_distribution <> rand100(0, 100); // [0, 99] 範囲の一様乱数 ptr = rand100(mt); printf("ptr=%d\n",ptr); mt_rand(); } void mt_rand(void) { int ptr=0; std::mt19937 mt(101); // メルセンヌツイスターの32ビット版、引数は初期シード std::uniform_int_distribution <> rand100(0, 100); // [0, 99] 範囲の一様乱数 ptr = rand100(mt); printf("%d\n",ptr); } void a_mt_rand(void) { int ptr=0; ptr = rand100(mt); printf("%d\n",ptr); }

  • 乱数発生ルーチンの使い方について

    数値計算において一様乱数を発生させるルーチンがいろいろあります。ソースが公開されているものやコンパイラが提供したりするものです。それらを利用する場合、乱数発生のシーズ（種）を与えてそれに応じて動作するというものが多いだろうと思います。そこで質問ですが、10000個の乱数を1回発生させる場合と100個の乱数を100回発生させる場合とで乱数の感じがかなり違います。いずれの場合も100×100の2次元データ（エクセルのシート状）として出力して作図したらその違いが簡単に分かります。この違いの原因はシーズの与え方が1回と100回という違いだろうと思います。100回のシーズの与え方にパターンが出来てしまうからだと思われます。例えば時間を使ってシーズを与えなおすことも考えられますが、今時のPCだとあっという間なのでシーズが同じだから、同じ乱数が100個できてしまいます。乱数を繰り返し発生させるときにその繰り返しの中でパターン化された乱数にならないように発生させる方法がないでしょうか。シーズが要らない乱数生成ルーチンとかですが。あるいはシーズをランダムに取得する方法が含まれたルーチン（シーズがないように見える）などです。あるいは本当にないものなど。メルセンヌツイスターはどうなのでしょうか。一応、フォートランでの利用を考えていますが、言語依存の問題ではないかもと思いますが。 よろしくお願いします。

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

    今、指数関数の積分をモンテカルロ法を用いて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)

  • n角形の重心を求めるアルゴリズム

    平面2次元のn角形の頂点のデータがあります。n点の座標ですから(x,y)がn個並んでいます。そのような図形の図心（重心）の座標を計算するアルゴリズムがないでしょうか。最終的にはプログラムとして離散的な処理をするため、１％ぐらいの誤差は許容範囲です。n角形と言ってもせいぜいn=3,4,5,6程度です。 欲を言うと、３次元も考えており、平面に含まれることが分かっているn個の点（３次元空間内）を平面の２次元空間に変換して重心を求め、それを３次元空間に引き戻せば３次元での重心となります。そのためにも２次元での重心の座標を求めるアルゴリズムが必要なのです。 よろしくお願いします。

  • 乱数の初期化について

    Cでモンテカルロシミュレーションを行っています。 乱数Merssenne twisterをつかって、各試行ごとに乱数を 時間で以下のように初期化しています。 init_genrand ((unsigned)time(NULL)) ; だいたい10000～100000回くらいシミュレーションを行う必要があって、 各試行は使用しているPCだと一瞬で終了します。 このときに上記の方法で種を初期化すると、 きちんと確認はしていないのですが、 試行時間が短すぎで種の時間が進んでいないような状態が起こります。 パソコンの時間で初期化する場合、間隔が短すぎると初期化種がかぶる ことはありますか。 もしあるなら、1回の試行時間が短いシミュレーションを各々初期化できる ような初期化の方法を教えていただけないでしょうか。 環境はcore2duo2.16、vista32、コンパイラはVisualC++2008です。 よろしくお願いします。

  • ベクトルの積分『∫(A→B)Fｄｒ』の計算の仕方を教えてください。

    皆様、こんにちは。 ベクトル積分『∫(A→B)Fｄｒ』の計算の仕方を教えて頂きたいのですが。 A、Bは空間の座標でそれぞれ（Xa、Ya、Za）、(Xb、Yb、Zc)。 F、ｒも３次元のベクトルで、 ｒは空間の位置のベクトルです。 Fの例として Fx(Fのｘ座標)＝Cx/(x^2+y^2+z^2)^(3/2) Fy(Fのy座標)＝Cy/(x^2+y^2+z^2)^(3/2) Fz(Fのz座標)＝Cz/(x^2+y^2+z^2)^(3/2) （Cはただの定数です。） で教えて頂けるとありがたいです。 （この例だと答えが分かっていますので。） 多分この例の答えは ーC/(Xa^2+Ya^2+Za^2)^(1/2)+C/(Xb^2+Yb^2+Zb^2)^(1/2)です。 少し数値が複雑そうですが、物理の保存力の基本問題です。 ベクトルの積分の計算方法をご存知の方やり方を教えて頂けないでしょうか。よろしくお願いします。（２次元の時は何となく分かった気になっていたのですが、３次元になって全く分からなくなってしまいました。）

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する