test_markov_transitions.py

#!/usr/bin/env python3
"""
マルコフ連鎖状態遷移の動作検証スクリプト
設備別カスタマイズ遷移行列によるマルコフ連鎖予測を詳細に検証
"""

import sys
sys.path.append('.')

import yaml
import numpy as np
from collections import defaultdict, Counter
from cbm_environment_v04 import MultiEquipmentCBMEnvironment
import matplotlib.pyplot as plt
import pandas as pd

def test_markov_transitions():
    """マルコフ連鎖による状態遷移の詳細検証"""
    
    print("🔍 マルコフ連鎖状態遷移検証開始")
    print("=" * 60)
    
    # 1. 環境初期化
    config_path = "config_hvac202_v04.yaml"
    with open(config_path, 'r', encoding='utf-8') as f:
        config = yaml.safe_load(f)
    
    env = MultiEquipmentCBMEnvironment(config_path)
    
    print(f"✅ 環境初期化完了")
    print(f"   実データベース遷移行列使用: {env.use_real_data_transitions}")
    print(f"   対象設備数: {env.n_equipment}")
    
    # 2. 各設備・各アクションの理論遷移確率を表示
    print("\n📊 設備別・アクション別理論遷移確率:")
    print("-" * 60)
    
    equipment_names = ['R-1-1', 'R-1-3', 'R-2-2', 'AHU-A-2', 'AHU-B-1', 'AHU-A-1']
    action_names = ['Do Nothing', 'Repair', 'Replace']
    
    theoretical_transitions = {}
    
    for eq_idx in range(env.n_equipment):
        print(f"\n🏭 設備 {eq_idx+1}: {equipment_names[eq_idx]}")
        theoretical_transitions[eq_idx] = {}
        
        for action in range(3):
            if env.use_real_data_transitions:
                trans_matrix = env._get_data_driven_transition(action, eq_idx)
            else:
                base_trans = env.DEFAULT_TRANSITIONS[action] 
                trans_matrix = env._get_age_adjusted_transition(base_trans, eq_idx)
            
            theoretical_transitions[eq_idx][action] = trans_matrix
            
            print(f"   {action_names[action]}:")
            print(f"     Normal → Normal:   {trans_matrix[0,0]:.3f}")
            print(f"     Normal → Anomalous: {trans_matrix[0,1]:.3f}")
            print(f"     Anomalous → Normal:   {trans_matrix[1,0]:.3f}")
            print(f"     Anomalous → Anomalous: {trans_matrix[1,1]:.3f}")
    
    # 3. 実際のマルコフ連鎖遷移を多回数実行して統計検証
    print("\n🎲 マルコフ連鎖遷移統計検証")
    print("-" * 60)
    
    n_trials = 10000  # 多回数実行
    transition_counts = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: Counter())))
    
    # 各設備・各アクション・各初期状態での遷移を記録
    for eq_idx in range(env.n_equipment):
        for action in range(3):
            for initial_state in [0, 1]:  # Normal, Anomalous
                
                state_transitions = []
                
                for trial in range(n_trials):
                    # 環境をリセットして特定の初期状態に設定
                    env.reset()
                    env.equipment_conditions[eq_idx] = initial_state
                    
                    # 特定設備のみにアクションを適用するマスクを作成
                    full_action = np.zeros(env.n_equipment, dtype=int)
                    full_action[eq_idx] = action
                    
                    # アクションをトリプレットインデックスに変換
                    action_idx = 0
                    for i, a in enumerate(full_action):
                        action_idx += a * (3 ** (env.n_equipment - 1 - i))
                    
                    # ステップ実行
                    old_state = env.equipment_conditions[eq_idx]
                    env.step(action_idx)
                    new_state = env.equipment_conditions[eq_idx]
                    
                    state_transitions.append((old_state, new_state))
                    transition_counts[eq_idx][action][initial_state][new_state] += 1
                
    # 4. 理論値と実測値の比較
    print("\n📈 理論値 vs 実測値比較:")
    print("-" * 60)
    
    verification_results = []
    
    for eq_idx in range(env.n_equipment):
        print(f"\n🏭 設備 {eq_idx+1}: {equipment_names[eq_idx]}")
        
        for action in range(3):
            print(f"\n   {action_names[action]}:")
            
            for initial_state in [0, 1]:
                state_name = "Normal" if initial_state == 0 else "Anomalous"
                
                # 実測遷移確率の計算
                total = sum(transition_counts[eq_idx][action][initial_state].values())
                if total > 0:
                    empirical_probs = {
                        next_state: count / total 
                        for next_state, count in transition_counts[eq_idx][action][initial_state].items()
                    }
                    
                    # 理論確率
                    theoretical_matrix = theoretical_transitions[eq_idx][action]
                    theoretical_probs = {
                        0: theoretical_matrix[initial_state, 0],
                        1: theoretical_matrix[initial_state, 1]
                    }
                    
                    print(f"     {state_name} → Normal:")
                    print(f"       理論値: {theoretical_probs[0]:.3f}")
                    print(f"       実測値: {empirical_probs.get(0, 0):.3f}")
                    print(f"       誤差: {abs(theoretical_probs[0] - empirical_probs.get(0, 0)):.3f}")
                    
                    print(f"     {state_name} → Anomalous:")
                    print(f"       理論値: {theoretical_probs[1]:.3f}")
                    print(f"       実測値: {empirical_probs.get(1, 0):.3f}")
                    print(f"       誤差: {abs(theoretical_probs[1] - empirical_probs.get(1, 0)):.3f}")
                    
                    # 検証結果記録
                    verification_results.append({
                        'equipment': eq_idx,
                        'action': action,
                        'initial_state': initial_state,
                        'transition_to_normal_theory': theoretical_probs[0],
                        'transition_to_normal_empirical': empirical_probs.get(0, 0),
                        'transition_to_anomalous_theory': theoretical_probs[1],
                        'transition_to_anomalous_empirical': empirical_probs.get(1, 0)
                    })
    
    # 5. 統計的精度の検証
    print("\n📊 マルコフ連鎖精度統計:")
    print("-" * 60)
    
    df_results = pd.DataFrame(verification_results)
    
    # Normal遷移の誤差分析
    normal_errors = abs(df_results['transition_to_normal_theory'] - df_results['transition_to_normal_empirical'])
    anomalous_errors = abs(df_results['transition_to_anomalous_theory'] - df_results['transition_to_anomalous_empirical'])
    
    print(f"Normal遷移誤差:")
    print(f"  平均誤差: {normal_errors.mean():.4f}")
    print(f"  最大誤差: {normal_errors.max():.4f}")
    print(f"  標準偏差: {normal_errors.std():.4f}")
    
    print(f"\nAnomalous遷移誤差:")
    print(f"  平均誤差: {anomalous_errors.mean():.4f}")
    print(f"  最大誤差: {anomalous_errors.max():.4f}")
    print(f"  標準偏差: {anomalous_errors.std():.4f}")
    
    # 6. 結論
    print("\n" + "=" * 60)
    max_error_threshold = 0.02  # 2%以内なら十分精密
    
    if normal_errors.max() < max_error_threshold and anomalous_errors.max() < max_error_threshold:
        print("🎉 マルコフ連鎖状態遷移検証成功！")
        print(f"   設備別カスタマイズ遷移行列による正確なマルコフ連鎖が動作中")
        print(f"   最大誤差 {max(normal_errors.max(), anomalous_errors.max()):.4f} < 閾値 {max_error_threshold}")
    else:
        print("⚠️ マルコフ連鎖に精度の問題があります")
        print(f"   最大誤差 {max(normal_errors.max(), anomalous_errors.max()):.4f} >= 閾値 {max_error_threshold}")
    
    print("=" * 60)
    
    return verification_results


if __name__ == "__main__":
    test_markov_transitions()