折れ線で定義した確率密度関数に従う乱数列を生成するプログラム

概要

モンテカルロ法を使う際、任意の確率密度関数に従う乱数列を生成したいことがあります。 C++ には、折れ線で定義した確率密度関数から乱数列を生成できる標準クラス piecewise_linear_distribution があるのですが、下記３点が気になりましたので自作してみました。

• Visual Studio 2012 は累積分布関数の逆関数を求めるところにバグがあるらしく、(折れ線を等間隔で定義していても)正しく計算できない。 参考サイト
• Visual Studio 2017、2019 にもバグが残っているらしく、折れ線を等間隔で定義していない場合は正しく計算できない・・・(標準関数でもこういうことがあるみたいなので、正しく動くかどうかは確認が必要ですね) 参考サイト
• もう少し処理速度を速めたい

処理速度の問題について

確率密度関数に従う乱数列を生成する方法の１つに「逆関数法」があり、piecewise_linear_distribution もその手法に基づいて実装されているようです。 確率密度関数を積分する⇒逆関数を求める⇒逆関数に一様乱数を入力したときに得られる値を乱数列として使用する、というシンプルな手順で実現できるため強力なのですが、 piecewise_linear_distribution では確率密度関数を折れ線で定義する都合上、逆関数を求める際に探索が必要となり、それが速度的にネックになるというデメリットも持っています。 標準クラスでは２分探索で実現しているようですので、２分探索よりも速く探索できる方法を考えることになります。 本プログラムでは、逆関数の横軸上でほぼ等間隔になるような逆引きインデックス列を生成して、インデックスで当たりをつけてから線形探索する、という方法を採ってみました。

検証結果

下記サンプルコードで1千万回乱数を生成したときの処理時間の比較です。 (単位はms:ミリ秒) サンプル数 33～ 65537 の範囲では自作版の方が速いことが確認できます。大体 3 ～ 4 倍くらいでしょうか？

＜検証環境＞

• コンパイラ： Visual Studio 2017、最適化オプションは /O2
• CPU：Core i7-6700
• メモリ：16GB

自作 piecewise_linear_distribution ＆ サンプルソースコード

標準クラスと同じ名前ですが、仕様は異なりますのでご注意ください。一様乱数の生成にはSIMDxorshift(Apache License)を使わせていただきました。 (速度を検証していたところ、C++標準のメルセンヌツイスターの処理が結構重いことが判明したので、もっと速いものを探したところ見つけた次第です。 これはこれでモンテカルロ法を高速に実行しようと思った際に欠かせないアイテムに思います)。

ライセンスは Boost v 1.0 とします。

なお、続編として折れ線関数の簡略化を書きました。piecewise_linear_distribution に機能追加もしてますので、ソースを参照される場合はそちらもご参照して頂くと良いかもしれません。

/*
 * inverse sample
 * Copylight (C) 2023 mocchi
 * License: Boost ver.1
 */

#include <cstdio>
#include <cstdint>
#include <vector>
#include <cmath>
#include <random>
#include <chrono>
#include "simdxorshift128plus.h"

#define USE_STL 0
#define SAMPLES_FROM_PFUNC 2048

struct xorshift_rnd_32bit {
	avx_xorshift128plus_key_t key;
	int cnt;
	__m256i rnd;
	xorshift_rnd_32bit(uint64_t key1, uint64_t key2) {
		avx_xorshift128plus_init(key1, key2, &key);
		cnt = 8;
	}
	double operator()() {
		if (cnt == 8) cnt = 0, rnd = avx_xorshift128plus(&key);
		return static_cast<double>(rnd.m256i_u32[cnt++]) / (4294967295.0);
	}
};

struct xorshift_rnd_32bit_stl {
	avx_xorshift128plus_key_t key;
	int cnt;
	__m256i rnd;
	xorshift_rnd_32bit_stl(uint64_t key1, uint64_t key2) {
		avx_xorshift128plus_init(key1, key2, &key);
		cnt = 8;
	}
	uint32_t operator()() {
		if (cnt == 8) cnt = 0, rnd = avx_xorshift128plus(&key);
		return rnd.m256i_u32[cnt++];
	}
	inline uint32_t min() const {
		return 0;
	}
	inline uint32_t max() const {
		return 4294967295;
	}
};


struct piecewise_linear_distribution {
	std::vector<double> v_ary, p_ary;

