隠れマルコフモデルにおける能動的推論の実装


# はじめに
- 隠れマルコフモデルにおける能動的推論の実装について解説する
- 観測データに基づく状態推定と、エージェントの能動的な行動選択に焦点を当てる

# 言葉の整理
能動的推論では変分自由エネルギーを最小化することで、隠れマルコフモデルの状態推定と行動選択を同時に最適化します。

それぞれの概念について説明します。


### 隠れマルコフモデル
直接観測できない隠れ状態と、それから生成される観測値の時系列を確率的にモデル化したもの
* 特徴
  * **マルコフ性**：次の状態は現在の状態のみに依存
  * **状態の不可観測性**：真の状態は直接観測できない
  * **確率的生成**：観測値は隠れ状態から確率的に生成される

### 能動的推論
エージェントが自身の行動を選択することで、環境の不確実性を積極的に低減していく推論プロセス
* 特徴：
  * **行動の主体性**：エージェントが能動的に行動を選択
  * **情報獲得**：不確実性を低減するための情報収集
  * **目的指向性**：特定の目標達成のための行動選択

### 変分自由エネルギー

システムの状態推定の質を評価する指標。精度(予測誤差項)とモデルの複雑さの項の和として表現する
* 構成要素：
  * 予測誤差項：モデルの予測と実際の観測の一致度
  * 複雑性項：モデルの複雑さによるペナルティ
* 数式表現：
    * F = -予測誤差項 + 複雑性項
        * 予測誤差項 = $E[log p(o|s)]$
        * 複雑性項 = $KL[q(s)||p(s)]$

複雑性項は、モデルが過度に複雑になることを防ぐために導入されており、事前分布からの逸脱を抑制しデータに過度にフィットすることを防止しています。また、より単純な表現モデルが優先されることを意味しています。

## プロセスフロー

処理フローに整理するとこのようになります。

```mermaid

flowchart LR
    %% ノードの定義
    %% 観測段階のノード
    O1[環境から観測データを取得]
    O2[現在の状態を推定]
    
    %% 行動選択段階のノード
    A1[可能な行動の評価]
    A2[VFEを最小化する行動を選択]
    
    %% 実行段階のノード
    E1[選択された行動の実行]
    E2[新しい観測データの取得]
    
    %% 更新段階のノード
    U1[モデルパラメータの更新]
    U2[状態推定の更新]

%%o[t]: 時刻tでの観測
%%P(o[t]|s[t]): 状態s[t]が与えられた時の観測確率
%%E[情報獲得量]: 期待される情報の増加量
%%コスト: 行動に伴うエネルギー消費など
%%s[t+1]: 次時刻の状態
%%s[t]: 現在の状態
%%a[t]: 選択された行動
%%Match(): 予測誤差の計算関数
%%データ: 実際の観測値
%%予測: モデルによる予測値
%%π: 選択される行動方策
%%argmin[F]: 自由エネルギーFを最小化する引数
%%o[t+1]: 次時刻での観測
%%P(o[t+1]|s[t+1]): 状態から観測への生成モデル
%%q(s): 状態の事後分布
%%P(s|o,π): 観測と行動が与えられた時の状態の条件付き確率
%%予測誤差: データと予測の一致度
%%複雑さ: モデルの複雑性によるペナルティ

    %% サブグラフの定義
    subgraph Observation["観測段階"]
        direction TB
        O1 --> |"o[t] = P(o[t]|s[t])<br>現在の状態から観測を生成"| O2
    end

    subgraph Action["行動選択段階"]
        direction TB
        A1 --> |"行動価値 = E[情報獲得量] - コスト<br>情報獲得と行動コストのバランスを評価"| A2
    end

    subgraph Execution["実行段階"]
        direction TB
        E1 --> |"s[t+1] = P(s[t+1]|s[t],a[t])<br>行動に基づく状態遷移"| E2
    end

    subgraph Update["更新段階"]
        direction TB
        U1 --> |"F = -(予測誤差) + (複雑さ)<br>モデルの評価と更新"| U2
    end

    %% サブグラフ間の関係定義
    Observation --> |"状態推定の精度 = Match(データ, 予測)<br>観測データと予測の一致度を評価"| Action
    
    Action --> |"π = argmin[F]<br>自由エネルギーを最小化する行動を選択"| Execution
    
    Execution --> |"o[t+1] = P(o[t+1]|s[t+1])<br>新しい状態から観測を生成"| Update
    
    Update --> |"q(s) = P(s|o,π)<br>観測と行動に基づき状態分布を更新"| Observation
```

