Boost で正規乱数の生成

シミュレーション実験を行っていると正規乱数が必要になることがよくあります．正規乱数はかなり基本的な乱数であるためさまざまなライブラリに実装されているのですが，私は Boost をよく使っているので，正規乱数でも Boost を使用して発生させています．
しかし Boost の実装は柔軟性を求め過ぎるあまり使いやすさを失ってしまっているように感じます（単に正規乱数を発生させたい人にとっては）．というのも，Boost で正規乱数を発生させるためにはまず一様分布を発生させるエンジンを選び，それを適切にコンストラクトする必要があります．エンジンによってはテンプレート化されたコンストラクタをもつため，エンジンを適切に初期化するためにはコンストラクタに渡す種(seed)の型を正確に一致させる必要があります（そうしなければテンプレート化されたコンストラクタが呼ばれてしまいます）．その後 normal_distribution と variable_generator を構築するという手順を踏み，ようやく正規乱数を得ることができるのです．
このような骨が折れる作業も，簡単なラッパークラスを作れば楽になります．以下にそのソースコードを示します．

gauss_rand.hpp

#ifndef GAUSS_RAND_HPP_20081010
#define GAUSS_RAND_HPP_20081010

#if defined(_MSC_VER) && (_MSC_VER >= 1020)
# pragma once
#endif


#include <boost/random.hpp>


namespace math {
    template <class RealType = double, class Engine = boost::mt19937>
    class gauss_rand {
    private:
        template <class T>
        struct seed_traits {
            typedef typename T::result_type seed_type;
        };

        template <class RT, int w, unsigned int p, unsigned int q>
        struct seed_traits<boost::random::lagged_fibonacci_01<RT, w, p, q> > {
            typedef boost::uint32_t seed_type;
        };

    public:
        typedef boost::variate_generator<Engine, boost::normal_distribution<RealType> > generator_type;

        typedef typename generator_type::engine_value_type engine_value_type;
        typedef typename generator_type::engine_type       engine_type;
        typedef typename generator_type::distribution_type distribution_type;
        typedef typename generator_type::result_type       result_type;

        typedef typename seed_traits<Engine>::seed_type    seed_type;

        explicit gauss_rand(seed_type seed, const result_type& mean = result_type(0), const result_type& sigma = result_type(1))
            : generator_(Engine(seed), boost::normal_distribution<RealType>(mean, sigma))
        {}

        void seed(seed_type value) { generator_.engine().seed(value); }

        result_type mean()  const { return generator_.distribution().mean();  }
        result_type sigma() const { return generator_.distribution().sigma(); }

        result_type operator()() { return generator_(); }

        engine_value_type&       engine()             { return generator_.engine();       }
        const engine_value_type& engine()       const { return generator_.engine();       }

        distribution_type&       distribution()       { return generator_.distribution(); }
        const distribution_type& distribution() const { return generator_.distribution(); }

        generator_type&          generator()          { return generator_;                }
        const generator_type&    generator()    const { return generator_;                }

    private:
        boost::variate_generator<Engine, boost::normal_distribution<RealType> > generator_;
    };


    template <class RealType>
    class gauss_rand<RealType, boost::ecuyer1988> {
    public:
        typedef boost::variate_generator<boost::ecuyer1988, boost::normal_distribution<RealType> > generator_type;

        typedef typename generator_type::engine_value_type engine_value_type;
        typedef typename generator_type::engine_type       engine_type;
        typedef typename generator_type::distribution_type distribution_type;
        typedef typename generator_type::result_type       result_type;

        typedef boost::ecuyer1988::result_type             seed_type;

        explicit gauss_rand(seed_type seed1, seed_type seed2, const result_type& mean = result_type(0), const result_type& sigma = result_type(1))
            : generator_(boost::ecuyer1988(seed1, seed2), boost::normal_distribution<RealType>(mean, sigma))
        {}

        void seed(seed_type value1, seed_type value2) { generator_.engine().seed(value1, value2); }

