27 - 数学到机器学习的过渡¶

import numpy as np

# 示例：计算矩阵的特征值和特征向量
A = np.array([[4, 2],
              [1, 3]])

# 计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(A)

print("矩阵 A:")
print(A)
print("\n特征值:")
print(eigenvalues)
print("\n特征向量（列）:")
print(eigenvectors)

import numpy as np

def pca_from_scratch(X, n_components=2):
    """
    从特征值分解实现PCA

    参数:
        X: 数据矩阵 (n_samples, n_features)
        n_components: 主成分数量

    返回:
        X_pca: 降维后的数据
        components: 主成分（特征向量）
        explained_variance: 解释的方差
    """
    # 1. 中心化数据
    X_centered = X - X.mean(axis=0)

    # 2. 计算协方差矩阵
    cov_matrix = np.cov(X_centered, rowvar=False)

    # 3. 特征值分解
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

    # 4. 选择最大的k个特征值对应的特征向量
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # 5. 投影数据
    components = eigenvectors[:, :n_components]
    X_pca = X_centered @ components

    return X_pca, components, eigenvalues[:n_components]

# 示例
np.random.seed(42)
# 生成相关数据
n_samples = 100
X = np.random.randn(n_samples, 5) @ np.random.randn(5, 3) + np.random.randn(1, 3)
X += np.array([10, 5, 3])  # 添加偏移

X_pca, components, ev = pca_from_scratch(X, n_components=2)
print(f"原始形状: {X.shape}")
print(f"降维后形状: {X_pca.shape}")
print(f"解释的方差比: {ev / ev.sum()}")

import numpy as np

# SVD示例
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

U, S, Vt = np.linalg.svd(A)

print("矩阵 A:")
print(A)
print("\n左奇异向量 U:")
print(U)
print("\n奇异值 S:")
print(S)
print("\n右奇异向量 V^T:")
print(Vt)

# 验证分解
print("\n重构矩阵 U @ Σ @ V^T:")
print(U @ np.diag(S) @ Vt)

def pca_via_svd(X, n_components=2):
    """
    使用SVD实现PCA

    对于中心化数据，PCA等价于截断SVD
    """
    # 中心化
    X_centered = X - X.mean(axis=0)

    # SVD（使用 econoded SVD 减少计算）
    U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)

    # 投影到主成分空间
    X_pca = U[:, :n_components] * S[:n_components]

    return X_pca, Vt[:n_components].T, S[:n_components]**2

# 比较两种方法
X_pca_eig, _, _ = pca_from_scratch(X, 2)
X_pca_svd, _, _ = pca_via_svd(X, 2)

print("特征值分解PCA结果:\n", X_pca_eig[:5])
print("\nSVD-PCA结果:\n", X_pca_svd[:5])

def truncated_svd_demo():
    """演示截断SVD降维"""
    # 生成一个秩为2的矩阵
    np.random.seed(123)
    u = np.random.randn(100, 2)
    v = np.random.randn(150, 2)
    A = u @ v.T

    print(f"原始矩阵形状: {A.shape}")
    print(f"原始矩阵秩: {np.linalg.matrix_rank(A)}")

    # 截断SVD
    k = 2
    U, S, Vt = np.linalg.svd(A, full_matrices=False)
    A_approx = U[:, :k] @ np.diag(S[:k]) @ Vt[:k]

    # 计算重构误差
    error = np.linalg.norm(A - A_approx, 'fro') / np.linalg.norm(A, 'fro')
    print(f"保留{k}个奇异值的重构误差: {error:.4f}")

    # 解释的方差
    total_var = S**2
    explained_ratio = total_var[:k].sum() / total_var.sum()
    print(f"解释的方差比例: {explained_ratio:.4f}")

truncated_svd_demo()

import numpy as np

