※学習メモのため、中身の検証等は行っていない点にご留意ください

今回は統計モデルの係数の決定方法の一つとしてのベイズ推定に関して述べる。
ベイズ推定とはそもそも何かというと、事前分布を仮定し、それに基づいた分析により尤度分布を導出、最後にベイズの定理に従い、事前分布と尤度分布により事後分布を導出するというものである。根本的な考え方としては、最終的に求めたい確率分布を推定するにあたって、元々の仮説を新しい情報を得るに従って更新するというという考え方である。これにより、通常例えば最尤推定や最小二乗法等による偏回帰係数の推定を行うようなモデルであっても、偏回帰係数等の推定が出来る。

MCMC法

現在では、ベイズ推定の方法としてはMCMC法（マルコフ連鎖モンテカルロ法）がメジャーであると認識している。これは特にアルゴリズムとしてメトロポリタン・ヘイスティング法を用いて、ランダムにパラメーターを指定した上で、少しでも事後確率が高くなるようなパラメーターを採用し続ける事により、当てはまりのいいパラメーターを推定するといった方法である。これにより、非常に歪で綺麗な確率分布を仮定しにくいデータに対しても、分析を行う事が出来る。

例えば以下のような形である

mcmc.cpp

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

// 事前確率分布を定義するための関数
double prior(double theta) {
    // この例では一様分布を使います
    return 1.0;  // 正規化は無視します
}

// データの尤度を定義するための関数
double likelihood(double theta, const std::vector<double>& x, const std::vector<double>& y) {
    double result = 1.0;
    double sigma = 1.0;  // データのノイズを仮定します
    for (int i = 0; i < x.size(); ++i) {
        double expected_y = theta * x[i];  // 線形モデルを仮定します
        double diff = y[i] - expected_y;
        result *= (1.0 / std::sqrt(2.0 * M_PI * sigma * sigma)) * std::exp(-(diff * diff) / (2.0 * sigma * sigma));  // 正規分布の確率密度関数
    }
    return result;
}

int main() {
    std::random_device rd;  // ランダムなシードを生成するための非決定的な乱数生成器
    std::mt19937 gen(rd());  // メルセンヌ・ツイスター法に基づく決定的な乱数生成器
    std::uniform_real_distribution<> dis(-10.0, 10.0);  // -10から10までの一様分布

    // データを準備します
    std::vector<double> x = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    std::vector<double> y = {2.3, 4.2, 6.1, 8.4, 10.1, 12.4, 14.2, 16.1, 18.3, 20.2};

    // パラメータの初期値を設定します
    double theta = dis(gen);  // ランダムな初期値

    // メトロポリス・ヘイスティングス法を実行します
    int num_iterations = 1000;
    for (int i = 0; i < num_iterations; ++i) {
        double new_theta = dis(gen);  // パラメータの新しい候補
        double prob = (likelihood(new_theta, x, y) * prior(new_theta)) / (likelihood(theta, x, y) * prior(theta));  // 受容確率
        std::uniform_real_distribution<> accept_dis(0.0, 1.0);
        if (accept_dis(gen) < prob) {
            theta = new_theta;  // パラメータを更新します
        }
    }

    std::cout << "Estimated parameter: " << theta << std::endl;

    return 0;
}

ここでポイントなのは、 if (accept_dis(gen) < prob) 　の部分。
直感的には、受容確率が1を超えた時（事後確率が事前確立よりも大きい）場合に、そのパラメーターを採用するといったものが理解しやすいと思われる。つまり、このif文は　if(1 < prob) 　であるべきなのでは？といった事を疑問として持ちやすいかと思われる。しかし、この場合はaccept_dis(gen)は0から1までの乱数であり、場合によっては事前確率の方が大きい場合でも採用してしまう事になる。これはどういうことか。

実はこれはマルコフ連鎖の性質と関連して、敢えてこのような方法を取っている。どういう事かというと、一定程度の確率で、事前確率＞事後確率となるようなパラメーターを採用する事で、全パラメーターを隈なく探索するといった事が可能になるからである。深い理由に関しては調査中であるが、これによりエルゴード性と詳細つり合い条件が保たれて、それにより確率分布が定常分布に収束する事が保証される事を目指すことが目的。ベイズ推定のゴールはパラメーターの確率分布を一つに絞る事だから、定常分布を仮定することは必須条件なのである。

※エルゴード性は、マルコフ連鎖が全ての状態を十分な時間が経てば確率1で訪問でき、そしてその状態から他の状態に移動できることを意味する。
詳細つり合い条件は、任意の2つのパラメータ値iとjについて、iからjへの遷移確率とjからiへの遷移確率が一定の比（具体的には、それぞれのパラメータ値における事後確率の比）を保つ事を要求する条件。

マルコフ連鎖と定常分布に関して理解するためには、確率行列を理解するとわかりやすい。そちらに関しても下記に示す。

確率行列と定常確率ベクトル

確率行列はマルコフ連鎖の遷移確率を表す正方行列である。これに対し主張されていることとしては、この確率行列の任意の要素が正であれば、定常確率ベクトル（左固有値が1となる固有ベクトル）がただ一つ存在するといったものである。そして、この定常確率ベクトルは確率行列を無限乗する事によって生成されるとされている。これはペロン・フロベニウスの定理によって証明されている。

参考ページ
https://www.headboost.jp/markov-chain-montecarlo/
https://math-fun.net/20210728/16780/