        result_type mean()  const { return generator_.distribution().mean();  }
        result_type sigma() const { return generator_.distribution().sigma(); }

        result_type operator()() { return generator_(); }

        engine_value_type&       engine()             { return generator_.engine();       }
        const engine_value_type& engine()       const { return generator_.engine();       }

        distribution_type&       distribution()       { return generator_.distribution(); }
        const distribution_type& distribution() const { return generator_.distribution(); }

        generator_type&          generator()          { return generator_;                }
        const generator_type&    generator()    const { return generator_;                }

    private:
        boost::variate_generator<boost::ecuyer1988, boost::normal_distribution<RealType> > generator_;
    };
}


#endif

gauss_rand.hpp に正規乱数を計算するファンクタ gauss_rand を実装しました．gauss_rand.hpp を include するだけで gauss_rand が使用可能になります．
gauss_rand の第一テンプレート引数に float ã‚„ double を渡すことにより戻り値の型を，第二テンプレート引数に minstd_rand ã‚„ mt19937 を渡すことにより一様分布のエンジンを指定することが出来ます．Boost の normal_distribution ではここで指定した一様分布を用いてボックス＝ミューラー法(Box-Muller transform)により正規乱数を計算することになります．

以下のテストコードで gauss_rand の動作を確認できます．

test.cpp

#include "gauss_rand.hpp"
#include <iostream>


int main()
{
    // 最も単純な使い方
    // 種=0x5eed, 平均=0, 標準偏差=10 の正規乱数
    // （戻り値は double, 一様乱数は mt19937 が使われます）
    {
        math::gauss_rand<> grand(0x5eed, 0, 10);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=5, 標準偏差=7 の正規乱数
    // 戻り値=float
    // （一様乱数は mt19937 が使われます）
    {
        math::gauss_rand<float> grand(0x5eed, 5, 7);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=10, 標準偏差=20 の正規乱数
    // 戻り値=double
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に ministd を用いる
    {
        math::gauss_rand<double, boost::minstd_rand> grand(0x5eed, 10, 20);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=1, 標準偏差=2 の正規乱数
    // 戻り値=long double
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に rand48 を用いる
    {
        math::gauss_rand<long double, boost::rand48> grand(0x5eed, 1, 2);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種1=0x5eed1, 種2=0x5eed2, 平均=11, 標準偏差=0 の正規乱数
    // 戻り値=double
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に ecuyer1988 を用いる
    {
        math::gauss_rand<double, boost::ecuyer1988> grand(0x5eed1, 0x5eed2, 11, 0);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=100, 標準偏差=5 の正規乱数
    // 戻り値=float
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に kreutzer1986 を用いる
    {
        math::gauss_rand<float, boost::kreutzer1986> grand(0x5eed, 100, 5);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=5, 標準偏差=100 の正規乱数
    // 戻り値=float
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に hellekalek1995 を用いる
    {
        math::gauss_rand<float, boost::hellekalek1995> grand(0x5eed, 5, 100);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=10, 標準偏差=12 の正規乱数
    // 戻り値=float
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に mt19937 を用いる
    {
        math::gauss_rand<float, boost::mt19937> grand(0x5eed, 10, 12);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    // 種=0x5eed, 平均=10, 標準偏差=5 の正規乱数
    // 戻り値=double
    // 正規乱数の発生アルゴリズムで用いる一様乱数の生成に lagged_fibonacci19937 を用いる
    {
        math::gauss_rand<double, boost::lagged_fibonacci19937> grand(0x5eed, 10, 5);

        for (int i = 0; i < 100; ++i) {
            std::cout << grand() << std::endl;
        }
    }

    return 0;
}

• VC9
• g++ (GCC) 3.4.4 (cygming special)
• g++ (GCC) 4.1.3

でのコンパイルを確認しています．
実行してみると，正規乱数が求められていることが分かります．
gauss_rand.hpp はここからダウンロードできます．

次回は OpenCV の画素アクセスについて紹介しようと思います．
それでは．