class MatrixFactorization:
    """矩阵分解推荐系统"""

    def __init__(self, n_factors=10, learning_rate=0.01, 
                 regularization=0.1, n_epochs=100):
        self.n_factors = n_factors
        self.lr = learning_rate
        self.reg = regularization
        self.n_epochs = n_epochs

    def fit(self, R, mask):
        """
        训练矩阵分解模型

        参数:
            R: 评分矩阵 (n_users, n_items)
            mask: 布尔矩阵，标识观测到的评分位置
        """
        np.random.seed(42)
        n_users, n_items = R.shape

        # 初始化隐因子矩阵
        self.U = np.random.randn(n_users, self.n_factors) * 0.1
        self.V = np.random.randn(n_items, self.n_factors) * 0.1

        # 获取观测到的评分位置
        user_indices, item_indices = np.where(mask)
        n_ratings = len(user_indices)

        for epoch in range(self.n_epochs):
            # 随机梯度下降
            indices = np.random.permutation(n_ratings)
            total_loss = 0

            for idx in indices:
                i = user_indices[idx]
                j = item_indices[idx]

                # 预测评分
                pred = self.U[i] @ self.V[j]

                # 计算误差
                error = R[i, j] - pred
                total_loss += error ** 2

                # 更新隐因子
                self.U[i] += self.lr * (error * self.V[j] - self.reg * self.U[i])
                self.V[j] += self.lr * (error * self.U[i] - self.reg * self.V[j])

            if (epoch + 1) % 20 == 0:
                rmse = np.sqrt(total_loss / n_ratings)
                print(f"Epoch {epoch+1}/{self.n_epochs}, RMSE: {rmse:.4f}")

        return self

    def predict(self):
        """返回完整的预测评分矩阵"""
        return self.U @ self.V.T

# 示例
np.random.seed(42)
n_users, n_items, n_factors = 50, 100, 5

# 生成真实隐因子
true_U = np.random.randn(n_users, n_factors)
true_V = np.random.randn(n_items, n_factors)

# 生成评分矩阵（加入噪声）
R_true = true_U @ true_V.T
noise = np.random.randn(n_users, n_items) * 0.5
R = R_true + noise

# 创建观测掩码（70%观测到）
mask = np.random.rand(n_users, n_items) < 0.7

# 训练模型
model = MatrixFactorization(n_factors=5, n_epochs=100)
model.fit(R, mask)

# 预测
R_pred = model.predict()

# 评估
mask_test = ~mask
rmse = np.sqrt(((R[mask_test] - R_pred[mask_test])**2).mean())
print(f"测试集RMSE: {rmse:.4f}")

class ALS:
    """交替最小二乘矩阵分解"""

    def __init__(self, n_factors=10, regularization=0.1, n_epochs=20):
        self.n_factors = n_factors
        self.reg = regularization
        self.n_epochs = n_epochs

    def fit(self, R, mask):
        n_users, n_items = R.shape

        # 初始化
        np.random.seed(42)
        self.U = np.random.randn(n_users, self.n_factors) * 0.1
        self.V = np.random.randn(n_items, self.n_factors) * 0.1

        user_indices, item_indices = np.where(mask)

        for epoch in range(self.n_epochs):
            # 更新U（固定V）
            for i in range(n_users):
                items_i = np.where(mask[i, :])[0]
                if len(items_i) == 0:
                    continue
                V_i = self.V[items_i]
                R_i = R[i, items_i]
                self.U[i] = np.linalg.solve(
                    V_i.T @ V_i + self.reg * np.eye(self.n_factors),
                    V_i.T @ R_i
                )

            # 更新V（固定U）
            for j in range(n_items):
                users_j = np.where(mask[:, j])[0]
                if len(users_j) == 0:
                    continue
                U_j = self.U[users_j]
                R_j = R[users_j, j]
                self.V[j] = np.linalg.solve(
                    U_j.T @ U_j + self.reg * np.eye(self.n_factors),
                    U_j.T @ R_j
                )

            # 计算损失
            pred = self.U @ self.V.T
            loss = np.sqrt(((R[mask] - pred[mask])**2).mean())
            if (epoch + 1) % 5 == 0:
                print(f"Epoch {epoch+1}, RMSE: {loss:.4f}")

        return self

# 使用ALS训练
als_model = ALS(n_factors=5, n_epochs=20)
als_model.fit(R, mask)

import numpy as np
from scipy import stats

# 示例1：正态分布的参数估计
np.random.seed(42)
true_mean, true_std = 5.0, 2.0
samples = np.random.normal(true_mean, true_std, 100)

# MLE估计
mu_mle = samples.mean()
sigma_mle = np.sqrt(((samples - mu_mle)**2).sum() / len(samples))

print(f"真实参数: μ={true_mean}, σ={true_std}")
print(f"MLE估计: μ={mu_mle:.4f}, σ={sigma_mle:.4f}")

