正規分布と正規乱数 - chmod 777 myknowledge

ノイズを付加する場合に確率密度関数として正規分布を用いることがある。
これを実際にやってみたくなったので実現方法を調べてみた。自分への覚え書きとして以下に示す。

正規分布の説明は以下を参照。
http://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83

正規分布を表す式は以下。

正規分布に従うノイズを発生させるためには、上式における x をノイズ量として x 軸上の各点におけるノイズを f(x) で表される確率で発生させればよい。
一方、多くのプログラミング言語で提供されているのは一様乱数(ある有限の区間を区切って、その区間内で全ての実数が同じ確率（濃度）で現れるような乱数 by Wikipedia)を発生させるメソッドである。C言語でいえば stdlib.h で提供される rand() が相当する。
0より大きく(0は含まない)1以下の、すなわち(0,1]の一様乱数を正規乱数(正規分布を持つ乱数)に変換する方法として、ボックス＝ミューラー法(Box-Muller transform)がある。

ボックス＝ミューラー法では、一様乱数(0,1]の要素αとβを次式で変換する。
(参照：http://ja.wikipedia.org/wiki/%E4%B9%B1%E6%95%B0#.E6.AD.A3.E8.A6.8F.E4.B9.B1.E6.95.B0)

こうして二つの直交するN(0,1)の正規乱数が得られる。
この正規乱数に標準偏差σを乗じて、平均μを加えることで、平均μ、分散 $\sigma^2$ の正規分布N(μ, $\sigma^2$ )に従う正規乱数が得られる。

ボックス＝ミューラー法で行っていることを以下に直観的に説明する（数式による厳密な証明についてはここでは言及しない）。以下を参考にしている。
http://en.literateprograms.org/Box-Muller_transform_(C)

単位円内の点は(0,1]の一様乱数αとβを用いてα・cos(2πβ)とα・sin(2πβ)で表現できる。一様乱数を用いているので単位円内の各点における存在確率に偏りは無い。ここでボックス＝ミュラー法の式を見ると、αではなく $sqrt{-2 \cdot log\alpha}$ となっている。これは直観的には下図の様に、原点を中心とする単位円内の各円を別の原点を中心とする円に射影することになる。射影先の円では原点に行くほど蜜で外側に行くほど疎になっている。これにより単位円内において存在確率に偏りが無かったものが、射影先では外側に行くほど存在確率が下がる円になる。射影先の円における存在確率はN(0,1)の正規分布に従う。

ボックス＝ミュラー法で正規乱数を生成するプログラムをC++で書いたので以下に掲載する。

box_muller_transform.h

#ifndef ___BOX_MULLER_TRANSFORM_H___
#define ___BOX_MULLER_TRANSFORM_H___

#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <time.h>

using namespace std;

/**
 * @brief 正規乱数(Normally Distributed Random Numbers)を生成するためのクラス
 *        Box-Muller transform アルゴリズムに基づく。
 */
class NormDistRand{
 private:
  double m_n1; // 2組の一様乱数 を Box-Muller transform した値1
  double m_n2; // 2組の一様乱数 を Box-Muller transform した値2

  bool   m_flg_n2_cached;

  double m_mean;   // 正規分布の平均
  double m_stddev; // 正規分布の標準偏差

  static const double PI;
 public:

  /// constructor
  NormDistRand(double mean, double stddev) : 
  m_mean(mean), m_stddev(stddev), m_flg_n2_cached(false)
    { srand(static_cast<unsigned int>(time(NULL))); }

  /// constructor
  NormDistRand(double mean, double stddev, unsigned int seed) : 
  m_mean(mean), m_stddev(stddev), m_flg_n2_cached(false)
    { srand(seed); }

  /**
   * @brief 正規乱数値を生成して返す。
   * @param [out] 正規乱数
   */
  double rand(void);

 private:

  /// inhibit using default constructor.
  NormDistRand();

  /**
   * @brief Box-Muller transform を実施して一様乱数を m_n1, m_n2 に変換する。
   */
  void _box_muller_transform(void);
};

#endif

box_muller_transform.cc

#include "box_muller_transform.h"

const double NormDistRand::PI = 3.141592653589793116;

// ------------------------------------------------------------------------------- //
double NormDistRand::rand(void) {
  if(!m_flg_n2_cached) {
    _box_muller_transform();
    m_flg_n2_cached = true;
    return m_n1 * m_stddev + m_mean;
  }
  else{
    m_flg_n2_cached = false;
    return m_n2 * m_stddev + m_mean;
  }
}

// ------------------------------------------------------------------------------- //
void NormDistRand::_box_muller_transform(void) {
  double alpha = 0.0;
  while(alpha == 0.0){
    alpha = static_cast<double>(std::rand()) / RAND_MAX;
  }
  double beta  = static_cast<double>(std::rand()) / RAND_MAX;
  m_n1 = sqrt(-2 * log(alpha)) * sin(2 * PI * beta);
  m_n2 = sqrt(-2 * log(alpha)) * cos(2 * PI * beta);
  return;
}

以上のプログラムを用いて、正規乱数を発生させてヒストグラムを作成するサンプルプログラムを以下に示す。

box_muller.cc

#include "box_muller_transform.h"
#include <string>
#include <vector>

void usage(const string& prgm) {
  cout << "\n"
       << "Usage : " << prgm << "  mean  stddev  hist_size  hist_step  gene_time\n"
       << "=============================================\n"
       << "  @brief 正規乱数を生成してヒストグラムを出力する。\n"
       << "=============================================\n"
       << "   mean      : 正規分布の平均\n"
       << "   stddev    : 正規分布の標準偏差\n"
       << "   hist_size : ヒストグラムは数直線を単位区分長hist_stepで区切って単位区分ごとの度数で表現する。\n"
       << "               単位区分の個数をhist_sizeで指定。\n"
       << "   hist_step : 上記単位区分長\n"
       << "   gene_time : サンプル数を指定\n"
       << endl;
}

int main(int argc, char* argv[]) {
  if(6 != argc) {
    usage(argv[0]);
    return EXIT_FAILURE;
  }

  int n = 0;
  double mean      = atof(argv[++n]);
  double stddev    = atof(argv[++n]);
  int    hist_size = atoi(argv[++n]);
  double hist_step = atof(argv[++n]);
  int    gene_time = atoi(argv[++n]);

  vector<int> hist(2 * hist_size + 1, 0);

  NormDistRand ndr(mean, stddev);
  double val;
  int index;
  for(int i = 0; i < gene_time; i++) {
    val = ndr.rand();
    index = static_cast<int>( (val + hist_size * hist_step) / hist_step + 0.5 );
    if(index < hist.size() && index >= 0){
      hist[index]++;
    }
  }

  double pos = -1 * hist_size * hist_step;
  for(int i = 0; i < hist.size(); i++) {
    cout << pos << "\t" << hist[i] << endl;
    pos += hist_step;
  }
  return EXIT_SUCCESS;
}

上記サンプルプログラムの実行ファイル box_muller を用いてヒストグラムファイルを以下のように作成する。

./box_muller 0.0 0.447213595 500 0.01 1000000 > mean0.0_var0.2.dat
./box_muller 0.0 1.0 500 0.01 1000000 > mean0.0_var1.0.dat
./box_muller 0.0 2.236067978 500 0.01 1000000 > mean0.0_var5.0.dat
./box_muller -2.0 0.707106781 500 0.01 1000000 > mean-2.0_var0.5.dat

ここで第２引数には分散ではなく、標準偏差を入力していることに注意。分散としては上から、0.2, 1.0, 5.0, 0.5 を指定しているのと同義。この値は Wikipedia での正規分布の例の図中に掲載されている分散にあわせてある(2008年5月11日現在)。

このように作成したデータを gnuplot でプロットする。

plot "mean0.0_var0.2.dat", "mean0.0_var1.0.dat", "mean0.0_var5.0.dat", "mean-2.0_var0.5.dat"

作成した図を以下に示す。Wikipedia
http://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83
の図と比較すると同じような分布を描いているのを確認できる。