ボックス＝ミュラー法による正規分布の乱数生成

namespace Tips_Random

{

    using System;

    using System.Linq;

    class Program

    {

        static void Main(string[] args)

        {

            #region ボックス＝ミュラー法による正規分布の乱数生成

            Console.WriteLine("ボックス＝ミュラー法による正規分布の乱数生成器を作る。");

            Console.WriteLine("指定された平均値と標準偏差を持つ分布の乱数を生成できる。");

            Console.WriteLine("output.csv に出力しているので Excel などでその分布が確認できる。");

            var rand2 = new BoxMullerRandom();

            var mu = 2.0;

            var sigma = 0.27;       // だいたい ±1.0 の分布になる

            using (var writer = new System.IO.StreamWriter("output.csv"))

            {

                var NN = 1000;

                var boxValues = Enumerable.Range(0, NN).Select(_ =>

                {

                    var value = rand2.Next(mu, sigma);

                    writer.WriteLine(value);

                    return value;

                }).ToArray();

                var boxAve = boxValues.Average();

                var boxStdEv = Math.Sqrt(boxValues.Select(x => x * x).Sum() / NN - boxAve * boxAve);

                Console.WriteLine("指定した値 : μ = {0}, σ = {1}", mu, sigma);

                Console.WriteLine("計算した値 : boxAve = {0}, σ = {1}", boxAve, boxStdEv);

                Console.WriteLine("実際の平均値と標準偏差が指定した値に近くなっている。");

            }

            Console.WriteLine("");

            Console.ReadKey();

            #endregion ボックス＝ミュラー法による正規分布の乱数生成

        }

    }

    /// <summary>

    /// ボックス＝ミュラー法による正規分布の乱数生成器を表します。

    /// </summary>

    public class BoxMullerRandom

    {

        /// <summary>

        /// 乱数生成器

        /// </summary>

        private Random _rand;

        /// <summary>

        /// 新しいインスタンスを生成します。

        /// </summary>

        public BoxMullerRandom()

            : this(Environment.TickCount)

        {

        }

        /// <summary>

        /// 新しいインスタンスを生成します。

        /// </summary>

        /// <param name="seed">乱数生成用のシードを指定します。</param>

        public BoxMullerRandom(int seed)

        {

            this._rand = new Random(seed);

        }

        /// <summary>

        /// 指定された分布にしたがった乱数を返します。

        /// </summary>

        /// <param name="mu">平均値を指定します。</param>

        /// <param name="sigma">標準偏差を指定します。</param>

        /// <param name="isGetCos">余弦波を用いる場合に true を指定します。</param>

        /// <returns>指定された条件にしたがった乱数を返します。</returns>

        public double Next(double mu = 0.0, double sigma = 1.0, bool isGetCos = true)

        {

            var rand = 0.0;

            while ((rand = this._rand.NextDouble()) == 0.0) ;

            var rand2 = this._rand.NextDouble();

            var normRand = GetNormRandValue(rand, rand2, isGetCos);

            return normRand * sigma + mu;

        }

        private double GetNormRandValue(double rand, double rand2, bool isGetCos = true)

        {

            return Math.Sqrt(-2.0 * Math.Log(rand)) * (isGetCos ? Math.Cos(2.0 * Math.PI * rand2) : Math.Sin(2.0 * Math.PI * rand2));

        }

        /// <summary>

        /// 指定された分布に従った乱数の組を返します。

        /// </summary>

        /// <param name="mu">平均値を指定します。</param>

        /// <param name="sigma">標準偏差を指定します。</param>

        /// <returns>指定された条件にしたがった乱数の組を返します。</returns>

        public double[] NextPair(double mu = 0.0, double sigma = 1.0)

        {

            var rand = 0.0;

            while ((rand = this._rand.NextDouble()) == 0.0) ;

            var rand2 = this._rand.NextDouble();

            return new double[]

            {

                GetNormRandValue(rand, rand2, true) * sigma + mu,

                GetNormRandValue(rand, rand2, false) * sigma + mu,

            };

        }

    }

}