# 示例2：伯努利分布的参数估计（硬币问题）
# 10次抛硬币，7次正面
n_trials, n_heads = 10, 7
p_mle = n_heads / n_trials

# 解析解 vs 数值解
from scipy.optimize import minimize

def neg_log_likelihood(p, k, n):
    """负对数似然（需要优化）"""
    if p <= 0 or p >= 1:
        return np.inf
    return -stats.binom.logpmf(k, n, p)

result = minimize(neg_log_likelihood, 0.5, args=(n_heads, n_trials), 
                  bounds=[(0.001, 0.999)])
p_numeric = result.x[0]

print(f"\n硬币问题:")
print(f"解析MLE: {p_mle:.4f}")
print(f"数值MLE: {p_numeric:.4f}")

import numpy as np
from scipy import stats

# 示例：伯努利分布的贝叶斯估计
# 先验：Beta(α, β)
# 后验：Beta(α + k, β + n - k)

alpha_prior, beta_prior = 2, 2  # 弱信息先验
n_trials, n_heads = 10, 7

# 后验参数
alpha_post = alpha_prior + n_heads
beta_post = beta_prior + n_trials - n_heads

# 后验均值（作为点估计）
p_bayes = alpha_post / (alpha_post + beta_post)

print("贝叶斯估计示例（硬币问题）:")
print(f"先验: Beta({alpha_prior}, {beta_prior})")
print(f"观测: {n_trials}次试验, {n_heads}次正面")
print(f"后验: Beta({alpha_post}, {beta_post})")
print(f"后验均值: {p_bayes:.4f}")
print(f"MLE: {n_heads/n_trials:.4f}")

# 可视化
import matplotlib.pyplot as plt
x = np.linspace(0, 1, 100)
plt.plot(x, stats.beta.pdf(x, alpha_prior, beta_prior), 
         label='Prior: Beta(2,2)', linestyle='--')
plt.plot(x, stats.beta.pdf(x, alpha_post, beta_post), 
         label=f'Posterior: Beta({alpha_post},{beta_post})')
plt.axvline(p_bayes, color='red', linestyle=':', label=f'Posterior mean: {p_bayes:.3f}')
plt.xlabel('p (硬币正面概率)')
plt.ylabel('概率密度')
plt.legend()
plt.title('贝叶斯估计：从先验到后验')
plt.savefig('bayesian_estimation.png', dpi=100, bbox_inches='tight')
plt.show()

def mle_vs_bayesian_comparison():
    """MLE与贝叶斯估计对比：小样本场景"""
    np.random.seed(42)

    # 真实参数
    true_p = 0.7

    # 小样本情况
    small_sample = np.random.binomial(1, true_p, 5)  # 5次观测
    n_heads_small = small_sample.sum()
    n_trials_small = len(small_sample)

    print(f"小样本实验 (n={n_trials_small}, k={n_heads_small}):")

    # MLE
    p_mle_small = n_heads_small / n_trials_small

    # 贝叶斯（Beta先验）
    alpha_prior, beta_prior = 2, 2
    alpha_post = alpha_prior + n_heads_small
    beta_post = beta_prior + n_trials_small - n_heads_small
    p_bayes_small = alpha_post / (alpha_post + beta_post)

    print(f"  MLE: {p_mle_small:.4f}")
    print(f"  贝叶斯: {p_bayes_small:.4f}")
    print(f"  真实值: {true_p:.4f}")

    # 大样本情况
    large_sample = np.random.binomial(1, true_p, 1000)
    n_heads_large = large_sample.sum()
    n_trials_large = len(large_sample)

    p_mle_large = n_heads_large / n_trials_large
    alpha_post_l = alpha_prior + n_heads_large
    beta_post_l = beta_prior + n_trials_large - n_heads_large
    p_bayes_large = alpha_post_l / (alpha_post_l + beta_post_l)

    print(f"\n大样本实验 (n={n_trials_large}, k={n_heads_large}):")
    print(f"  MLE: {p_mle_large:.4f}")
    print(f"  贝叶斯: {p_bayes_large:.4f}")
    print(f"  真实值: {true_p:.4f}")

    print("\n结论：大样本时两者趋同，小样本时贝叶斯更稳定")

mle_vs_bayesian_comparison()

import numpy as np

