【ベイズ統計】実装を見て学ぶMCMC法 メトロポリスヘイスティング法（MH法）とハミルトニアンモンテカルロ法（HMC法）

　ベイズ統計で躓き易いのは、恐らく以下の二点かと思われます。

• 分布からのサンプリングシミュレーション法（MCMC法)
• サンプリングシミュレーション法を用いた事後分布の推定

　ここでは理論の詳細は他書に任せ、実装を見て、あぁこんなもんか、と思える事を目指します。
　理論については、豊田秀樹著　基礎からのベイズ統計学　ハミルトニアンモンテカルロ法による実践的入門がおすすめです。
　しかし、理論の説明は端的でとても良いのですが、アルゴリズムの説明が非常に不親切で、行間を読み、それを補完する知識のない初学者には結構難しいと感じたのでこの記事を書きます。
　
　本記事では、

• メトロポリスヘイスティング法（一番シンプル？用途は限られる）
• HMC法による分布からのサンプリング
• HMC法による多次元事後分布の母数ベクトル推定

のみを扱います。

　プログラムはC++で記述されており、依存する外部ライブラリは以下の通りです。

• Eigen
  • 【環境】Windows11、Visual Studio 2022、vcpkg
  • Eigenは他環境でも簡単に導入できると思います

　ではプログラムを書きます。
　後述するソースコード内では、事前分布や事後分布並びにサンプリング対象分布が固定されたものですので、別の分布を用いたい場合は別途実装する必要があります。

• メトロポリスヘイスティング法

#include <iostream>
#include <fstream>
#include <cmath>
#include <random>

#define M_PI 3.14159265358979323846

// ターゲット分布（正規分布）
// ここでは平均0、分散1の正規分布
double target_dist(double _x) {
	return std::exp(-0.5 * _x * _x) / std::sqrt(2 * M_PI);
}

// メトロポリスヘイスティング法
void metropolis_hastings(int _num_samples, double _sigma, double _theta) {
	std::random_device rd;
	std::mt19937 gen(rd());
	std::normal_distribution<> proposal_dist(0.0, _sigma); // 提案分布
	std::uniform_real_distribution<> uniform_dist(0.0, 1.0);

	std::ofstream samples("samples.csv");
	double theta_t = _theta;
	for (int i = 0; i < _num_samples; ++i) {
		double a = theta_t + proposal_dist(gen);// 提案された新しいサンプル
		double r = target_dist(a) / target_dist(theta_t);// 受容確率の計算
		// 新しいサンプルを受容するかどうか
		if (uniform_dist(gen) < r)//r未満の領域の乱数が出る確率はrだからサンプルを更新
			theta_t = a;//サンプルを更新
		// サンプルを出力
		// 出力されたサンプルはExcelやグラフライブラリで「ヒストグラム表示」しないと所望の形状を得られない
		samples << theta_t << std::endl;
	}
}

int main() {
    int num_samples = 1000;   // サンプル数
    double sigma = 1.0;       // 提案分布の標準偏差
    double initial_x = 0.0;   // 初期値
    metropolis_hastings(num_samples, sigma, initial_x);
    return 0;
}

• ハミルトニアン・モンテカルロ法（分布からのサンプリング）

#include <iostream>
#include <cmath>
#include <random>
#include <vector>
#include <Eigen/Dense>

#define M_PI 3.14159265358979323846

// ランダムエンジンの設定
std::mt19937 rng(std::random_device{}());
std::normal_distribution<> normal_dist(0.0, 1.0);
std::uniform_real_distribution<> uniform_dist(0.0, 1.0);

// 目的のポテンシャルエネルギー関数
double potential(const Eigen::VectorXd& x) {
	// 多次元正規分布を例としたエネルギー関数
	return 0.5 * x.squaredNorm();
}

// ポテンシャルエネルギーの勾配（xについての偏微分）
Eigen::VectorXd potential_gradient(const Eigen::VectorXd& x) {
	return x;  // 多次元正規分布の場合、勾配は単純に x になります
}

// ハミルトニアンモンテカルロ法のステップ
Eigen::VectorXd hmc_step(
	const Eigen::VectorXd& current_x
	, double step_size
	, int num_steps
) {
	int dim = current_x.size();
	// 運動量の初期化（標準正規分布からサンプリング）
	Eigen::VectorXd p(dim);
	for (int i = 0; i < dim; ++i) p(i) = normal_dist(rng);
	// Leapfrog法を用いた更新
	Eigen::VectorXd x = current_x;
	Eigen::VectorXd p_half = p - 0.5 * step_size * potential_gradient(x);
	for (int i = 0; i < num_steps; ++i) {
		x += step_size * p_half;
		if (i != num_steps - 1) {
			p_half -= step_size * potential_gradient(x);
		}
	}
	p_half -= 0.5 * step_size * potential_gradient(x);
	// 提案された位置と運動量でハミルトニアンを計算
	double current_H = potential(current_x) + 0.5 * p.squaredNorm();
	double proposed_H = potential(x) + 0.5 * p_half.squaredNorm();
	// メトロポリス・ヘイスティングス法による遷移の採択
	if (uniform_dist(rng) < std::exp(current_H - proposed_H)) {
		return x;  // 提案位置を採択
	}
	else {
		return current_x;  // 現在の位置を保持
	}
}