	std::vector<double> p_accum_ary;
	std::vector<int> p_accum_revlookup;
	double p_sum;
	piecewise_linear_distribution() : p_sum(0) {}
	piecewise_linear_distribution(double *v_ary_, double *p_ary_, int data_count, int lookup_count) {
		init(v_ary_, p_ary_, data_count, lookup_count);
	}
	// v_ary_, p_ary_ を新たに指定して計算用データを作成する
	void init(double *v_ary_, double *p_ary_, int data_count, int lookup_count) {
		v_ary.resize(data_count), p_ary.resize(data_count);
		for (int i = 0; i < data_count; ++i) {
			v_ary[i] = v_ary_[i], p_ary[i] = p_ary_[i];
		}
		init(lookup_count);
	}
	// 入力済の v_ary, p_ary を基に計算用データを作成する
	void init(int lookup_count) {
		int data_count = static_cast<int>(p_ary.size());
		p_accum_ary.resize(data_count - 1);
		// 累積配列の作成
		double *v_ary_ = v_ary.data();
		double *p_ary_ = p_ary.data();
		double *p_accum_ary_ = p_accum_ary.data();
		p_sum = 0;
		for (int i = 0; i < data_count - 1; ++i) {
			double v_delta = v_ary_[i + 1] - v_ary_[i];
			double p_accum = p_sum + (p_ary_[i + 1] + p_ary_[i]) * 0.5 * v_delta;
			p_accum_ary_[i] = p_accum;
			p_sum = p_accum;
		}
		double *p_accum_begin = p_accum_ary_;
		double *p_accum_end = p_accum_ary_ + data_count - 1;

		// 逆引きインデックス配列の作成
		p_accum_revlookup.resize(lookup_count);
		int *p_accum_revlookup_ = p_accum_revlookup.data();
		double *p_accum_iter = p_accum_begin;
		for (int i = 0; i < lookup_count; ++i) {
			double pa = static_cast<double>(i) * p_sum / static_cast<double>(lookup_count - 1);
			p_accum_iter = std::lower_bound(p_accum_iter, p_accum_end, pa);
			p_accum_revlookup_[i] = static_cast<int>(p_accum_iter - p_accum_begin);
		}
	}
	template <typename F> double operator() (F &rnd) {
		double a_dist = rnd();
		double rv_idx_i = std::floor(a_dist * static_cast<double>(p_accum_revlookup.size() - 1));
		double t = rnd();
		return sample_detail(a_dist, rv_idx_i, t);
	}
	double sample(double p) {
		double a_dist = p;
		double rv_idx_i;
		double t = std::modf(a_dist * static_cast<double>(p_accum_revlookup.size() - 1), &rv_idx_i);
		return sample_detail(a_dist, rv_idx_i, t);
	}
private:
	inline double sample_detail(double a_dist, double rv_idx_i, double t) {
		double a_p = a_dist * p_sum;
		int rv_idx = static_cast<int>(rv_idx_i);
		int rv = p_accum_revlookup[rv_idx];
		while (p_accum_ary[rv] < a_p && rv < p_accum_ary.size() - 1) ++rv;
		double *p_ary_ = p_ary.data();
		double p0 = p_ary_[rv], p1 = p_ary_[rv + 1];
		if (p0 != p1) {
			t = (std::sqrt(p0 * p0 * (1.0 - t) + p1 * p1 * t) - p0) / (p1 - p0);
		}
		double *v_ary_ = v_ary.data();
		double v0 = v_ary_[rv], v1 = v_ary_[rv + 1];
		return v0 + t * (v1 - v0);
	}
};

int main() {
	auto &func_p = [](double x) {return (x < 0 || x>1) ? 0 : 1.0 - std::cos(x*2.0*3.14159265358979); };
	std::vector<double> v_ary, p_ary;
	std::printf("probability density\n");
	std::printf("index\tx\ty\n");

	v_ary.clear(), p_ary.clear();
	// 確率密度配列を作成
	for (double x = 0; x <= 1.0; x += 1.0/static_cast<double>(SAMPLES_FROM_PFUNC)) {
			v_ary.push_back(x);
			p_ary.push_back(func_p(x));
		}
	for (size_t i = 0; i < v_ary.size(); ++i) {
		std::printf("%zu\t%f\t%f\n", i, v_ary[i], p_ary[i]);
	}

#if USE_STL == 0
	piecewise_linear_distribution dist(v_ary.data(), p_ary.data(), static_cast<int>(v_ary.size()), static_cast<int>(v_ary.size()));
	xorshift_rnd_32bit rnd(344, 2583);
#else
	std::piecewise_linear_distribution<> dist(v_ary.data(), v_ary.data() + v_ary.size(), p_ary.data());
	xorshift_rnd_32bit_stl rnd(344, 2583);
#endif
	// 検証
	int test_count = 10000000;
	std::vector<double> histogram(500);
	auto time_start = std::chrono::system_clock::now();
	for (size_t i = 0; i < test_count; ++i) {
		double v_idx = dist(rnd);
		int h_idx = static_cast<int>(std::floor(v_idx * static_cast<double>(histogram.size() - 1) + 0.5));
		histogram[h_idx]++;
	}
	auto time_end = std::chrono::system_clock::now();

	int64_t duration = std::chrono::duration_cast<std::chrono::milliseconds>(time_end - time_start).count();
	std::printf("elapsed_time:%lld [msec]\n", duration);
	std::printf("elapsed time per sample:%g [msec/sample]\n", static_cast<double>(duration) / static_cast<double>(test_count));

	std::printf("montecarlo test\n");
	for (size_t i = 0; i < histogram.size(); ++i) {
		std::printf("%f\t%f\n", static_cast<double>(i) / static_cast<double>(histogram.size() - 1), histogram[i] / static_cast<double>(test_count));
	}
}