class EM_GaussianMixture:
    """EM算法估计高斯混合模型参数"""

    def __init__(self, n_components=2, max_iter=100, tol=1e-6):
        self.n_components = n_components
        self.max_iter = max_iter
        self.tol = tol

    def fit(self, X):
        n_samples, n_features = X.shape

        # 初始化参数
        np.random.seed(42)

        # 随机选择初始均值
        idx = np.random.choice(n_samples, self.n_components, replace=False)
        self.means = X[idx].copy()

        # 初始化协方差（单位矩阵）
        self.covs = np.array([np.eye(n_features) for _ in range(self.n_components)])

        # 初始化混合系数
        self.weights = np.ones(self.n_components) / self.n_components

        # EM迭代
        for iteration in range(self.max_iter):
            # E步：计算每个样本属于每个成分的后验概率
            responsibilities = self._e_step(X)

            # M步：更新参数
            self._m_step(X, responsibilities)

            # 计算对数似然
            log_likelihood = self._compute_log_likelihood(X)

            if iteration > 0 and abs(log_likelihood - prev_ll) < self.tol:
                print(f"收敛于第 {iteration} 次迭代，对数似然: {log_likelihood:.4f}")
                break
            prev_ll = log_likelihood

        return self

    def _e_step(self, X):
        """E步：计算责任度"""
        n_samples = X.shape[0]
        responsibilities = np.zeros((n_samples, self.n_components))

        for k in range(self.n_components):
            responsibilities[:, k] = self.weights[k] * self._gaussian_pdf(
                X, self.means[k], self.covs[k]
            )

        # 归一化
        responsibilities /= responsibilities.sum(axis=1, keepdims=True) + 1e-10
        return responsibilities

    def _m_step(self, X, responsibilities):
        """M步：更新参数"""
        n_samples = X.shape[0]
        N_k = responsibilities.sum(axis=0)  # 每个成分的有效样本数

        # 更新混合系数
        self.weights = N_k / n_samples

        # 更新均值和协方差
        for k in range(self.n_components):
            # 均值
            self.means[k] = (responsibilities[:, k:k+1] * X).sum(axis=0) / N_k[k]

            # 协方差
            diff = X - self.means[k]
            self.covs[k] = (responsibilities[:, k:k+1] * diff).T @ diff / N_k[k]

    def _gaussian_pdf(self, X, mean, cov):
        """多元高斯分布的概率密度"""
        diff = X - mean
        cov_inv = np.linalg.inv(cov + 1e-6 * np.eye(cov.shape[0]))
        cov_det = np.linalg.det(cov + 1e-6 * np.eye(cov.shape[0]))

        exponent = -0.5 * np.sum(diff @ cov_inv * diff, axis=1)
        return np.exp(exponent) / np.sqrt((2 * np.pi) ** X.shape[1] * cov_det)

    def _compute_log_likelihood(self, X):
        """计算对数似然"""
        log_likelihood = np.zeros(X.shape[0])
        for k in range(self.n_components):
            log_likelihood += self.weights[k] * self._gaussian_pdf(
                X, self.means[k], self.covs[k]
            )
        return np.log(log_likelihood + 1e-10).sum()

# 示例：使用EM算法估计GMM
np.random.seed(42)

# 生成混合高斯数据
n1 = 200
n2 = 150
X1 = np.random.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], n1)
X2 = np.random.multivariate_normal([5, 5], [[1, -0.5], [-0.5, 1]], n2)
X = np.vstack([X1, X2])

# 使用EM算法拟合
em = EM_GaussianMixture(n_components=2, max_iter=100)
em.fit(X)

print("\n估计的参数:")
for k in range(2):
    print(f"\n成分 {k+1}:")
    print(f"  权重: {em.weights[k]:.4f}")
    print(f"  均值: {em.means[k]}")
    print(f"  协方差:\n{em.covs[k]}")

