超一様分布列をBoost.Randomを用いて実装する

ご存知の通り、C++のライブラリには高性能な擬似乱数生成エンジンBoost.Randomが存在している。Boost.Randomは擬似乱数を生成するエンジン部分と、一様分布、ガウス分布、ベルーヌーイ分布など好みの分布を生成する出力部分の2つで構成されている。

しかし、ことモンテカルロ法においては主に収束速度を早める目的として、擬似乱数ではなく低食い違い量列(Low-Discrepancy Sequence, LDS)と呼ばれる数列が用いられる場合がある。これを用いると低次元での積分計算の際に収束速度を大きく改善することができ、具体的には従来のモンテカルロ法で(1/√N)だった収束速度が約(1/N)くらいにまで落とすことができる。

さて、このような低食い違い量列をBoost.Randomで用いるためにはエンジン部分を製作すればいいだけで、実際には以下のような形で実装すれば良い。今回は簡単なHalton列を用いて実装を行った:

	#pragma once
	#include <boost/cstdint.hpp>
	#include <boost/config.hpp>
	#include <cmath>
	namespace aoba {
	struct halton_sequence_engine {
	public:
	typedef double result_type;
	halton_sequence_engine(): m_index(1), m_base(2) {}
	halton_sequence_engine(uint32_t index): m_index(index), m_base(2) {}
	halton_sequence_engine(uint32_t index, uint32_t base)
	: m_index(index), m_base(base) {}
	result_type operator()() {
	result_type result = 0;
	result_type f = 1 / static_cast<result_type>(m_base);
	uint32_t i = m_index;
	while(i > 0) {
	const uint32_t div_i_base = i / m_base;
	const uint32_t mod_i_base = i % m_base;
	result = result + f*mod_i_base;
	i = div_i_base;
	f = f / m_base;
	}
	++m_index;
	return result;
	}
	void seed(uint32_t index) { m_index = index; }
	static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return 0; }
	static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return 1; }
	private:
	uint32_t m_index;
	uint32_t m_base;
	};
	} // namespace aoba

view raw halton_sequence_engine.hpp hosted with ❤ by GitHub

具体的には、対象のファンクタがUniform Random Number Generatorのコンセプトを満たすためには以下の記述が必要となる:

• X::result_typeにoperator()が返す型を指定
• operator()にエンジン部分を記述
• min(), max()にエンジンが生成する分布の最小値、最大値を指定

さらに、対象のファンクタが上述の内容に加えPseudo-Random Number Generatorのコンセプトを満たすためには:

• X()にデフォルトの状態で生成されるようなコンストラクタを指定
• X(...)にユーザ指定の状態で生成されるコンストラクタを指定
• seed(...)にエンジンの状態を変更するためのメンバ関数を指定

するだけで良い。後はBoost側が良きに計らってくれる。非常に便利だと思う。

余談
C++0xには同等のstd::randomが付属しているため、これを利用すればいいのではという声もあるとは思うが、std::uniform_real_distribution<>にこのLDSを食わせたところ、どうみてもLDSではない値が出力される結果となった。恐らくresult_typeがdoubleという一般の擬似乱数ではない型を指定しているために起こる現象だと思うが、ソースを見ていないので一先ずは従来のBoost.Randomの手法を用いることにする。