from enum import Enum
from dataclasses import dataclass
import numpy as np

class RobotState(Enum):
    """机器人状态空间"""
    IDLE = "空闲"      # 静止等待
    MOVING = "移动"    # 正在移动
    CHARGING = "充电"  # 正在充电

@dataclass
class MarkovChain:
    """一阶马尔科夫链实现"""
    
    # 状态转移矩阵
    # P[next_state | current_state]
    transition_matrix = {
        RobotState.IDLE: {
            RobotState.IDLE: 0.5,     # 保持空闲
            RobotState.MOVING: 0.4,   # 开始移动
            RobotState.CHARGING: 0.1  # 开始充电
        },
        RobotState.MOVING: {
            RobotState.IDLE: 0.3,     # 停止移动
            RobotState.MOVING: 0.5,   # 继续移动
            RobotState.CHARGING: 0.2  # 开始充电
        },
        RobotState.CHARGING: {
            RobotState.IDLE: 0.8,     # 充满电，空闲
            RobotState.MOVING: 0.1,   # 充满电，开始移动
            RobotState.CHARGING: 0.1  # 继续充电
        }
    }
    
    def next_state(self, current: RobotState) -> RobotState:
        """基于当前状态生成下一个状态"""
        import random
        
        # 获取当前状态的所有可能转移
        transitions = self.transition_matrix[current]
        
        # 随机选择（按概率权重）
        states = list(transitions.keys())
        weights = list(transitions.values())
        
        return random.choices(states, weights=weights)[0]

import numpy as np

# 状态顺序：[空闲, 移动, 充电]
# 行：当前状态，列：下一状态
P = np.array([
    [0.5, 0.4, 0.1],  # 空闲 → [空闲, 移动, 充电]
    [0.3, 0.5, 0.2],  # 移动 → [空闲, 移动, 充电]
    [0.8, 0.1, 0.1]   # 充电 → [空闲, 移动, 充电]
])

# 关键性质：每行和为1（概率归一化）
print("行和验证:", np.sum(P, axis=1))  # [1., 1., 1.]

# 一阶模型：P(下雨|今天=晴) = 0.2
# 问题：连续10天晴天后，下雨概率还是0.2吗？

# 三阶模型更准确：
weather_probs = {
    # (前前天, 前天, 今天) → 明天天气概率
    ('晴', '晴', '晴'): {'雨': 0.6, '晴': 0.4},  # 长期晴天后更可能下雨
    ('雨', '晴', '晴'): {'雨': 0.3, '晴': 0.7},
    ('晴', '雨', '晴'): {'雨': 0.4, '晴': 0.6},
}

class HigherOrderMarkovChain:
    """τ阶马尔科夫链通用实现"""
    
    def __init__(self, order: int):
        self.order = order  # 记忆长度
        self.memory = []    # 存储最近order个状态
        
        # 转移概率表：P(X_t | X_{t-τ}, ..., X_{t-1})
        self.transitions = {}
    
    def add_transition(self, history: tuple, next_state: str, prob: float):
        """添加转移概率"""
        if history not in self.transitions:
            self.transitions[history] = {}
        self.transitions[history][next_state] = prob
    
    def predict(self) -> str:
        """基于历史预测下一个状态"""
        if len(self.memory) < self.order:
            return None
        
        # 获取最近的order个状态作为历史
        recent_history = tuple(self.memory[-self.order:])
        
        if recent_history in self.transitions:
            # 根据概率随机选择
            probs = self.transitions[recent_history]
            import random
            return random.choices(list(probs.keys()), 
                                 weights=probs.values())[0]
        return None

def convert_to_first_order(states: list, high_order_probs: dict, order: int):
    """
    将高阶马尔科夫链转换为等价的一阶链
    
    思想：将长度为order的历史序列视为一个"复合状态"
    例如：二阶链的(S,M)视为一个新状态
    """
    
    # 1. 创建所有可能的复合状态
    composite_states = []
    from itertools import product
    
    # 生成所有长度为order的状态序列
    for combo in product(states, repeat=order):
        composite_states.append(combo)
    
    # 2. 构建一阶转移矩阵
    n_composite = len(composite_states)
    P_first_order = np.zeros((n_composite, n_composite))
    
    # 3. 填充转移概率
    composite_index = {cs: i for i, cs in enumerate(composite_states)}
    
    for history_tuple, next_probs in high_order_probs.items():
        i = composite_index[history_tuple]
        
        for next_state, prob in next_probs.items():
            # 新历史：(移出最旧状态，加入新状态)
            # 例如：(S,M) + M → (M,M)
            new_history = history_tuple[1:] + (next_state,)
            j = composite_index[new_history]
            P_first_order[i, j] = prob
    
    return composite_states, P_first_order

# 示例：二阶链转换
states = ['S', 'M', 'C']
second_order_probs = {
    ('S', 'S'): {'S': 0.6, 'M': 0.3, 'C': 0.1},
    ('S', 'M'): {'S': 0.2, 'M': 0.6, 'C': 0.2},
    ('M', 'S'): {'S': 0.4, 'M': 0.5, 'C': 0.1},
    # ... 其他组合
}

composite_states, P_1st = convert_to_first_order(states, second_order_probs, order=2)
print(f"原始状态数: {len(states)}")
print(f"复合状态数: {len(composite_states)}")  # 3² = 9