def em_convergence_demo():
    """演示EM算法的收敛性"""
    import matplotlib.pyplot as plt

    np.random.seed(42)

    # 生成双峰数据
    X1 = np.random.normal(0, 1, 300)
    X2 = np.random.normal(5, 1, 200)
    X = np.concatenate([X1, X2]).reshape(-1, 1)

    # 记录每次迭代的似然
    log_likelihoods = []

    class EM_with_tracking(EM_GaussianMixture):
        def _compute_log_likelihood(self, X):
            ll = super()._compute_log_likelihood(X)
            log_likelihoods.append(ll)
            return ll

    em = EM_with_tracking(n_components=2, max_iter=50)
    em.fit(X)

    # 可视化收敛过程
    plt.figure(figsize=(10, 4))

    plt.subplot(1, 2, 1)
    plt.hist(X.flatten(), bins=50, density=True, alpha=0.5, label='数据')
    x_range = np.linspace(-4, 10, 200).reshape(-1, 1)
    for k in range(2):
        pdf = em.weights[k] * stats.norm.pdf(x_range, em.means[k], 
                                                np.sqrt(em.covs[k, 0, 0]))
        plt.plot(x_range, pdf, label=f'成分 {k+1}')
    plt.xlabel('x')
    plt.ylabel('概率密度')
    plt.title('EM估计的GMM')
    plt.legend()

    plt.subplot(1, 2, 2)
    plt.plot(log_likelihoods, 'b-', linewidth=2)
    plt.xlabel('迭代次数')
    plt.ylabel('对数似然')
    plt.title('EM算法收敛曲线')
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('em_convergence.png', dpi=100, bbox_inches='tight')
    plt.show()

em_convergence_demo()

import numpy as np

# 示例：简单贝叶斯网络的联合分布
# 场景：天气 -> 是否出门 -> 是否带伞
# 节点: W(天气), O(出门), U(带伞)

# 条件概率表
# P(W=好) = 0.7, P(W=坏) = 0.3
p_weather_good = 0.7

# P(O=出门|W=好) = 0.9, P(O=出门|W=坏) = 0.2
p_outdoor_given_good = 0.9
p_outdoor_given_bad = 0.2

# P(U=带伞|O=出门) = 0.1, P(U=带伞|O=不出门) = 0.8
p_umbrella_given_out = 0.1
p_umbrella_given_home = 0.8

def compute_joint_probability(weather, outdoor, umbrella):
    """计算联合概率"""
    # P(W)
    p_w = p_weather_good if weather == 'good' else 1 - p_weather_good

    # P(O|W)
    if outdoor == 'out':
        p_o_given_w = p_outdoor_given_good if weather == 'good' else p_outdoor_given_bad
    else:
        p_o_given_w = (1 - p_outdoor_given_good) if weather == 'good' else (1 - p_outdoor_given_bad)

    # P(U|O)
    if outdoor == 'out':
        p_u_given_o = p_umbrella_given_out
    else:
        p_u_given_o = p_umbrella_given_home

    if umbrella == 'yes':
        p_u = p_u_given_o
    else:
        p_u = 1 - p_u_given_o

    return p_w * p_o_given_w * p_u

# 计算所有可能的联合概率
print("联合分布:")
print("天气  出门  带伞  概率")
for weather in ['good', 'bad']:
    for outdoor in ['out', 'home']:
        for umbrella in ['yes', 'no']:
            p = compute_joint_probability(weather, outdoor, umbrella)
            print(f"{weather:5s} {outdoor:5s} {umbrella:5s} {p:.4f}")

# 边缘概率
print("\n边缘概率:")
p_out = sum(compute_joint_probability(w, 'out', u) 
            for w in ['good', 'bad'] for u in ['yes', 'no'])
print(f"P(出门) = {p_out:.4f}")

# 条件概率推理
# P(带伞|天气好) = ?
p_umbrella_given_weather_good = (
    compute_joint_probability('good', 'out', 'yes') + 
    compute_joint_probability('good', 'home', 'yes')
) / p_weather_good
print(f"P(带伞|天气好) = {p_umbrella_given_weather_good:.4f}")

# 示例：简单的马尔可夫随机场 - 图像去噪
# 使用Ising模型进行二值图像去噪

import numpy as np

def ising_denoising(noisy_image, beta=1.0, n_iter=10):
    """
    使用Ising模型的图像去噪

    Ising模型: P(x) ∝ exp(β * Σ_{(i,j)相邻} x_i * x_j)
    """
    image = noisy_image.copy()
    n_rows, n_cols = image.shape

    # 定义邻居（上下左右）
    neighbors = [(-1, 0), (1, 0), (0, -1), (0, 1)]

    for iteration in range(n_iter):
        for i in range(n_rows):
            for j in range(n_cols):
                # 计算当前像素的局部场
                local_field = 0
                for di, dj in neighbors:
                    ni, nj = i + di, j + dj
                    if 0 <= ni < n_rows and 0 <= nj < n_cols:
                        local_field += image[ni, nj]

                # 吉布斯采样
                prob_plus = 1 / (1 + np.exp(-2 * beta * local_field))
                image[i, j] = 1 if np.random.rand() < prob_plus else -1

    return image

