• ベストアンサー

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

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

noname#40366
noname#40366

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

  • ベストアンサー
noname#108554
noname#108554
回答No.3

>​http://phi.med.gunma-u.ac.jp/swtips/Rsim/​ >作成したグラフは上記ページの二つめの図(の2次元版) >のようになります。 2つ目の図は「ほぼ均一にプロットされているように見える。」と書いてあるとおり、 うまくいった例ですから、このぐらい均一なら問題ないのでは? これ以上、均一な分布が欲しいのなら、格子にしてしまうか、 超一様分布(集中講義で聞いた覚えがありますが、用途は主に数値積分か。ただし、本当に一様分布より均一なのかは不明。生成アルゴリズムも不明。) ぐらいしかないと思います。

noname#40366
質問者

お礼

何度も御手数おかけします。 現在まさに数値積分(またたは収束の早い平均値の計算)を 目指しています。 1) Monte Carloでできるだけ精度を上げた計算 (MT使用) 2) Quasi-Monte Carloを用いた計算 の二通りを進めており、1)の計算方法に問題がないかが 心配する部分でした。やりとりの間で、極端に間違った ことをしていないということがわかり少し安心です。 2)の方法もプログラムができており、比較した一様性は No.2のお礼の所で書きました図と同様の結果がえられています。 色々ありがとうございました。

その他の回答 (2)

noname#108554
noname#108554
回答No.2

こちらの環境はC++(OSはlinux)です。 再度やってみました。 やはりほとんど一様に見えます。 具体的にどのように偏っているのでしょうか? 例えば、原点付近にあまり点がないとか、ムラが大きいとか。 できれば、一様性の検定して「ほーらやっぱり一様でない」と 確認できれば一番いいのですが、それは面倒ですかね・・・ >同じ乱数生成器grnd()をxとyに同時に使用している点 それは問題ないはずです。

noname#40366
質問者

お礼

ご確認ありがとうございます。 一様性の計算の仕方は勉強不足のためすぐにできそう にはありません。 また、画像が貼れないため、かわりにwebページで似たような 結果を探しました。 http://phi.med.gunma-u.ac.jp/swtips/Rsim/ 作成したグラフは上記ページの二つめの図(の2次元版) のようになります。 確かに変な周期性は見えないのですが、必ずしも均質な 分布でないように思えます。黒が濃く見える部分と薄く 見える部分ができています。 Monte Carlo+MTとQuasi-Monte Carloの2次元プロットにて 計算した結果では、下記ページにある左と右の図のように 明確な違いが見えます。 http://public.lanl.gov/kmh/uncertainty/meetings/warnpix.gif Monte Carlo+MTを間違って使って、悪い分布になってしまって いるということでないかを知りたかったのですが、ここまで の所では問題はなさそうなのですね。

noname#108554
noname#108554
回答No.1

問題ないはずです。 「メルセンヌツイスタによって生成された2つの乱数を 2次元プロットしてみた」ということでよいと思いますが、 それは以前やってみたことがあります。 よければ、プログラムを見せてください。 だめなら、「1)で用いた乱数生成を利用して」の部分を もう少し詳しく説明してください。

noname#40366
質問者

お礼

回答ありがとうございます。 プログラム(Fortran)は単純に下記のようにしています。ここで grnd()がMTから乱数を読み取る部分です。 do i=1,10000 x = grnd() y = grnd() print *, x,y end do 同じ乱数生成器grnd()をxとyに同時に使用している点 など問題ないでしょうか?