int test_hmc(void)
{
	int dim = 2;  // 次元数
	int num_samples = 1000;
	double step_size = 0.1;
	int num_steps = 10;
	// 初期位置
	Eigen::VectorXd x = Eigen::VectorXd::Zero(dim);
	// サンプリングの実行
	for (int i = 0; i < num_samples; ++i) {
		x = hmc_step(x, step_size, num_steps);
		std::cout << "Sample " << i + 1 << ": " << x.transpose() << std::endl;
	}
	return 0;
}

• ハミルトニアン・モンテカルロ法による多次元事後分布の母数推定

// ランダムエンジンの設定
#include <iostream>
#include <cmath>
#include <random>
#include <vector>
#include <Eigen/Dense>

#define M_PI 3.14159265358979323846

std::mt19937 rng(std::random_device{}());
std::normal_distribution<> normal_dist(0.0, 1.0);
std::uniform_real_distribution<> uniform_dist(0.0, 1.0);

// 観測データ（例として生成されたデータ）
std::vector<double> observed_data = { 
	5.1, 5.5, 5.3, 5.8, 5.2 
};

// ログ尤度関数
double log_likelihood(const Eigen::VectorXd& params) {
	double mu = params(0);
	double sigma_sq = params(1);
	double log_likelihood = 0.0;
	for (double x : observed_data)
		log_likelihood += -0.5 * std::log(2 * M_PI * sigma_sq) - (x - mu) * (x - mu) / (2 * sigma_sq);
	return log_likelihood;
}

// ログ事前分布
double log_prior(const Eigen::VectorXd& params) {
	double mu = params(0);
	double sigma_sq = params(1);
	// μ ~ N(0, 10^2)
	double log_prior_mu = -0.5 * (mu * mu) / 100.0 - 0.5 * std::log(2 * M_PI * 100.0);
	// σ^2 ~ Inverse-Gamma(α=2, β=2)
	if (sigma_sq <= 0) return -INFINITY; // σ^2 > 0 の条件をチェック
	double alpha = 2.0;// 2.0;
	double beta = 2.0;// 2.0;
	double log_prior_sigma_sq = -(alpha + 1) * std::log(sigma_sq) - beta / sigma_sq;

	return log_prior_mu + log_prior_sigma_sq;
}

// ハミルトニアン
double hamiltonian(const Eigen::VectorXd& params, const Eigen::VectorXd& momentum) {
	return -log_likelihood(params) - log_prior(params) + 0.5 * momentum.squaredNorm();
}

// ポテンシャルエネルギーの勾配＝h'()
/// 平均と分散だけが母数なので多次元ベクトルを扱えるようには実装されていない
Eigen::VectorXd grad(const Eigen::VectorXd& params) {
	double mu = params(0);
	double sigma_sq = params(1);
	Eigen::VectorXd grad(2);
	// μの勾配
	double grad_mu = 0.0;
	for (double x : observed_data)
		grad_mu += (x - mu) / sigma_sq;
	// 事前分布による勾配
	grad_mu -= mu / 100.0;
	// σ^2の勾配
	double grad_sigma_sq = 0.0;
	for (double x : observed_data)
		grad_sigma_sq += -0.5 / sigma_sq + (x - mu) * (x - mu) / (2 * sigma_sq * sigma_sq);
	// 事前分布による勾配
	grad_sigma_sq -= (2.0 + 1) / sigma_sq + 2.0 / (sigma_sq * sigma_sq);
	//
	grad(0) = grad_mu;
	grad(1) = grad_sigma_sq;
	return grad;
}