# 测试
np.random.seed(42)

# 原始图像（二值：+1 或 -1）
original = np.array([[1, 1, 1, -1, -1],
                     [1, 1, 1, -1, -1],
                     [1, 1, 1, -1, -1],
                     [-1, -1, -1, 1, 1],
                     [-1, -1, -1, 1, 1]])

# 加入噪声
noise_level = 0.2
noisy = original.copy()
for i in range(5):
    for j in range(5):
        if np.random.rand() < noise_level:
            noisy[i, j] *= -1

print("原始图像:")
print(original)
print("\n噪声图像:")
print(noisy)

# 去噪
denoised = ising_denoising(noisy, beta=1.0, n_iter=20)
print("\n去噪后图像:")
print(denoised)

# 计算准确率
accuracy = (denoised == original).mean()
print(f"\n去噪准确率: {accuracy:.2%}")

# 示例：因子图上的信念传播 - 简单链状结构
# X1 -- f1 -- X2 -- f2 -- X3

def belief_propagation_chain():
    """
    链状结构的信念传播

    消息从左向右传递，然后从右向左传递
    """
    # 定义变量
    n_states = 2  # 二值变量

    # 势函数（因子）
    # f1(x1, x2): 相邻变量的兼容性
    f1 = np.array([[0.9, 0.1],  # x1=0时
                   [0.1, 0.9]]) # x1=1时

    f2 = np.array([[0.9, 0.1],  # x2=0时
                   [0.1, 0.9]]) # x2=1时

    # 一元势函数（来自观测）
    phi1 = np.array([0.8, 0.2])  # P(x1) 的观测
    phi2 = np.array([0.3, 0.7])  # P(x2) 的观测
    phi3 = np.array([0.6, 0.4])  # P(x3) 的观测

    # 从左向右传递消息
    # m1_2(x2) = Σ_x1 [phi1(x1) * f1(x1,x2)]
    m1_2 = np.zeros(n_states)
    for x2 in range(n_states):
        for x1 in range(n_states):
            m1_2[x2] += phi1[x1] * f1[x1, x2]
    m1_2 /= m1_2.sum()  # 归一化

    # m2_3(x3) = Σ_x2 [m1_2(x2) * phi2(x2) * f2(x2,x3)]
    m2_3 = np.zeros(n_states)
    for x3 in range(n_states):
        for x2 in range(n_states):
            m2_3[x3] += m1_2[x2] * phi2[x2] * f2[x2, x3]
    m2_3 /= m2_3.sum()

    # 计算边缘分布
    # P(x1) ∝ phi1(x1)
    p_x1 = phi1 / phi1.sum()

    # P(x2) ∝ phi2(x2) * m1_2(x2)
    p_x2 = phi2 * m1_2
    p_x2 /= p_x2.sum()

    # P(x3) ∝ phi3(x3) * m2_3(x3)
    p_x3 = phi3 * m2_3
    p_x3 /= p_x3.sum()

    print("信念传播结果:")
    print(f"P(x1=0)={p_x1[0]:.3f}, P(x1=1)={p_x1[1]:.3f}")
    print(f"P(x2=0)={p_x2[0]:.3f}, P(x2=1)={p_x2[1]:.3f}")
    print(f"P(x3=0)={p_x3[0]:.3f}, P(x3=1)={p_x3[1]:.3f}")

    return p_x1, p_x2, p_x3

belief_propagation_chain()

import numpy as np
import matplotlib.pyplot as plt

def convex_examples():
    """凸函数与凹函数示例"""
    x = np.linspace(-3, 3, 100)

    # 凸函数示例
    f1 = x**2  # 二次函数（凸）
    f2 = np.exp(x)  # 指数函数（凸）
    f3 = np.abs(x)  # 绝对值（凸）

    # 非凸函数示例
    f4 = x**3  # 三次函数（非凸）
    f5 = np.sin(x)  # 正弦函数（非凸）

    fig, axes = plt.subplots(2, 3, figsize=(12, 8))

    functions = [
        (x**2, '凸函数: $f(x)=x^2$'),
        (np.exp(x), '凸函数: $f(x)=e^x$'),
        (np.abs(x), '凸函数: $f(x)=|x|$'),
        (x**3, '非凸: $f(x)=x^3$'),
        (np.sin(x), '非凸: $f(x)=\sin(x)$'),
        (x**4 - 2*x**2, '非凸: $f(x)=x^4-2x^2$')
    ]

    for ax, (f, title) in zip(axes.flat, functions):
        ax.plot(x, f, 'b-', linewidth=2)
        ax.axhline(y=0, color='k', linestyle='--', alpha=0.3)
        ax.set_title(title)
        ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('convex_functions.png', dpi=100, bbox_inches='tight')
    plt.show()

