マルコフ連鎖における特定時刻の状態到達確率を求めるC++プログラム

本記事では、マルコフ連鎖（Markov Chain）において、初期状態から出発し、指定された時間が経過した時点で特定の状態へ到達する確率を求めるC++プログラムについて解説します。

マルコフ連鎖とは

マルコフ連鎖とは、複数の「状態」と、ある状態から別の状態へ遷移する際の「遷移確率」から構成される確率過程（ランダムプロセス）です。状態間の遷移には単位時間が1単位かかるものとします。マルコフ連鎖の重要な特徴は、「次の状態」が現在の状態のみによって決まり、それ以前の履歴に依存しないというマルコフ性を持つ点です。

マルコフ連鎖は有向グラフとして表現できます。この問題を解くには、与えられたマルコフ連鎖を行列形式に変換します。この行列において、位置 (a, b) の要素は、状態「a」から状態「b」へ遷移する確率を表します。

遷移行列を用いると、時刻 t における確率分布は次の漸化式を使って再帰的に求められます。

P(t) = M × P(t-1)

したがって、遷移行列 M の T 乗を計算すれば、その各要素から T 単位時間後の状態間の到達確率を読み取ることができます。

C++による実装

以下のプログラムでは、行列の積を計算する multiply()、繰り返し二乗法によって行列のべき乗を高速に求める matrix_power()、そして初期状態から最終状態への到達確率を算出する calc_prob() の3つの関数を実装しています。

#include <bits/stdc++.h>
using namespace std;
#define float_vec vector<float>

// 2つの行列の積を計算する
vector<float_vec > multiply(vector<float_vec > A, vector<float_vec > B, int N) {
    vector<float_vec > C(N, float_vec(N, 0));
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            for (int k = 0; k < N; ++k)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

// 繰り返し二乗法により行列のべき乗を計算する
vector<float_vec > matrix_power(vector<float_vec > M, int p, int n) {
    vector<float_vec > A(n, float_vec(n, 0));
    for (int i = 0; i < n; ++i)
        A[i][i] = 1;
    while (p) {
        if (p % 2)
            A = multiply(A, M, n);
        M = multiply(M, M, n);
        p /= 2;
    }
    return A;
}

// 初期状態から最終状態への到達確率を計算する
float calc_prob(vector<float_vec > M, int N, int F, int S, int T) {
    vector<float_vec > matrix_t = matrix_power(M, T, N);
    return matrix_t[F - 1][S - 1];
}

int main() {
    vector<float_vec > G{
        { 0, 0.08, 0, 0, 0, 0 },
        { 0.33, 0, 0, 0, 0, 0.62 },
        { 0, 0.06, 0, 0, 0, 0 },
        { 0.77, 0, 0.63, 0, 0, 0 },
        { 0, 0, 0, 0.65, 0, 0.38 },
        { 0, 0.85, 0.37, 0.35, 1.0, 0 }
    };
    // 利用可能な状態の数
    int N = 6;
    int S = 4, F = 2, T = 100;
    cout << "Probability of reaching: " << F << " in time " << T << " after starting from: " << S << " is " << calc_prob(G, N, F, S, T);
    return 0;
}

出力

Probability of reaching: 2 in time 100 after starting from: 4 is 0.271464

この実行結果は、状態4から出発して時間100が経過した時点で、状態2へ到達している確率が約 0.2715（約27.15%）であることを示しています。

アルゴリズムのポイント

このプログラムの効率性と正確性を支えているのは、次の点です。

・遷移行列の要素 (a, b) には、状態 a から状態 b への遷移確率を格納します。
・行列のべき乗を素朴に計算すると O(T × N³) の計算量が必要ですが、繰り返し二乗法を用いることで O(N³ log T) まで削減できます。
・これにより、T = 100 のような大きな時間ステップでも効率的に確率を求められます。
・プログラム内では状態番号が1始まりで扱われるため、配列参照時にインデックスを1つずらして調整しています。