## 実装

状態推定のメカニズムと能動的行動選択の仕組みをPython で構築してみましょう。

```python
import numpy as np
from scipy.special import softmax
from scipy.stats import dirichlet
import matplotlib.pyplot as plt

class ActiveInferenceAgent:
    def __init__(self, n_states, n_observations, n_actions):
        """
        能動的推論エージェントの初期化
        
        Parameters:
        -----------
        n_states : int
            隠れ状態の数
        n_observations : int
            観測可能な状態の数
        n_actions : int
            可能な行動の数
        """
        self.n_states = n_states
        self.n_observations = n_observations
        self.n_actions = n_actions
        
        # モデルパラメータの初期化
        self.A = self._initialize_matrix((n_observations, n_states))  # 観測モデル
        self.B = self._initialize_transition_matrices(n_states, n_actions)  # 状態遷移モデル
        self.C = np.ones(n_observations) / n_observations  # 選好分布
        self.D = np.ones(n_states) / n_states  # 初期状態分布
        
        # 信念状態の初期化
        self.beliefs = np.ones(n_states) / n_states
        
        # ハイパーパラメータ
        self.alpha = 0.1  # 学習率
        self.gamma = 2.0  # 時間割引率
        
    def _initialize_matrix(self, shape):
        """確率行列の初期化"""
        matrix = np.random.rand(*shape)
        return matrix / matrix.sum(axis=0)
    
    def _initialize_transition_matrices(self, n_states, n_actions):
        """行動ごとの状態遷移行列の初期化"""
        B = np.zeros((n_states, n_states, n_actions))
        for a in range(n_actions):
            B[:,:,a] = self._initialize_matrix((n_states, n_states))
        return B
    
    def calculate_free_energy(self, observation, action):
        """
        変分自由エネルギーの計算
        
        Parameters:
        -----------
        observation : int
            現在の観測
        action : int
            評価する行動
            
        Returns:
        --------
        float
            変分自由エネルギー
        """
        # 予測誤差項（accuracy）
        predicted_obs = self.A @ self.beliefs
        accuracy = -np.sum(predicted_obs * np.log(self.C))
        
        # 複雑性項（complexity）
        complexity = np.sum(self.beliefs * np.log(self.beliefs / self.D))
        
        # 行動に依存する状態遷移の評価
        next_state = self.B[:,:,action] @ self.beliefs
        transition_cost = np.sum(next_state * np.log(next_state / self.beliefs))
        
        return accuracy + complexity + transition_cost
    
    def update_beliefs(self, observation):
        """
        観測に基づく信念の更新
        
        Parameters:
        -----------
        observation : int
            現在の観測
        """
        # ベイズ更新
        likelihood = self.A[observation,:]
        self.beliefs = self.beliefs * likelihood
        self.beliefs /= self.beliefs.sum()  # 正規化
    
    def select_action(self, observation):
        """
        行動の選択
        
        Parameters:
        -----------
        observation : int
            現在の観測
            
        Returns:
        --------
        int
            選択された行動
        """
        # 各行動の自由エネルギーを計算
        action_values = np.array([
            self.calculate_free_energy(observation, a)
            for a in range(self.n_actions)
        ])
        
        # ソフトマックス方策で行動を選択
        action_probs = softmax(-self.gamma * action_values)
        return np.random.choice(self.n_actions, p=action_probs)
    
    def update_model(self, observation, action, next_observation):
        """
        モデルパラメータの更新
        
        Parameters:
        -----------
        observation : int
            現在の観測
        action : int
            実行された行動
        next_observation : int
            次の観測
        """
        # 観測モデルの更新
        self.A[observation,:] += self.alpha * self.beliefs
        self.A /= self.A.sum(axis=0)
        
        # 状態遷移モデルの更新
        next_beliefs = self.B[:,:,action] @ self.beliefs
        self.B[:,:,action] += self.alpha * np.outer(next_beliefs, self.beliefs)
        self.B[:,:,action] /= self.B[:,:,action].sum(axis=0)
        
        # 選好分布の更新
        self.C[observation] += self.alpha
        self.C /= self.C.sum()

```