関連するQ&A

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

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

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

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

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

    メルセンヌツイスターによる乱数を生成したいのですが、下記コードにすると、mt_rand関数を呼ぶたびにシードが初期化され他乱数が生成されます。 やりたいことは、main関数で一度シード101で初期化した後は、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); }

  • 多次元正規乱数が正しく生成されているか確認するには?

    多次元正規乱数に従うサンプルを生成したく,以下の資料を参考にし,プログラムを作成しました。生成したサンプルが本当に要求した多次元正規乱数にしたがっているか確認したいのですが,以下の手順で正しいといえるでしょうか? 例: 平均μ=[1; 2] 分散協分散行列 Σ=[1,0;    0,1]; に従う多次元正規乱数y=[y1,y2]を全部でn=100サンプル生成した. このとき,n=100個のy1の平均が1,n=100個のy2の平均が2となり, n=100個のy1の分散が1,n=100個のy2の分散が1,y1とy2の共分散が0となれば正しく生成できているといってよいでしょうか? この方法で確認したのですが,生成したサンプルが元の平均,共分散から0.1~0.2のずれがあり正しいのか不安です。 参考にした資料: 「多次元正規乱数の作成法」 http://ci.nii.ac.jp/naid/110001186428/ja 「共分散の計算の仕方」 http://www.ge.fukui-nct.ac.jp/~nagamizu/f-3-s.pdf

  • excelの乱数を用いて円が重ならない座標を選び出

    excelの乱数を用いて円が重ならない座標を選び出す方法 2次元平面(0<x<100,0<y<100)に同じ半径rの円を重ならないように配置したいのですが、困ってます。 考えとしては 1.まず1個の円の中心座標を決める(R1(x1,y1)とする) 2.2個目の円を、1個目と重ならないよう配置できる中心座標R2(x2,y2)を乱数で決める。 ((x1-x2)^2+(y1-y2)^2>4r^2) 3.3個目の円を、1個目と2個目と重ならないように配置できる中心座標R3(x3,y3)を決める。 ((x1-x3)^2+(y1-y3)^2>4r^2かつ(x2-x3)^2+(y2-y3)^2>4r^2) これを繰り返す。 このような座標をexcelの乱数機能を用いて出すことは可能でしょうか? 皆様のお力をお貸しください。

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

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

  • 乱数の平均

    十分大きいxに対して、0~1のx個の乱数をつくります。 そして、この値をラインプロットします。 次にこの値の(隣り合うy個の)移動平均を書きます。すると、 yが小さいとき乱数を少しなめした構造になると考えられます。 yが十分大きい値の場合0.5になると思われます。 しかしその中間状態などはよくわかりません。 数学的にはいったいどうなっているのでしょうか?? 数学にうといのでわかりやすく教えてください。 お願いします。

  • 乱数

    乱数x、y(0〈x〈1)、(0〈y〈1)をn回発生させて、これらを座標とする点をx‐y座標平面上にプロットしてこの点と原点との距離が1以下となるときの回数をカウントしm回としてn:m=(長さ1の正方形の面積):(円の4分の1の面積)となる。これを利用して円周率を求めるプログラムなんてどう作るんですか?

  • 3次元座標XYZを視覚化したい。

    3次元座標(X,Y,Z)の値がいくつかあるのですが、 それを3Dでグラフ化するようなソフトはないでしょうか?

  • こんな条件を満たす乱数生成関数教えてください

    1.任意の周期を指定できる 2.種を指定できる(直前の生成値を引数にとる) 3.逆関数が定義できる 4.生成された乱数 x、y の距離を(定数時間で)求められる   つまり y = f(x) ならxとyの距離は1、 y = f( f( f(x) ) ) なら距離3、というように 乱数としての質(均等に分布していること)はあまり重視しません。 ビット幅は32~128bitくらい(任意ならベスト)であればいいと思っています。 以下のような感じにしたいです。  int rand(x, p);   // 戻り値 y = f(x)、pは周期、xは直前の乱数値  int inv(y, p);   // 戻り値 x = f^-1(y)  int distance(x, y) // y = f(f(x)) のとき、distance(x, y) = 2 で distance(y, x) = -2 一応以下の関数が条件1~4を満たすのですが、残念ながら乱数としての性質が皆無なので使えないです。  int rand(x, p) { return (x+1) % p; }  int inv(y, p) { return y ? y - 1 : p-1; }  int distance(x, y) { return y - x; } よろしくご教授お願いします。