convex_examples()

def gradient_descent_convergence():
    """演示梯度下降在不同函数上的收敛"""

    def gradient_descent(f, grad_f, x0, lr=0.1, tol=1e-6, max_iter=1000):
        """梯度下降"""
        x = x0
        trajectory = [x]

        for _ in range(max_iter):
            gradient = grad_f(x)
            x_new = x - lr * gradient

            if abs(x_new - x) < tol:
                break

            x = x_new
            trajectory.append(x)

        return x, trajectory

    # 凸函数：x^2
    f_convex = lambda x: x**2
    grad_convex = lambda x: 2*x

    # 非凸函数：x^4 - 2x^2
    f_nonconvex = lambda x: x**4 - 2*x**2
    grad_nonconvex = lambda x: 4*x**3 - 4*x

    # 凸函数优化
    x_opt_c, traj_c = gradient_descent(f_convex, grad_convex, x0=5, lr=0.1)
    print(f"凸函数 x^2:")
    print(f"  最优解: x* = {x_opt_c:.6f}")
    print(f"  迭代次数: {len(traj_c)}")

    # 非凸函数优化（从不同起点）
    starts = [-2, -1, 0, 1, 2]
    print(f"\n非凸函数 x^4 - 2x^2:")
    for x0 in starts:
        x_opt, traj = gradient_descent(f_nonconvex, grad_nonconvex, x0=x0, lr=0.01)
        print(f"  起点 x0={x0}: 收敛到 x* = {x_opt:.6f}, 迭代 {len(traj)} 次")

gradient_descent_convergence()

import numpy as np

def lagrange_multiplier_example():
    """
    示例：在约束 x^2 + y^2 = 1 下最小化 f(x,y) = x + y

    几何解释：在单位圆上找到使 x+y 最大的点
    """
    from scipy.optimize import minimize

    # 定义目标函数和约束
    def objective(x):
        return x[0] + x[1]  # 最小化 x + y

    def constraint(x):
        return x[0]**2 + x[1]**2 - 1  # = 0

    # 使用SLSQP求解器
    from scipy.optimize import minimize

    # 方法1：使用等式约束
    result = minimize(objective, x0=[0.5, 0.5], 
                     constraints={'type': 'eq', 'fun': constraint})

    print("拉格朗日乘数法示例:")
    print(f"  最优解: x* = {result.x}")
    print(f"  目标函数值: {result.fun:.4f}")
    print(f"  约束满足度: {constraint(result.x):.6f}")

    # 解析验证
    # 使用拉格朗日乘数法：
    # ∇f = (1, 1), ∇h = (2x, 2y)
    # 1 + λ*2x = 0 => x = -1/(2λ)
    # 1 + λ*2y = 0 => y = -1/(2λ)
    # 代入约束: 2*(1/(2λ)^2) = 1 => λ = -1/√2 或 λ = 1/√2
    # 最小化 x+y，需要 λ = 1/√2, x=y=-1/√2
    print(f"\n  解析解: x* = y* = {-1/np.sqrt(2):.4f}")

lagrange_multiplier_example()

def kkt_conditions_example():
    """
    示例：理解KKT条件

    问题：min x^2 s.t. x >= 1
    即 min x^2 s.t. 1 - x <= 0
    """
    from scipy.optimize import minimize

    def objective(x):
        return x**2

    def inequality_constraint(x):
        return 1 - x  # g(x) = 1 - x <= 0

    # 求解
    result = minimize(objective, x0=0, 
                     constraints={'type': 'ineq', 'fun': lambda x: x - 1})

    x_opt = result.x[0]
    print("KKT条件示例: min x² s.t. x ≥ 1")
    print(f"\n最优解: x* = {x_opt:.4f}")

    # 验证KKT条件
    # 1. 原始可行性: x >= 1 => g(x) = 1 - x <= 0
    g = 1 - x_opt
    print(f"  g(x*) = {g:.4f} <= 0 ✓")

    # 2. 平稳性: 2x* - λ = 0 => λ = 2x*
    # 由于约束是激活的（g(x*) = 0），λ > 0
    # 从求解器获取λ
    print(f"  约束边界处，最优解位于约束边界 x=1")

    # 解析解
    print(f"\n解析解:")
    print(f"  当 x* >= 1 时，最小化 x²，x* = 1")
    print(f"  此时 λ = 2，满足 λ >= 0 和 λ * g(x*) = 0")

