【備忘録】Kleinbergのバースト分析がだいぶ理解できた

SNSのバースト分析などでも頻繁に登場するKleinberg博士のバースト分析のロジック、論文見てもさっぱりでした。

いろんなサイト見たら、隠れマルコフモデルとビタビアルゴリズムがポイントということはわかったんですが、今度はビタビアルゴリズムが理解できない。。。

こんなダメダメ状態でしたが、この本とChat-gpt先生のお陰でやっとだいぶ理解できました。

このシリーズはトピックモデルも読んでいますが、数式の解説が丁寧で、文系出身には大変ありがたい本です。ありがたやありがたや。黄色本にはまだ手が出ないので、このシリーズで基礎力を底上げしていきたいと思います。こんな状態で４月から博士課程大丈夫か俺、いやだいじょばない。

アルゴリズムの理解のためChat-gpt先生にコードを出力してもらいましたが、V[{}] で早速つまずきました。でもこれがビタビテーブルで、ビタビアルゴリズムの肝なのですね。

簡潔に言うと、時間 × 状態の2次元構造の確率を記録するテーブルで、その確率値はすべて「その状態に最も高確率で到達するパスの確率」を記録したものです。

t(時刻)	Rainy（雨）	Sunny（晴れ）
0（初期）	0.06	0.24
1	0.048	0.0432
2	0.0168	0.003888

これを再帰的に更新していき、最後にそれぞれ時刻におけるの確率値の大きい方をかけていくことで、最後の状態に至る最尤確率とその状態に至るパスが分かると理解しました。

以下、Chat-gpt先生に出力してもらったコードです。最後の方でまだ理解できない部分がありますが、とりあえずはなんとなく理解できたので満足。

viterbi_algorism.py

import numpy as np

# 状態（hidden states）
states = ['Rainy', 'Sunny']

# 観測値（observations）
observations = ['walk', 'shop', 'clean']

# 初期確率（start probabilities）
start_probability = {'Rainy': 0.6, 'Sunny': 0.4}

# 状態遷移確率（transition probabilities）
transition_probability = {
    'Rainy': {'Rainy': 0.7, 'Sunny': 0.3},
    'Sunny': {'Rainy': 0.4, 'Sunny': 0.6}
}

# 出力（観測）確率（emission probabilities）
emission_probability = {
    'Rainy': {'walk': 0.1, 'shop': 0.4, 'clean': 0.5},
    'Sunny': {'walk': 0.6, 'shop': 0.3, 'clean': 0.1}
}

def viterbi(obs, states, start_p, trans_p, emit_p):
    """
    Viterbi アルゴリズムの実装
    obs: 観測列（例: ['walk', 'shop', 'clean']）
    states: 状態のリスト（例: ['Rainy', 'Sunny']）
    start_p: 初期状態確率（辞書）
    trans_p: 状態遷移確率（辞書の辞書）
    emit_p: 観測確率（辞書の辞書）
    """
    V = [{}]  # Viterbiテーブル
    path = {}

    # 初期化
    for s in states:
        V[0][s] = start_p[s] * emit_p[s][obs[0]]
        path[s] = [s]

    # 再帰ステップ
    for t in range(1, len(obs)):
        V.append({})
        new_path = {}

        for curr_state in states:
            (prob, prev_state) = max(
                (V[t - 1][prev_state] * trans_p[prev_state][curr_state] * emit_p[curr_state][obs[t]], prev_state)
                for prev_state in states
            )
            V[t][curr_state] = prob
            new_path[curr_state] = path[prev_state] + [curr_state]

        path = new_path

    # 終了ステップ
    n = len(obs) - 1
    (prob, state) = max((V[n][s], s) for s in states)
    return prob, path[state]

# 実行
prob, path = viterbi(observations, states, start_probability, transition_probability, emission_probability)

print("最尤確率:", prob)
print("最も確からしい状態列:", path)