読者です 読者をやめる 読者になる 読者になる

真実の楽譜(フルスコア)

プログラム関係の忘備録になるはず

C#でボックスミューラー法による正規分布に従う乱数生成

.NET Frameworkでは一様乱数を生成するクラスとして
System.Randomクラスが提供されています。
Random クラス (System)

データにノイズを加えるときに一様乱数よりも、
正規分布に従うノイズ(ホワイトガウス)を加えたいというケースがありますが、
.NET Frameworkではそこまでは提供してくれていないので、
ボックスミューラー法による正規分布に従う乱数生成クラスを作ってみました。

ボックスミューラー法

ボックスミューラー法のアルゴリズムは単純で、
(0~1]の範囲を取る一様乱数を2つ用意し(R1,R2)、
以下の式に当てはめるだけで、
N(0.0,1.0)の正規分布に従う乱数が1組生成できます。

{ \displaystyle
Norm_1 = \sqrt{-2 \log R_1} \cos{2 \pi R_2}
}
{ \displaystyle
Norm_2 = \sqrt{-2 \log R_1} \sin{2 \pi R_2}
}

また平均値μ、分散σ^2に従う正規分布N(μ,σ^2)での乱数は、
得られた値にσを乗算して、平均値を加えることで求まります。

{ \displaystyle
Norm_1 = \sigma \sqrt{-2 \log R_1} \cos{2 \pi R_2} + \mu
}
{ \displaystyle
Norm_2 =  \sigma \sqrt{-2 \log R_1} \sin{2 \pi R_2} + \mu
}

参考
http://ja.wikipedia.org/wiki/乱数列
http://ja.wikipedia.org/wiki/ボックス=ミュラー法

プログラムソース

RandomBoxMuller.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace RandomBoxMuller
{
  public class RandomBoxMuller
  {
    private Random random;

    public RandomBoxMuller()
    {
      random = new Random(Environment.TickCount);
    }

    public RandomBoxMuller(int seed)
    {
      random = new Random(seed);
    }

    public double next(double mu = 0.0, double sigma = 1.0, bool getCos = true)
    {
      if (getCos)
      {
        double rand = 0.0;
        while ((rand = random.NextDouble()) == 0.0) ;
        double rand2 = random.NextDouble();
        double normrand = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Cos(2.0 * Math.PI * rand2);
        normrand = normrand * sigma + mu;
        return normrand;
      }
      else
      {
        double rand;
        while ((rand = random.NextDouble()) == 0.0) ;
        double rand2 = random.NextDouble();
        double normrand = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Sin(2.0 * Math.PI * rand2);
        normrand = normrand * sigma + mu;
        return normrand;
      }
    }

    public double[] nextPair(double mu = 0.0, double sigma = 1.0)
    {
      double[] normrand = new double[2];
      double rand = 0.0;
      while ((rand = random.NextDouble()) == 0.0) ;
      double rand2 = random.NextDouble();
      normrand[0] = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Cos(2.0 * Math.PI * rand2);
      normrand[0] = normrand[0] * sigma + mu;
      normrand[1] = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Sin(2.0 * Math.PI * rand2);
      normrand[1] = normrand[1] * sigma + mu;
      return normrand;
    }
  }
}

実行例

作成したクラスはRandom関数と同様に、
コンストラクタにシード値を指定させてオブジェクトを初期化し、
(未指定の場合はシステムの起動時刻をシード値として利用)
next関数の引数に正規分布のμとσを指定することで乱数が求まります。

Program.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace RandomBoxMuller
{
  class Program
  {
    static void Main(string[] args)
    {
      var random = new RandomBoxMuller();

      using (var sr = new System.IO.StreamWriter("output.csv"))
      {
        for (int i = 0; i < 10000; i++)
        {
          sr.WriteLine(random.next());
        }
      }
    }
  }
}


上記プログラムではN(0.0,1.0)の乱数を1万個生成しています。
得られたデータのヒストグラムは以下のようになりました。

f:id:s_sikisya:20140603132934p:plain

今回-1.0~1.0のデータは6816個得られており、
Z=1.0のときの正規分布データ点数は6826個(0.3413)
なのでほぼ近い値が得られていることが確認できました。

参考
http://www.koka.ac.jp/morigiwa/sjs/standard_normal_distribution.htm