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

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

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

ボックスミューラー法

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

${ \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());
        }
      }
    }
  }
}