import matplotlib.pyplot as plt
import networkx as nx
from matplotlib.patches import FancyBboxPatch

def visualize_markov_chains():
    """改进的可视化：二阶马尔科夫链转换为一阶链"""
    
    # 创建图形
    fig, axes = plt.subplots(1, 2, figsize=(16, 8))
    fig.suptitle('Transforming Second-Order Markov Chain to First-Order', 
                fontsize=16, fontweight='bold', y=0.95)
    
    # ========== 左图：二阶马尔科夫链 ==========
    ax1 = axes[0]
    ax1.set_title('Second-Order Markov Chain\n(Needs memory of past 2 days)', 
                 fontsize=14, fontweight='bold', pad=20)
    ax1.set_xlim(-1, 11)
    ax1.set_ylim(-1, 11)
    ax1.axis('off')
    
    # 时间线标签
    time_labels = ['Day -2', 'Day -1', 'Today']
    times = [2, 5, 8]
    
    # 绘制时间线
    ax1.plot([1.5, 8.5], [8, 8], 'k-', linewidth=2, alpha=0.7)
    
    for i, (time, label) in enumerate(zip(times, time_labels)):
        ax1.text(time, 8.3, label, ha='center', fontsize=11, 
                fontweight='bold', color='darkblue')
        # 时间点标记
        ax1.plot(time, 8, 'ko', markersize=10)
    
    # 天气示例：Sunny, Sunny, Cloudy
    weather_sequence = ['Sunny', 'Sunny', 'Cloudy']
    weather_icons = {'Sunny': '☀️', 'Cloudy': '☁️', 'Rainy': '🌧️'}
    
    for i, (time, weather) in enumerate(zip(times, weather_sequence)):
        ax1.text(time, 7.5, weather_icons[weather], fontsize=40, ha='center')
        ax1.text(time, 7.0, weather, ha='center', fontsize=12, 
                fontweight='bold', color='darkblue')
    
    # 依赖箭头
    ax1.arrow(times[0], 6.8, times[1]-times[0]-0.3, -0.7, 
              head_width=0.15, head_length=0.2, fc='red', ec='red', alpha=0.7)
    ax1.arrow(times[1], 6.8, times[2]-times[1]-0.3, -0.7, 
              head_width=0.15, head_length=0.2, fc='red', ec='red', alpha=0.7)
    
    # 状态空间说明
    state_space_box = FancyBboxPatch((1, 3), 8, 2,
                                    boxstyle="round,pad=0.3",
                                    facecolor="lightcoral", alpha=0.2,
                                    edgecolor="red", linewidth=1.5)
    ax1.add_patch(state_space_box)
    
    ax1.text(5, 4.5, 'State Space Size Problem', 
            ha='center', fontsize=12, fontweight='bold', color='darkred')
    ax1.text(5, 3.8, '3 weather types → 3² = 9 possible history pairs', 
            ha='center', fontsize=10, color='darkred')
    ax1.text(5, 3.3, 'P(Today | Yesterday, Day-before-yesterday)', 
            ha='center', fontsize=10, color='darkred', style='italic')
    
    # 内存需求
    memory_box = FancyBboxPatch((1, 0.5), 8, 1.5,
                               boxstyle="round,pad=0.3",
                               facecolor="lightblue", alpha=0.2,
                               edgecolor="blue", linewidth=1.5)
    ax1.add_patch(memory_box)
    
    ax1.text(5, 1.7, 'Memory Requirement', 
            ha='center', fontsize=12, fontweight='bold', color='darkblue')
    ax1.text(5, 1.0, 'Need to remember last 2 states', 
            ha='center', fontsize=10, color='darkblue')
    
    # ========== 右图：转换为一阶链 ==========
    ax2 = axes[1]
    ax2.set_title('Equivalent First-Order Markov Chain\n(Composite states)', 
                 fontsize=14, fontweight='bold', pad=20)
    ax2.set_xlim(-1, 11)
    ax2.set_ylim(-1, 11)
    ax2.axis('off')
    
    # 复合状态的概念
    composite_box = FancyBboxPatch((1, 8), 8, 2,
                                  boxstyle="round,pad=0.3",
                                  facecolor="lightgreen", alpha=0.2,
                                  edgecolor="green", linewidth=2)
    ax2.add_patch(composite_box)
    
    ax2.text(5, 9.5, 'Composite State = Memory Encoded in State Name', 
            ha='center', fontsize=12, fontweight='bold', color='darkgreen')
    
    ax2.text(5, 8.8, 'Instead of: P(Weather_today | Weather_yesterday, Weather_day-before)', 
            ha='center', fontsize=9, color='darkgreen')
    ax2.text(5, 8.2, 'We use: P(Composite_today | Composite_yesterday)', 
            ha='center', fontsize=9, color='darkgreen')
    
    # 复合状态示例
    ax2.text(3, 7.0, 'Composite State Example:', ha='left', fontsize=11, fontweight='bold')
    
    # 状态分解图示
    # 复合状态
    comp_state = FancyBboxPatch((3, 5.5), 4, 1,
                               boxstyle="round,pad=0.3",
                               facecolor="lightblue", alpha=0.3)
    ax2.add_patch(comp_state)
    ax2.text(5, 6.0, '"Sunny-Sunny"', ha='center', fontsize=12, 
            fontweight='bold', color='darkblue')
    ax2.text(5, 5.5, 'represents: Sunny yesterday + Sunny day-before', 
            ha='center', fontsize=9, color='blue')
    
    # 箭头到天气
    ax2.arrow(5, 5.2, 0, -1, head_width=0.2, head_length=0.15, 
              fc='purple', ec='purple', alpha=0.7, linestyle='--')
    
    # 对应天气
    weather_box = FancyBboxPatch((3, 3.5), 4, 1,
                                boxstyle="round,pad=0.3",
                                facecolor="yellow", alpha=0.2)
    ax2.add_patch(weather_box)
    ax2.text(5, 4.0, 'Actual Weather Today', ha='center', fontsize=10, fontweight='bold')
    ax2.text(5, 3.5, 'Sunny', ha='center', fontsize=14, fontweight='bold', color='darkorange')
    ax2.text(5, 3.5, '☀️', fontsize=30, ha='center')
    
    # 状态转移示例
    ax2.text(3, 2.5, 'State Transition Example:', ha='left', fontsize=11, fontweight='bold')
    
    # 转移箭头
    ax2.plot([2, 8], [2, 2], 'k-', linewidth=1, alpha=0.5)
    
    # 从状态 (Sunny, Sunny)
    state1_box = FancyBboxPatch((1.5, 1.2), 2, 1,
                               boxstyle="round,pad=0.3",
                               facecolor="lightblue", alpha=0.4)
    ax2.add_patch(state1_box)
    ax2.text(2.5, 1.7, 'State: (S,S)', ha='center', fontsize=10, fontweight='bold')
    
    # 转移箭头
    ax2.arrow(3.5, 1.7, 2, 0, head_width=0.15, head_length=0.2, 
              fc='red', ec='red', alpha=0.7)
    
    # 到状态 (Sunny, Cloudy)
    state2_box = FancyBboxPatch((5.5, 1.2), 2, 1,
                               boxstyle="round,pad=0.3",
                               facecolor="lightblue", alpha=0.4)
    ax2.add_patch(state2_box)
    ax2.text(6.5, 1.7, 'State: (S,C)', ha='center', fontsize=10, fontweight='bold')
    
    # 转移概率说明
    prob_text = 'Transition Probability:\nP((S,C) | (S,S)) = 0.4'
    ax2.text(6.5, 0.5, prob_text, ha='center', fontsize=9, 
            bbox=dict(boxstyle="round,pad=0.2", facecolor="yellow", alpha=0.3))
    
    # 关键优势
    advantage_box = FancyBboxPatch((1, -0.5), 8, 1.5,
                                  boxstyle="round,pad=0.3",
                                  facecolor="gold", alpha=0.2,
                                  edgecolor="orange", linewidth=1.5)
    ax2.add_patch(advantage_box)
    
    ax2.text(5, 0.0, 'Key Advantage: Standard Markov techniques apply', 
            ha='center', fontsize=11, fontweight='bold', color='darkorange')
    ax2.text(5, -0.5, 'State space: 9 composite states instead of complex memory', 
            ha='center', fontsize=9, color='darkorange')
    
    plt.tight_layout()
    plt.show()

