乱数

fortran90においての乱数の発生方法です。

乱数を1から作るのは、乱数を研究するとかでない限り止めましょう。
ここではfortran90で乱数を発生させるためにどうすればいいか?を記述します。

私は乱数について詳しくありませんが、調べた限りメルセンヌツイスタという乱数発生法が良い乱数を発生するようです。
メルセンヌツイスタの詳細はMersenne Twister Home Page にあるのでご参照ください。
fortran90コードはhttps://gist.github.com/ykonishi/5569005にあります。

上記のリンクから、メルセンヌツイスタのfortran90コード(mt19937.f90)をダウンロードします。
そして、メインプログラム

program main
  use mt19937
  implicit none
  integer::i
  double precision::d

  call sgrnd(2) ! seed
 
  do i=1,10
     write(6,*)grnd()
  enddo
 
  stop
end program main

をmt19937.f90と一緒に

gfortran -fno-range-check mt19937.f90 main2.f90

というオプションを付けてコンパイル(※1)します。実行すると、

$ ./a.out
   6.4690865343436599E-002
  0.40996560384519398    
  0.14283204451203346    
  0.74749888177029788    
  0.59957117028534412    
  0.53676494420506060    
  0.38914293772540987    
  0.94303490873426199    
  0.15882158395834267    
  0.15933134709484875    
$

という結果を得ることが出来ます。

(※1)
このオプションが、型が扱える範囲を超えていても無視するためのオプションです。
何故かgfortranコンパイラ(ver4.8.4で確認)で出ます。
そして、これはコンパイラの問題だと私は結論付けました。
fortranの整数は
-2147483648から2147483647までで定義されているのにもかかわらず、
gfortran(ver4.8.4)は
-2147483647から2147483647の範囲でなければエラーを出すようです。
しかし、実際に整数に-2147483648を代入し、このエラーを無視する上のオプションを付けて実行すると、ちゃんと値が入ります。なので、コンパイラに問題があると判断しました。


スポンサーリンク


任意の関数に従う乱数(von Neumannの棄却法)


このページのサムネイル画像の右は等方的2次元調和振動子励起状態の密度
\(
\begin{align}
f(x,y)=\frac{4}{\pi}x^2y^2e^{-(x^2+y^2)}
\end{align}
\)
を点の濃さで表現するためのプログラムです。

これは、上の関数に従う乱数分布を書くことにほかなりません。
一様乱数を任意の関数に従う乱数に変換する方法として、
von Neumannの棄却法と呼ばれる方法があります。

以下のようにすると一様乱数を関数\(f(x)\)に従う乱数に変えることが出来ます。

\(x\)の範囲\([x_a,x_b]\)で定義された関数\(f(x)\)に従う乱数を発生させたい。
区間\([x_a,x_b]\)で関数は\([f_a,f_b]\)の値を持つとすると、
 1. 範囲\([x_a,x_b]\)の一様乱数\(x_r\)を生成する。
 2. 範囲\([f_a,f_b]\)の一様乱数\(f_r\)を生成する。
 3. もしも\(f_r\le f(x_r)\)ならば、その\(x_r\)は\(f(x)\)に従う乱数である。そうでなければ、発生させた乱数は捨てる。
 4. 1-3を繰り返す。
という方法です。アルゴリズムを見て分かる通り、この方法は無駄が多いです。
発生させた乱数を捨てる操作は出来ればしたくありません。
ですが、この方法は非常に簡単に実装できるのでとりあえず実装したいという場合には良い方法です。

以下は、上の考えの二次元版です。

program main
  use mt19937
  implicit none
  integer::i,count,N
  double precision::xp,yp,fp
  double precision::a,b,fmin,fmax
  double precision,external::f
 
  call sgrnd(2)

  N=2000
 
  ! fmin :: minimum value of function f(x,...)
  ! fmax :: maxmum value of function f(x,...)
  fmin=0.d0
  fmax=1.d0

  ! range of random number [a,b]
  a=-5d0
  b= 5d0
 
  open(13,file="./density.d")
  count=0
  do while(count.lt.N)
     xp=grnd()*(b-a)+a
     yp=grnd()*(b-a)+a
     fp=grnd()*(fmax-fmin)+fmin
     
     if(fp.le.f(xp,yp))then
        write(13,*)xp,yp
        count=count+1
     endif
  enddo
  close(13)

  stop
end program main

function f(x,y)
  implicit none
  double precision::f
  double precision,intent(in)::x,y

  double precision::pi=acos(-1.d0)
 
  f=abs(sqrt((1.d0/pi))*2d0*x*y*exp(-0.5d0*(x*x+y*y)))**2
 
  return
end function f

カラーマップと3次元での表現はこうなります。
20150524-210920_1_c

乱数を用いた密度表示(下のプログラム)ではこう表現されます。
nxny11density_c

スポンサーリンク

参考

任意分布の乱数生成 https://www.slideshare.net/shogoosawa581/ss-65259938

以降はfortranの標準組み込み関数を用いる場合のものです。
消すのは忍びないため、一応おいておきます。


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です