La mémoire de mathématiques

数学めも by Müde

Box-Muller法による正規分布列生成

大抵のプログラミング言語には、擬似乱数を生成するための関数が用意されています。そして、よくあるパターンとしては、その擬似乱数は0.0以上1.0以下の適当な値を返します。ヒストグラムを描くと分かりますが、大抵の場合はどの区間も同じ程度抽出されます(一様分布)。

さて、世の中には正規分布という分布が有ります。正規分布についての説明は端折りますが、正規分布に従う適当な値を出してくれるような便利な関数が常に用意されているとは限りません。一様分布から標準正規分布に変換する方法として、Box-Muller法というのが知られています。それを実際に実装してみましょう。

WikipediaのBox-Muller法を見れば分かりますが、方法としては(0,1)の一様分布から2つの確率変数(すなわち、乱数)を用意して、示されている式に従ってZ1もしくはZ2を求めれば良さそうです。

例えばZ1であれば、{Z_1 = \sqrt{-2 \log X} \cos 2 \pi Y}と書かれているので、これを素直にJavaのプログラムで表現して、z1 = Math.sqrt(-2.0 * Math.log(x)) * Math.cos(2.0 * Math.pi * y); とすれば良いですね。ということで、実際のJavaのコードは次の通りです。


BoxMuller.java … Box-Muller法によって正規分布列を生成するクラス

import java.util.Random;

public class BoxMuller {
 private Random random;
 private double mu;
 private double sigma;
 private boolean useCosine;
 
 /**
  * Constructor
  * @param mu
  * @param sigma
  * @param useCosine
  */
 public BoxMuller(double mu, double sigma, boolean useCosine){
  random = new Random(System.currentTimeMillis());
 
  this.mu = mu;
  this.sigma = sigma;
  this.useCosine = useCosine;
 }

 /**
  * get next value
  * @return
  */
 public double next(){
  double x = random.nextDouble();
  double y = random.nextDouble();
  double result = 0.0;
 
  if(useCosine){
   result = Math.sqrt(-2.0 * Math.log(x)) * Math.cos(2.0 * Math.PI * y);
  } else {
   result = Math.sqrt(-2.0 * Math.log(x)) * Math.sin(2.0 * Math.PI * y);
  }
  return result * sigma + mu;
 }
}

BoxMullerTest.java … Box-Mullerクラステスト用プログラム

public class BoxMullerTest {
 public static void main(String[] args) {
  BoxMuller boxMuller = new BoxMuller(0,1,true);
  for(int i=0;i<10000;++i){
   System.out.println(boxMuller.next());
  }
 }
}

10000回繰り返して値を取り出し、プロットしたグラフが次のとおりです。なんだか良さそうに見えますね。

f:id:mu-de:20140911104345p:plain