visualize_markov_chains()

import random

def estimate_pi(num_samples):
    count_inside = 0
    for _ in range(num_samples):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        if x**2 + y**2 <= 1:
            count_inside += 1
    return 4 * count_inside / num_samples

# 测试
for n in [100, 1000, 10000, 100000]:
    pi_est = estimate_pi(n)
    print(f"N={n}: π ≈ {pi_est}, 误差={abs(pi_est - 3.1415926535)}")

import random
import math

def mc_integral(func, a, b, num_samples):
    total = 0.0
    for _ in range(num_samples):
        x = random.uniform(a, b)
        total += func(x)
    return (b - a) * total / num_samples

def f(x):
    return math.exp(-x**2)

# 积分精确值可以用特殊函数算，这里仅用蒙特卡洛估计
for n in [1000, 10000, 100000]:
    result = mc_integral(f, 0, 1, n)
    print(f"N={n}: 积分 ≈ {result}")

import numpy as np
import matplotlib.pyplot as plt
from collections import defaultdict

class CliffWalkingEnv:
    """简化的悬崖行走环境"""
    
    def __init__(self):
        self.width = 12
        self.height = 4
        self.start = (3, 0)
        self.goal = (3, 11)
        self.cliff = [(3, i) for i in range(1, 11)]
        
        self.actions = [0, 1, 2, 3]  # 上下左右
        self.reset()
    
    def reset(self):
        self.state = self.start
        return self.state
    
    def step(self, action):
        x, y = self.state
        
        # 移动
        if action == 0: x = max(x-1, 0)      # 上
        elif action == 1: y = min(y+1, 11)   # 右
        elif action == 2: x = min(x+1, 3)    # 下
        elif action == 3: y = max(y-1, 0)    # 左
        
        self.state = (x, y)
        
        # 奖励
        if self.state in self.cliff:
            return self.start, -100, True
        elif self.state == self.goal:
            return self.state, 0, True
        else:
            return self.state, -1, False
    
    def show(self, title=""):
        """显示环境"""
        grid = []
        for i in range(4):
            row = []
            for j in range(12):
                if (i, j) == self.start:
                    row.append('S')
                elif (i, j) == self.goal:
                    row.append('G')
                elif (i, j) in self.cliff:
                    row.append('C')
                else:
                    row.append('.')
            grid.append(' '.join(row))
        print(f"\n{title}")
        print('\n'.join(grid))