// HMCにおける一回の母数ベクトル更新処理
/// 基礎からのベイズ統計学　豊田秀樹　p.111参照
/// 実装自体は簡単なので注釈の多さにビビらないように
Eigen::VectorXd hmc_step(
	// リープ・フロッグ法で探索する母数ベクトル（推定したい母数）
	/// θ
	const Eigen::VectorXd& _theta
	// リープフロッグ法による母数更新時に使う小さな値
	, double _ep
	//　リープフロッグ法のステップ数
	, int _L
) {
	int dim = _theta.size();
	// 運動量の初期化（標準正規分布からサンプリング）
	Eigen::VectorXd ps(dim);
	for (int i = 0; i < dim; ++i)
		// ps = p, 母数ごとに異なる運動量pを用いて更新するために標準正規分布から値を生成している
		ps(i) = normal_dist(rng);
	// リープフロッグ法を用いた更新
	/// 与えられた初期母数ベクトルから始める
	Eigen::VectorXd th = _theta;
	/// １．p(τ + 0.5) = p(τ) - 0.5 * ε * h'(θ(τ)) 式(5.31)
	Eigen::VectorXd p_half = ps - 0.5 * _ep * grad(th);
	for (int i = 0; i < _L; ++i) {
		/// ２．θ(τ + 1) = θ(τ) + ε * p(τ + 0.5) 式(5.32)
		th += _ep * p_half;
		// (_L - 1)番目のステップ以外の運動量の更新
		// 最終ステップの運動量更新はループを出た後に行う -> (＠)
		if (i != _L - 1)
			/// ３．p(τ + 1) = p(τ + 0.5) ‐ 0.5 * ε * h'(θ(τ + 1)) 式(5.33)
			/// ここでは、式(5.33)と式(5.31)を同時に行う為、0.5をかけずに更新している
			/// ループするので、式(5.33)のθ(τ+1)とp(τ+1)は、次のステップにおける式(5.31)
			/// のθ(τ)とp(τ)になることから、計算を一つにまとめることが出来る
			/// なので、0.5をかけていなくても誤りではない
			p_half -= _ep * grad(th);
	}
	// （＠）ループ内で生じた運動量と位置のズレを補正　
	// p.111の図5.4を参照
	// 
	// この実装の更新ループはp(τ + 0.5)から始まり、最終ステップ以外でθ、p共
	// に1単位時間ごとに更新されるので、最終ステップも含めて同様に更新すると、
	// ループ終了時点で、運動量だけが0.5単位時間分余計に進んでしまう
	// そこで、最終ステップだけ運動量の更新を行わないでおくと、運動量が0.5
	// 単位時間分だけ遅れ、その分をループを抜けた時点で計算する
	// 
	// 少し具体的に示すと、
	// 
	// １．ループに入る手前で、運動量の0.5単位時間更新　式(5.31)
	// ２．ループ内で、位置（母数ベクトル）の１単位時間更新　式(5.32)
	// ３．ループ内の最終ステップ以外で、運動量の１単位時間更新
	// 　　・式(5.33）と式(5.31)を同時に計算してしまっている
	// 　　・１の更新から始まる更新なので、運動量が位置より0.5単位時間常に早い
	// 　　・最終ステップではこの更新を行わないので、ループを抜けた時点で運動量
	// 　　　が位置に対して0.5単位時間遅れる
	// 
	// という順序で行われているため、ループを抜けた段階では最終ステップ
	// の運動量更新が行われないので、運動量だけが0.5単位時間遅れている
	// その運動量の0.5単位時間遅れのみをここで計算している
	p_half -= 0.5 * _ep * grad(th);
	// ハミルトニアンを計算して遷移を決定
	/// todo : 受理確率の計算処理を入れたい
	double current_H = hamiltonian(_theta, ps);
	double proposed_H = hamiltonian(th, p_half);
	if (uniform_dist(rng) < std::exp(current_H - proposed_H))
		return th;// 提案された母数ベクトルを返す。推定母数ベクトルを更新するということ
	else
		return _theta;// 引数の（現在の）母数ベクトルをそのまま帰す。推定母数ベクトルを更新しないということ
}

int test_hmc(void)
{
	int num_samples = 1000;
	double step_size = 0.005;
	int num_steps = 30;
	// 初期パラメータ (μ, σ^2)
	/// 調整しないとうまく受理されず探索できない
	/// 当初の、0.0、1.0、という値では無理だった
	/// 事前分布のパラメータは変更しなくても探索はしたが、場合によっては必要らしい
	Eigen::VectorXd params(2);
	params << 5.38, 1.0;
	// サンプリングの実行
	for (int i = 0; i < num_samples; ++i) {
		//paramsはこのループ内で更新され続ける(HMC法で受容されなかった場合は更新されない）
		params = hmc_step(params, step_size, num_steps);
		std::cout << "Sample " << i + 1 << ": μ = " << params(0) << ", σ^2 = " << params(1) << std::endl;
	}
	return 0;
}

　どうでしょうか？

　プログラム自体はとても簡単ですね。

　注釈が長い部分もありますが、これは、当記事の目的が、MCMC法の実装の学習であるため、シンプルなプログラムに至るまでにどの様な思考を経ているのか説明するためです。
　これらはベイズ統計に必要なMCMC法とそれを用いた母数推定法に過ぎず、実際にベイズ統計を用いるには、サンプリング結果から事後統計量を計算し（単純な統計量計算に帰着します。文献[1]を参照）、なんらかの結論を導き出します。
　そのあたりの解説は、冒頭で紹介した豊田先生の書籍（参考文献[1]）を参照してください。
　また、ここに掲載したプログラムが正しい保証はありませんので、その点だけご注意ください。

• 参考文献
  • [1]　豊田秀樹著　基礎からのベイズ統計学　ハミルトニアンモンテカルロ法による実践的入門

本記事著者のブログはこちらです。気が向いたらいらしてください。