kkt_conditions_example()

def svm_dual_problem():
    """
    SVM的KKT条件解释

    原始问题：
    min_{w,b} ½||w||² s.t. y_i(w·x_i + b) >= 1

    对偶问题：
    max_α Σα_i - ½ Σ_{i,j} α_i α_j y_i y_j x_i·x_j
    s.t. α_i >= 0, Σα_i y_i = 0
    """
    from sklearn.svm import SVC
    import numpy as np

    # 生成线性可分数据
    np.random.seed(42)
    X1 = np.random.randn(50, 2) + np.array([2, 2])
    X2 = np.random.randn(50, 2) + np.array([-2, -2])
    X = np.vstack([X1, X2])
    y = np.array([1]*50 + [-1]*50)

    # 训练SVM
    svm = SVC(kernel='linear', C=1.0)
    svm.fit(X, y)

    print("SVM中的KKT条件:")
    print(f"  支持向量数量: {sum(svm.n_support_)}")
    print(f"  权重 w: {svm.coef_[0]}")
    print(f"  偏置 b: {svm.intercept_[0]:.4f}")

    # KKT条件验证
    # 对于支持向量：y_i(w·x_i + b) = 1，α_i > 0
    # 对于非支持向量：y_i(w·x_i + b) > 1，α_i = 0
    decision_values = svm.decision_function(X)

    print(f"\nKKT条件验证:")
    support_idx = svm.support_
    non_support_idx = np.setdiff1d(np.arange(len(X)), support_idx)

    print(f"  支持向量处 y·f(x) ≈ 1:")
    for i in support_idx[:3]:
        print(f"    样本{i}: y={y[i]}, f(x)={decision_values[i]:.4f}, y·f(x)={y[i]*decision_values[i]:.4f}")

    print(f"\n  非支持向量处 y·f(x) > 1:")
    for i in non_support_idx[:3]:
        print(f"    样本{i}: y={y[i]}, f(x)={decision_values[i]:.4f}, y·f(x)={y[i]*decision_values[i]:.4f}")

svm_dual_problem()

def primal_dual_example():
    """
    示例：原始-对偶问题

    原始问题：
    min x^2 + y^2
    s.t. x + y = 1

    对偶问题：
    max_λ -¼λ² + λ
    """
    import numpy as np
    from scipy.optimize import minimize

    # 原始问题求解
    def primal_objective(x):
        return x[0]**2 + x[1]**2

    def equality_constraint(x):
        return x[0] + x[1] - 1

    result_primal = minimize(primal_objective, x0=[0, 0],
                            constraints={'type': 'eq', 'fun': equality_constraint})

    print("原始问题: min x² + y² s.t. x + y = 1")
    print(f"  最优解: x* = {result_primal.x}")
    print(f"  最优值: f(x*) = {result_primal.fun:.4f}")

    # 对偶函数
    def lagrangian(x, lam):
        return x[0]**2 + x[1]**2 + lam * (x[0] + x[1] - 1)

    # 对偶问题求解（固定λ，求最小化L，然后最大化）
    def dual_function(lam):
        # 对给定的λ，最小化拉格朗日函数
        result = minimize(lambda x: lagrangian(x, lam), x0=[0, 0])
        return -result.fun  # 对偶函数是拉格朗日的下界

    # 网格搜索找最优λ
    lambdas = np.linspace(-2, 4, 100)
    dual_values = [dual_function(lam) for lam in lambdas]
    lam_opt = lambdas[np.argmax(dual_values)]

    print(f"\n对偶问题:")
    print(f"  最优λ: {lam_opt:.4f}")
    print(f"  对偶最优值: {max(dual_values):.4f}")
    print(f"\n强对偶性验证: {result_primal.fun:.4f} ≈ {max(dual_values):.4f}")

primal_dual_example()