class SarsaAgent:
    """SARSA智能体"""
    
    def __init__(self):
        self.Q = defaultdict(lambda: np.zeros(4))
        self.alpha = 0.1  # 学习率
        self.gamma = 0.9  # 折扣因子
        self.epsilon = 0.1  # 探索率
    
    def choose_action(self, state):
        """ε-greedy策略"""
        if np.random.random() < self.epsilon:
            return np.random.randint(4)  # 探索
        else:
            return np.argmax(self.Q[state])  # 利用
    
    def update(self, s, a, r, s_next, a_next):
        """SARSA更新"""
        target = r + self.gamma * self.Q[s_next][a_next]
        self.Q[s][a] += self.alpha * (target - self.Q[s][a])


class QLearningAgent:
    """Q-Learning智能体"""
    
    def __init__(self):
        self.Q = defaultdict(lambda: np.zeros(4))
        self.alpha = 0.1
        self.gamma = 0.9
        self.epsilon = 0.1
    
    def choose_action(self, state):
        """ε-greedy策略"""
        if np.random.random() < self.epsilon:
            return np.random.randint(4)
        else:
            return np.argmax(self.Q[state])
    
    def update(self, s, a, r, s_next):
        """Q-Learning更新"""
        max_q = np.max(self.Q[s_next])
        target = r + self.gamma * max_q
        self.Q[s][a] += self.alpha * (target - self.Q[s][a])


def train(env, agent, episodes=200):
    """训练函数"""
    rewards = []
    
    for episode in range(episodes):
        s = env.reset()
        
        if isinstance(agent, SarsaAgent):
            a = agent.choose_action(s)  # SARSA需要先选动作
        
        total_reward = 0
        done = False
        
        while not done:
            if isinstance(agent, QLearningAgent):
                a = agent.choose_action(s)  # Q-Learning在循环内选动作
            
            s_next, r, done = env.step(a)
            total_reward += r
            
            if isinstance(agent, SarsaAgent):
                a_next = agent.choose_action(s_next)
                agent.update(s, a, r, s_next, a_next)
                a = a_next  # 使用实际要执行的动作
            else:
                agent.update(s, a, r, s_next)
                a = agent.choose_action(s_next)  # 重新选择动作
            
            s = s_next
            
            if done:
                break
        
        rewards.append(total_reward)
        
        # 逐渐减少探索
        if episode % 100 == 0 and episode > 0:
            agent.epsilon *= 0.9
    
    return rewards


def compare():
    """基础对比SARSA和Q-Learning"""
    print("悬崖行走环境对比")
    print("S:起点, G:终点, C:悬崖")
    print("-" * 50)
    
    env = CliffWalkingEnv()
    env.show("悬崖行走环境")
    
    # 训练两个智能体
    print("\n训练SARSA...")
    sarsa_agent = SarsaAgent()
    sarsa_rewards = train(env, sarsa_agent, 500)
    
    print("训练Q-Learning...")
    qlearning_agent = QLearningAgent()
    qlearning_rewards = train(env, qlearning_agent, 500)
    
    # 可视化结果
    plt.figure(figsize=(12, 4))
    
    # 原始奖励曲线
    plt.subplot(1, 3, 1)
    plt.plot(sarsa_rewards, label='SARSA', alpha=0.6, linewidth=0.5)
    plt.plot(qlearning_rewards, label='Q-Learning', alpha=0.6, linewidth=0.5)
    plt.xlabel('Episode')
    plt.ylabel('Reward')
    plt.title('Raw Rewards')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    # 滑动平均（更平滑）
    plt.subplot(1, 3, 2)
    window = 20
    sarsa_smooth = np.convolve(sarsa_rewards, np.ones(window)/window, mode='valid')
    qlearning_smooth = np.convolve(qlearning_rewards, np.ones(window)/window, mode='valid')
    plt.plot(sarsa_smooth, label='SARSA', linewidth=2)
    plt.plot(qlearning_smooth, label='Q-Learning', linewidth=2)
    plt.xlabel('Episode')
    plt.ylabel('Smoothed Reward')
    plt.title(f'Moving Average (window={window})')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    # 最后100个episode的分布
    plt.subplot(1, 3, 3)
    plt.boxplot([sarsa_rewards[-100:], qlearning_rewards[-100:]], 
                labels=['SARSA', 'Q-Learning'])
    plt.ylabel('Reward')
    plt.title('Reward Distribution (Last 100 Episodes)')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # 打印统计信息
    print("\n" + "="*50)
    print("结果对比:")
    print(f"SARSA 平均奖励（最后100次）: {np.mean(sarsa_rewards[-100:]):.1f}")
    print(f"Q-Learning 平均奖励（最后100次）: {np.mean(qlearning_rewards[-100:]):.1f}")
    
    print(f"\nSARSA 最低奖励（最后100次）: {np.min(sarsa_rewards[-100:]):.1f}")
    print(f"Q-Learning 最低奖励（最后100次）: {np.min(qlearning_rewards[-100:]):.1f}")
    
    # 计算掉下悬崖的次数
    sarsa_falls = sum(1 for r in sarsa_rewards[-100:] if r <= -100)
    qlearning_falls = sum(1 for r in qlearning_rewards[-100:] if r <= -100)
    
    print(f"\nSARSA 掉下悬崖次数（最后100次）: {sarsa_falls}")
    print(f"Q-Learning 掉下悬崖次数（最后100次）: {qlearning_falls}")
    
    # 分析原因
    print("\n" + "="*50)
    print("核心区别分析:")
    print("1. SARSA (On-Policy):")
    print("   - 更新时使用实际执行的下一个动作")
    print("   - 考虑探索风险，学习安全路径")
    print("   - 平均奖励较高且稳定")
    
    print("\n2. Q-Learning (Off-Policy):")
    print("   - 更新时使用理论最优的下一个动作")
    print("   - 学习最短但危险的路径")
    print("   - 探索时容易掉下悬崖（奖励-100）")
    print("   - 最低奖励较低")