## 使用例

```python
# 使用例
def run_simulation(n_steps=100):
    # エージェントの初期化
    agent = ActiveInferenceAgent(n_states=3, n_observations=4, n_actions=2)
    
    # 結果を保存するリスト
    observations = []
    actions = []
    beliefs = []
    
    # シミュレーション
    observation = np.random.randint(agent.n_observations)
    for t in range(n_steps):
        # 行動の選択と実行
        action = agent.select_action(observation)
        
        # 環境からの新しい観測（ここでは簡単のためにランダム）
        next_observation = np.random.randint(agent.n_observations)
        
        # エージェントの更新
        agent.update_beliefs(observation)
        agent.update_model(observation, action, next_observation)
        
        # 結果の記録
        observations.append(observation)
        actions.append(action)
        beliefs.append(agent.beliefs.copy())
        
        # 次のステップへ
        observation = next_observation
    
    return observations, actions, beliefs

# 結果の可視化
def plot_results(observations, actions, beliefs):
    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 8))
    
    # 観測の推移
    ax1.plot(observations)
    ax1.set_title('Observations')
    ax1.set_ylabel('Observation')
    
    # 行動の推移
    ax2.plot(actions)
    ax2.set_title('Actions')
    ax2.set_ylabel('Action')
    
    # 信念状態の推移
    beliefs_array = np.array(beliefs)
    for i in range(beliefs_array.shape[1]):
        ax3.plot(beliefs_array[:,i], label=f'State {i}')
    ax3.set_title('Beliefs')
    ax3.set_ylabel('Probability')
    ax3.legend()
    
    plt.tight_layout()
    plt.show()

# シミュレーションの実行
observations, actions, beliefs = run_simulation()
plot_results(observations, actions, beliefs)
```

エージェントが環境から観測を受け取り、それに基づいて行動を選択し、その結果として信念を更新する一連のプロセスを可視化しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/348593/2f4885ba-b7cc-4a40-a492-05467f347420.png)

特に信念の変化から、エージェントが学習を通じて特定の状態<State 1（オレンジ）>を排除し、残りの状態間で確率を分配していく様子が分かります。
State 0（青）とState 2（緑）は0.4-0.6の範囲で拮抗。学習の進行に伴い信念が安定化していきます。

# まとめ

能動的推論の基本的な実装を学習しました。


# 用語

| 用語 | 説明 |
|---------|------|
| シグモイド関数 | 入力を0から1の範囲に変換する非線形関数 |
| クロネッカー積 | 2つの行列から得られる積の特殊な形式 |
| ベイズ推論 | 事前確率と観測データから事後確率を推定する手法 |
| one-hot vector | 1つの要素のみが1で他が0のベクトル |
| 遷移行列 | 状態間の遷移確率を表す行列 |
| マッピング行列 | 入力から出力への変換を表現する行列 |
| 隠れマルコフモデル | 観測不可能な状態を含む確率過程モデル |
| 状態空間モデル | システムの内部状態と観測の関係を表現するモデル |
| 事前分布 | データ観測前の確率分布 |
| 事後分布 | データ観測後に更新された確率分布 |
| 対数尤度 | データの生起確率の対数 |
| カテゴリカル分布 | 複数のカテゴリのうち1つが生起する確率分布 |
| ディリクレ分布 | 複数のカテゴリの確率を表現する確率分布 |
| ベータ分布 | 2値の確率を表現する確率分布 |
| 累積分布関数 | 確率変数が特定値以下となる確率を表す関数 |