def compare_detailed():
    """更详细的对比（英文标签）"""
    print("悬崖行走环境对比 - Cliff Walking Environment Comparison")
    print("S: Start, G: Goal, C: Cliff")
    print("-" * 60)
    
    env = CliffWalkingEnv()
    env.show("Cliff Walking Environment")
    
    # 运行多次实验以获得更稳定的结果
    n_runs = 5
    all_sarsa_rewards = []
    all_qlearning_rewards = []
    
    print(f"\nRunning {n_runs} experiments for each algorithm...")
    
    for run in range(n_runs):
        print(f"Run {run+1}/{n_runs}...")
        
        # SARSA
        sarsa_agent = SarsaAgent()
        sarsa_rewards = train(env, sarsa_agent, 200)
        all_sarsa_rewards.append(sarsa_rewards)
        
        # Q-Learning
        qlearning_agent = QLearningAgent()
        qlearning_rewards = train(env, qlearning_agent, 200)
        all_qlearning_rewards.append(qlearning_rewards)
    
    # 转换为numpy数组便于计算
    all_sarsa_rewards = np.array(all_sarsa_rewards)
    all_qlearning_rewards = np.array(all_qlearning_rewards)
    
    # 计算平均值和标准差
    sarsa_mean = np.mean(all_sarsa_rewards, axis=0)
    qlearning_mean = np.mean(all_qlearning_rewards, axis=0)
    sarsa_std = np.std(all_sarsa_rewards, axis=0)
    qlearning_std = np.std(all_qlearning_rewards, axis=0)
    
    # 创建更专业的图表
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. 平均奖励曲线
    ax1 = axes[0, 0]
    episodes = range(len(sarsa_mean))
    ax1.plot(episodes, sarsa_mean, 'b-', label='SARSA Mean', linewidth=2)
    ax1.fill_between(episodes, sarsa_mean - sarsa_std, sarsa_mean + sarsa_std, 
                    alpha=0.2, color='blue')
    ax1.plot(episodes, qlearning_mean, 'r-', label='Q-Learning Mean', linewidth=2)
    ax1.fill_between(episodes, qlearning_mean - qlearning_std, qlearning_mean + qlearning_std, 
                    alpha=0.2, color='red')
    ax1.set_xlabel('Episode')
    ax1.set_ylabel('Average Reward')
    ax1.set_title('Learning Curves with Standard Deviation')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 2. 滑动平均对比
    ax2 = axes[0, 1]
    window = 20
    sarsa_smooth = np.convolve(sarsa_mean, np.ones(window)/window, mode='valid')
    qlearning_smooth = np.convolve(qlearning_mean, np.ones(window)/window, mode='valid')
    ax2.plot(sarsa_smooth, 'b-', label='SARSA', linewidth=2)
    ax2.plot(qlearning_smooth, 'r-', label='Q-Learning', linewidth=2)
    ax2.set_xlabel('Episode')
    ax2.set_ylabel('Smoothed Reward')
    ax2.set_title(f'Smoothed Learning Curves (Moving Average, window={window})')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 3. 最终性能分布
    ax3 = axes[1, 0]
    # 修复：使用np.mean计算平均值
    final_sarsa = [np.mean(run[-50:]) for run in all_sarsa_rewards]
    final_qlearning = [np.mean(run[-50:]) for run in all_qlearning_rewards]
    
    positions = [1, 2]
    ax3.boxplot([final_sarsa, final_qlearning], positions=positions, 
                labels=['SARSA', 'Q-Learning'], widths=0.6)
    ax3.set_ylabel('Average Reward (Last 50 Episodes)')
    ax3.set_title('Final Performance Distribution')
    ax3.grid(True, alpha=0.3)
    
    # 添加均值标记
    ax3.text(1, np.mean(final_sarsa), f'Mean: {np.mean(final_sarsa):.1f}', 
            ha='center', va='bottom', fontweight='bold')
    ax3.text(2, np.mean(final_qlearning), f'Mean: {np.mean(final_qlearning):.1f}', 
            ha='center', va='bottom', fontweight='bold')
    
    # 4. 算法特性对比
    ax4 = axes[1, 1]
    properties = ['Update Method', 'Policy Type', 'Risk Level', 'Convergence Speed']
    sarsa_scores = [1.0, 1.0, 0.8, 0.7]  # 相对评分
    qlearning_scores = [1.0, 0.6, 0.4, 0.9]
    
    x = np.arange(len(properties))
    width = 0.35
    
    ax4.bar(x - width/2, sarsa_scores, width, label='SARSA', color='blue', alpha=0.7)
    ax4.bar(x + width/2, qlearning_scores, width, label='Q-Learning', color='red', alpha=0.7)
    
    ax4.set_xlabel('Property')
    ax4.set_ylabel('Score (Relative)')
    ax4.set_title('Algorithm Properties Comparison')
    ax4.set_xticks(x)
    ax4.set_xticklabels(properties, rotation=15)
    ax4.set_ylim(0, 1.2)
    ax4.legend()
    ax4.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.show()
    
    # 打印中文分析
    print("\n" + "="*60)
    print("详细分析结果:")
    print(f"SARSA 最终平均奖励（最后50步）: {np.mean(final_sarsa):.1f} ± {np.std(final_sarsa):.1f}")
    print(f"Q-Learning 最终平均奖励（最后50步）: {np.mean(final_qlearning):.1f} ± {np.std(final_qlearning):.1f}")
    
    print("\n关键观察:")
    print("1. SARSA的奖励更稳定（标准差较小）")
    print("2. Q-Learning在某些运行中可能表现更好，但方差较大")
    print("3. 在收敛速度上，Q-Learning通常更快")
    print("4. SARSA在安全性方面表现更好")


# 运行对比
if __name__ == "__main__":
    print("选择对比模式:")
    print("1. 基础对比 (Basic Comparison)")
    print("2. 详细对比 (Detailed Comparison)")
    choice = input("请输入选择 (1或2): ").strip()
    
    if choice == "2":
        compare_detailed()
    else:
        compare()

--------------------------------------------------

悬崖行走环境
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
S C C C C C C C C C C G

训练SARSA...
训练Q-Learning...
<ipython-input-1-e8b610651e2d>:193: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  plt.boxplot([sarsa_rewards[-100:], qlearning_rewards[-100:]],

==================================================
结果对比:
SARSA 平均奖励（最后100次）: -15.6
Q-Learning 平均奖励（最后100次）: -27.7

SARSA 最低奖励（最后100次）: -24.0
Q-Learning 最低奖励（最后100次）: -115.0

SARSA 掉下悬崖次数（最后100次）: 0
Q-Learning 掉下悬崖次数（最后100次）: 15

==================================================
核心区别分析:
1. SARSA (On-Policy):
   - 更新时使用实际执行的下一个动作
   - 考虑探索风险，学习安全路径
   - 平均奖励较高且稳定

2. Q-Learning (Off-Policy):
   - 更新时使用理论最优的下一个动作
   - 学习最短但危险的路径
   - 探索时容易掉下悬崖（奖励-100）
   - 最低奖励较低

import numpy as np
import matplotlib.pyplot as plt

class QLearningWithBias:
    """带有过估计的Q-Learning实现"""
    
    def __init__(self, n_actions, alpha=0.1, gamma=0.9, epsilon=0.1, 
                 estimation_noise=0.5, use_double_q=False):
        self.n_actions = n_actions
        self.alpha = alpha
        self.gamma = gamma
        self.epsilon = epsilon
        self.estimation_noise = estimation_noise  # 估计噪声大小
        self.use_double_q = use_double_q
        
        # Q表
        self.Q = {}
        if use_double_q:
            self.Q2 = {}  # Double Q-Learning需要两个Q表
    
    def get_q_value(self, state, action):
        """获取Q值，加入模拟的估计噪声"""
        key = (state, action)
        if key not in self.Q:
            self.Q[key] = 0.0
        
        # 模拟估计噪声
        true_q = self.Q[key]
        if self.estimation_noise > 0:
            noise = np.random.normal(0, self.estimation_noise)
            return true_q + noise
        return true_q
    
    def choose_action(self, state):
        """ε-greedy策略"""
        if np.random.random() < self.epsilon:
            return np.random.randint(self.n_actions)
        else:
            q_values = [self.get_q_value(state, a) for a in range(self.n_actions)]
            return np.argmax(q_values)
    
    def update(self, state, action, reward, next_state, done):
        """更新Q值"""
        current_q = self.get_q_value(state, action)
        
        if done:
            target = reward
        else:
            if self.use_double_q:
                # Double Q-Learning更新
                # 1. 用Q1选择动作
                q1_values = [self.get_q_value(next_state, a) for a in range(self.n_actions)]
                best_action = np.argmax(q1_values)
                
                # 2. 用Q2评估动作
                if (next_state, best_action) not in self.Q2:
                    self.Q2[(next_state, best_action)] = 0.0
                target = reward + self.gamma * self.Q2[(next_state, best_action)]
            else:
                # 标准Q-Learning更新（容易过估计）
                q_values = [self.get_q_value(next_state, a) for a in range(self.n_actions)]
                target = reward + self.gamma * np.max(q_values)
        
        # 更新真正的Q值（不加噪声）
        key = (state, action)
        if key not in self.Q:
            self.Q[key] = 0.0
        self.Q[key] += self.alpha * (target - self.Q[key])
        
        # 如果是Double Q-Learning，也更新第二个Q表
        if self.use_double_q:
            if np.random.random() < 0.5:  # 随机选择一个Q表更新
                self.update_q2(state, action, reward, next_state, done)
    
    def update_q2(self, state, action, reward, next_state, done):
        """更新第二个Q表（用于Double Q-Learning）"""
        current_q = self.Q2.get((state, action), 0.0)
        
        if done:
            target = reward
        else:
            # 用Q2选择动作
            q2_values = [self.Q2.get((next_state, a), 0.0) for a in range(self.n_actions)]
            best_action = np.argmax(q2_values)
            
            # 用Q1评估动作
            target = reward + self.gamma * self.Q.get((next_state, best_action), 0.0)
        
        self.Q2[(state, action)] = current_q + self.alpha * (target - current_q)
    
    def get_true_q_values(self, state):
        """获取无噪声的真实Q值（用于分析）"""
        return [self.Q.get((state, a), 0.0) for a in range(self.n_actions)]

def compare_overestimation():
    """比较Q-Learning和Double Q-Learning的过估计"""
    print("\n=== Q-Learning vs Double Q-Learning 过估计对比 ===")
    print("修正版本：更好地展示Double Q-Learning的优势")
    
    # 创建一个更有结构的环境设置
    n_states = 50  # 增加状态数量
    n_actions = 6  # 增加动作数量，更容易产生过估计
    
    # 运行多次实验
    n_experiments = 50  # 增加实验次数
    qlearning_overestimations = []
    doubleq_overestimations = []
    
    for exp in range(n_experiments):
        # Q-Learning - 增加噪声以放大过估计效果
        ql_agent = QLearningWithBias(n_actions, estimation_noise=2.0, use_double_q=False)
        
        # Double Q-Learning - 相同的噪声设置
        dql_agent = QLearningWithBias(n_actions, estimation_noise=2.0, use_double_q=True)
        
        # 更多的训练步骤，让差异更明显
        for _ in range(500):
            state = np.random.randint(n_states)
            action = np.random.randint(n_actions)
            # 使用更结构化的奖励：有些状态-动作对有更高的期望奖励
            if (state + action) % 7 == 0:  # 创造一些"好"的状态-动作对
                reward = np.random.normal(1.0, 0.5)  # 较高的期望奖励
            else:
                reward = np.random.normal(-0.2, 0.8)  # 较低的期望奖励
            next_state = np.random.randint(n_states)
            done = np.random.random() < 0.05  # 降低episode结束概率
            
            ql_agent.update(state, action, reward, next_state, done)
            dql_agent.update(state, action, reward, next_state, done)
        
        # 分析过估计
        ql_over = analyze_overestimation(ql_agent, n_states, n_actions)
        dql_over = analyze_overestimation(dql_agent, n_states, n_actions)
        
        qlearning_overestimations.append(ql_over)
        doubleq_overestimations.append(dql_over)
    
    # Visualization results
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 1. Overestimation comparison
    axes[0, 0].bar(['Q-Learning', 'Double Q-Learning'], 
                   [np.mean(qlearning_overestimations), np.mean(doubleq_overestimations)],
                   yerr=[np.std(qlearning_overestimations), np.std(doubleq_overestimations)],
                   capsize=10, alpha=0.7)
    axes[0, 0].set_ylabel('Average Overestimation')
    axes[0, 0].set_title('Q-Learning vs Double Q-Learning Overestimation Comparison')
    axes[0, 0].grid(True, alpha=0.3)
    
    # 2. Distribution comparison
    axes[0, 1].boxplot([qlearning_overestimations, doubleq_overestimations], 
                       labels=['Q-Learning', 'Double Q-Learning'])
    axes[0, 1].set_ylabel('Overestimation')
    axes[0, 1].set_title('Overestimation Distribution')
    axes[0, 1].grid(True, alpha=0.3)
    
    # 3. Q-value distribution example
    axes[1, 0].hist(qlearning_overestimations, bins=15, alpha=0.7, label='Q-Learning', density=True)
    axes[1, 0].hist(doubleq_overestimations, bins=15, alpha=0.7, label='Double Q-Learning', density=True)
    axes[1, 0].set_xlabel('Overestimation')
    axes[1, 0].set_ylabel('Density')
    axes[1, 0].set_title('Overestimation Distribution Density')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    # 4. Algorithm principle diagram
    axes[1, 1].axis('off')
    explanation = """Double Q-Learning Principle:
    
    Standard Q-Learning:
    target = R + γ * max_a' Q(s', a')
    Problem: max operation uses the same noisy Q estimate
    
    Double Q-Learning:
    1. Use Q1 to select action: a* = argmax_a' Q1(s', a')
    2. Use Q2 to evaluate action: target = R + γ * Q2(s', a*)
    
    Effects:
    • Decouples "selection" and "evaluation"
    • Reduces overestimation bias
    • Improves stability"""
    
    axes[1, 1].text(0.1, 0.5, explanation, fontsize=10, 
                   verticalalignment='center', family='monospace')
    
    plt.tight_layout()
    plt.show()
    
    print(f"\n统计结果:")
    print(f"Q-Learning 平均过估计量: {np.mean(qlearning_overestimations):.3f} ± {np.std(qlearning_overestimations):.3f}")
    print(f"Double Q-Learning 平均过估计量: {np.mean(doubleq_overestimations):.3f} ± {np.std(doubleq_overestimations):.3f}")
    print(f"过估计减少比例: {(np.mean(qlearning_overestimations) - np.mean(doubleq_overestimations)) / np.mean(qlearning_overestimations) * 100:.1f}%")
    
    print(f"\n结论:")
    print("1. Q-Learning由于max操作，系统性地高估Q值")
    print("2. Double Q-Learning通过两个独立的Q函数，有效减少了过估计")
    print("3. 在深度强化学习中，这个思想被用于Double DQN")

def analyze_overestimation(agent, n_states, n_actions):
    """分析智能体的过估计程度"""
    overestimations = []
    
    for state in range(n_states):
        # 获取真实Q值
        true_q = agent.get_true_q_values(state)
        
        if len(true_q) > 0:
            # 计算估计的Q值（带噪声）
            estimated_q = [agent.get_q_value(state, a) for a in range(n_actions)]
            
            # 找到估计的最大值对应的动作
            best_estimated_action = np.argmax(estimated_q)
            
            # 计算过估计：估计的最大值 - 该动作的真实值
            overestimation = estimated_q[best_estimated_action] - true_q[best_estimated_action]
            overestimations.append(overestimation)
    
    return np.mean(overestimations) if overestimations else 0.0

def demonstrate_maximization_bias():
    """专门演示最大化偏差问题"""
    print("\n=== 最大化偏差演示 ===")
    print("展示Double Q-Learning如何解决max操作的过估计问题")
    
    # 简化的例子：10个动作，真实Q值都相同
    n_actions = 10
    true_q_value = 0.0  # 所有动作的真实Q值都是0
    noise_std = 1.0     # 估计噪声标准差
    
    # 模拟1000次估计
    n_trials = 1000
    standard_max_values = []
    double_q_values = []
    
    np.random.seed(42)
    
    for _ in range(n_trials):
        # 标准Q-Learning：用同一组带噪声的估计做max选择和评估
        noisy_estimates = np.random.normal(true_q_value, noise_std, n_actions)
        max_action = np.argmax(noisy_estimates)
        standard_max_value = noisy_estimates[max_action]  # 用同一个估计值
        standard_max_values.append(standard_max_value)
        
        # Double Q-Learning模拟：用一组估计选择，用另一组估计评估
        q1_estimates = np.random.normal(true_q_value, noise_std, n_actions)
        q2_estimates = np.random.normal(true_q_value, noise_std, n_actions)
        
        max_action_q1 = np.argmax(q1_estimates)  # 用Q1选择动作
        double_q_value = q2_estimates[max_action_q1]  # 用Q2评估该动作
        double_q_values.append(double_q_value)
    
    # 分析结果
    standard_mean = np.mean(standard_max_values)
    double_mean = np.mean(double_q_values)
    bias_reduction = standard_mean - double_mean
    
    print(f"\n真实Q值: {true_q_value:.3f}")
    print(f"标准max操作平均值: {standard_mean:.3f} (过估计: {standard_mean - true_q_value:.3f})")
    print(f"Double Q平均值: {double_mean:.3f} (过估计: {double_mean - true_q_value:.3f})")
    print(f"偏差减少量: {bias_reduction:.3f}")
    print(f"偏差减少比例: {bias_reduction/standard_mean*100:.1f}%")
    
    # 可视化
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # 分布对比
    ax1.hist(standard_max_values, bins=30, alpha=0.7, label='Standard Q-Learning (max)', density=True)
    ax1.hist(double_q_values, bins=30, alpha=0.7, label='Double Q-Learning', density=True)
    ax1.axvline(true_q_value, color='red', linestyle='--', linewidth=2, label=f'True Value ({true_q_value})')
    ax1.axvline(standard_mean, color='blue', linestyle='--', alpha=0.8, label=f'Standard Mean ({standard_mean:.2f})')
    ax1.axvline(double_mean, color='orange', linestyle='--', alpha=0.8, label=f'Double Q Mean ({double_mean:.2f})')
    ax1.set_xlabel('Estimated Q-value')
    ax1.set_ylabel('Density')
    ax1.set_title('Maximization Bias Comparison')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 偏差对比
    methods = ['Standard\nQ-Learning', 'Double\nQ-Learning']
    biases = [standard_mean - true_q_value, double_mean - true_q_value]
    colors = ['red', 'green']
    
    bars = ax2.bar(methods, biases, color=colors, alpha=0.7)
    ax2.axhline(0, color='black', linestyle='-', linewidth=1)
    ax2.set_ylabel('Overestimation Bias')
    ax2.set_title('Bias Reduction by Double Q-Learning')
    ax2.grid(True, alpha=0.3)
    
    # 添加数值标签
    for bar, bias in zip(bars, biases):
        ax2.text(bar.get_x() + bar.get_width()/2, bias + 0.01, f'{bias:.3f}', 
                ha='center', va='bottom' if bias > 0 else 'top')
    
    plt.tight_layout()
    plt.show()
    
    print(f"\nDouble Q-Learning优化的核心问题:")
    print("1. 最大化偏差：max(Q(s,a)) 系统性地高估真实价值")
    print("2. 选择-评估耦合：同一个估计既用来选择又用来评估")
    print("3. Double Q通过独立的选择器和评估器解决这个问题")

# 运行两个演示
demonstrate_maximization_bias()
print("\n" + "="*60)
compare_overestimation()