00 - 数学基础¶

import numpy as np
import matplotlib.pyplot as plt

# 可视化线性变换
A = np.array([[2, 1], [1, 2]])  # 对称矩阵
eigenvalues, eigenvectors = np.linalg.eig(A)  # np.linalg线性代数运算

print(f"矩阵 A:\n{A}")
print(f"\n特征值: {eigenvalues}")
print(f"特征向量:\n{eigenvectors}")

# 几何解释：特征向量是变换中方向不变的向量
# 特征值表示在该方向上的缩放因子

# PCA的数学本质：特征值分解
def pca_eig(X, n_components=2):
    """基于特征值分解的PCA"""
    # 中心化
    X_centered = X - np.mean(X, axis=0)

    # 协方差矩阵
    cov_matrix = np.cov(X_centered.T)

    # 特征值分解（协方差矩阵是对称矩阵，用eigh更稳定且保证实数结果）
    eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

    # 按特征值排序
    idx = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # 投影
    components = eigenvectors[:, :n_components]
    X_transformed = X_centered @ components

    return X_transformed, eigenvalues, eigenvectors

# SVD低秩近似
def low_rank_approximation(A, k):
    """SVD低秩近似"""
    U, s, Vt = np.linalg.svd(A, full_matrices=False)

    # 保留前k个奇异值
    U_k = U[:, :k]
    s_k = s[:k]
    Vt_k = Vt[:k, :]

    # 重构
    A_approx = U_k @ np.diag(s_k) @ Vt_k

    # 计算近似误差
    error = np.linalg.norm(A - A_approx, 'fro') / np.linalg.norm(A, 'fro')

    return A_approx, error, s

# 示例：图像压缩
from sklearn.datasets import load_sample_image
china = load_sample_image("china.jpg")
gray = np.mean(china, axis=2)

A_approx, error, singular_values = low_rank_approximation(gray, k=50)
print(f"压缩率: {50 * (gray.shape[0] + gray.shape[1]) / (gray.shape[0] * gray.shape[1]):.2%}")
print(f"近似误差: {error:.4f}")

# 矩阵求导在机器学习中的应用：线性回归的闭式解

def linear_regression_closed_form(X, y):
    """
    最小二乘法的闭式解
    目标：min ||Xw - y||^2
    求导：2X^T(Xw - y) = 0
    解得：w = (X^T X)^(-1) X^T y
    """
    return np.linalg.inv(X.T @ X) @ X.T @ y

# 梯度推导验证
def verify_gradient():
    """验证梯度计算"""
    X = np.random.randn(100, 5)
    w_true = np.random.randn(5)
    y = X @ w_true + 0.1 * np.random.randn(100)

    # 数值梯度
    def loss(w):
        return np.sum((X @ w - y) ** 2)

    def numerical_gradient(w, eps=1e-5):
        grad = np.zeros_like(w)
        for i in range(len(w)):
            w_plus = w.copy()
            w_plus[i] += eps
            w_minus = w.copy()
            w_minus[i] -= eps
            grad[i] = (loss(w_plus) - loss(w_minus)) / (2 * eps)
        return grad

    w = np.random.randn(5)
    numerical_grad = numerical_gradient(w)
    analytical_grad = 2 * X.T @ (X @ w - y)

    print(f"数值梯度: {numerical_grad}")
    print(f"解析梯度: {analytical_grad}")
    print(f"误差: {np.linalg.norm(numerical_grad - analytical_grad)}")

verify_gradient()

# Hessian矩阵在优化中的应用

def compute_hessian(f, x, eps=1e-5):
    """数值计算Hessian矩阵"""
    n = len(x)
    hessian = np.zeros((n, n))

    for i in range(n):
        for j in range(n):
            x_pp = x.copy()
            x_pp[i] += eps
            x_pp[j] += eps

            x_pm = x.copy()
            x_pm[i] += eps
            x_pm[j] -= eps

            x_mp = x.copy()
            x_mp[i] -= eps
            x_mp[j] += eps

            x_mm = x.copy()
            x_mm[i] -= eps
            x_mm[j] -= eps

            hessian[i, j] = (f(x_pp) - f(x_pm) - f(x_mp) + f(x_mm)) / (4 * eps**2)

    return hessian

# 示例：二次函数的Hessian
f = lambda x: x[0]**2 + 2*x[1]**2 + 3*x[0]*x[1]
x = np.array([1.0, 1.0])
hessian = compute_hessian(f, x)
print(f"Hessian矩阵:\n{hessian}")
print(f"特征值: {np.linalg.eigvals(hessian)}")

import numpy as np

# 张量基础操作
A = np.random.randn(3, 4, 5)  # 3阶张量
B = np.random.randn(4, 5, 6)  # 3阶张量

# 张量缩并（在维度1和2上）
C = np.tensordot(A, B, axes=([1, 2], [0, 1]))
print(f"A形状: {A.shape}, B形状: {B.shape}")
print(f"缩并后C形状: {C.shape}")

# Einstein求和约定
def batch_matrix_multiply(A, B):
    """批量矩阵乘法: (batch, m, n) @ (batch, n, p) -> (batch, m, p)"""
    return np.einsum('bmn,bnp->bmp', A, B)

A_batch = np.random.randn(10, 3, 4)
B_batch = np.random.randn(10, 4, 5)
C_batch = batch_matrix_multiply(A_batch, B_batch)
print(f"批量矩阵乘法结果形状: {C_batch.shape}")

from tensorly.decomposition import parafac, tucker
import tensorly as tl

# CP分解
X = np.random.randn(10, 20, 30)
factors = parafac(X, rank=5)
print("CP分解因子形状:")
for i, factor in enumerate(factors.factors):  # enumerate同时获取索引和元素
    print(f"  因子{i}: {factor.shape}")

# Tucker分解
core, factors = tucker(X, ranks=[5, 10, 15])
print(f"\nTucker分解核心张量形状: {core.shape}")
for i, factor in enumerate(factors):
    print(f"  因子矩阵{i}: {factor.shape}")

def matrix_completion_svt(M, mask, rank=5, max_iter=100, tau=1.0):
    """
    奇异值阈值算法进行矩阵补全

    M: 观测矩阵（含缺失值设为0）
    mask: 观测掩码（1表示观测到，0表示缺失）
    """
    X = M.copy()

    for iteration in range(max_iter):
        # SVD
        U, s, Vt = np.linalg.svd(X, full_matrices=False)

        # 软阈值（Soft Thresholding）
        s_threshold = np.maximum(s - tau, 0)

        # 重构
        X_new = U @ np.diag(s_threshold) @ Vt

        # 保持观测值不变
        X = mask * M + (1 - mask) * X_new

        # 检查收敛
        if iteration % 10 == 0:
            error = np.linalg.norm(mask * (X - M), 'fro')
            print(f"迭代 {iteration}: 误差 = {error:.6f}")

    return X

# 示例
M_true = np.random.randn(50, 40) @ np.random.randn(40, 60)  # 秩40矩阵
mask = np.random.rand(50, 60) > 0.7  # 70%缺失
M_observed = M_true * mask

M_completed = matrix_completion_svt(M_observed, mask, rank=10)
error = np.linalg.norm(M_completed - M_true, 'fro') / np.linalg.norm(M_true, 'fro')
print(f"\n补全相对误差: {error:.4f}")

def robust_pca(M, lambda_=0.01, max_iter=100):
    """
    鲁棒PCA（Principal Component Pursuit）
    """
    L = np.zeros_like(M)
    S = np.zeros_like(M)
    Y = np.zeros_like(M)
    mu = 1.25 / np.linalg.norm(M, 2)

    for iteration in range(max_iter):
        # 更新L（低秩部分）
        U, s, Vt = np.linalg.svd(M - S + Y/mu, full_matrices=False)
        s_threshold = np.maximum(s - 1/mu, 0)
        L = U @ np.diag(s_threshold) @ Vt

        # 更新S（稀疏部分）
        temp = M - L + Y/mu
        S = np.sign(temp) * np.maximum(np.abs(temp) - lambda_/mu, 0)

        # 更新对偶变量
        Y = Y + mu * (M - L - S)

        # 增大mu
        mu = min(mu * 1.5, 1e6)

        if iteration % 10 == 0:
            error = np.linalg.norm(M - L - S, 'fro') / np.linalg.norm(M, 'fro')
            print(f"迭代 {iteration}: 相对误差 = {error:.6f}")

    return L, S

# 示例：去除图像中的噪声
M_clean = np.random.randn(100, 100) @ np.random.randn(100, 100)
noise = np.random.randn(100, 100) * 0.1
# 添加稀疏异常值
outliers = np.random.rand(100, 100) < 0.05
noise[outliers] = np.random.randn(np.sum(outliers)) * 5
M_noisy = M_clean + noise

L, S = robust_pca(M_noisy)
print(f"\n低秩部分秩: {np.linalg.matrix_rank(L, tol=1e-3)}")
print(f"稀疏部分非零元素比例: {np.mean(S != 0):.2%}")

import scipy.sparse as sp
from scipy.sparse.linalg import eigsh

def compute_graph_laplacian(adj_matrix, norm_type='sym'):
    """
    计算图拉普拉斯矩阵

    adj_matrix: 邻接矩阵（可以是稀疏矩阵）
    norm_type: 'sym'（对称归一化）或'rw'（随机游走）
    """
    # 度矩阵
    degrees = np.array(adj_matrix.sum(axis=1)).flatten()
    D = sp.diags(degrees)

    # 未归一化拉普拉斯
    L = D - adj_matrix

    if norm_type == 'sym':
        # 对称归一化
        D_inv_sqrt = sp.diags(1.0 / np.sqrt(degrees + 1e-10))
        L_norm = D_inv_sqrt @ L @ D_inv_sqrt
    elif norm_type == 'rw':
        # 随机游走归一化
        D_inv = sp.diags(1.0 / (degrees + 1e-10))
        L_norm = D_inv @ L
    else:
        L_norm = L

    return L_norm

# 示例：圆环图的拉普拉斯
n = 20
adj = sp.diags([np.ones(n-1), np.ones(n-1), np.ones(n), np.ones(n)],
               offsets=[-1, 1, 1-n, n-1], shape=(n, n))
L = compute_graph_laplacian(adj, norm_type='sym')

# 计算特征值
eigenvalues, eigenvectors = eigsh(L, k=5, which='SM')
print(f"最小的5个特征值: {eigenvalues}")

def graph_fourier_transform(signal, eigenvectors):
    """
    图傅里叶变换

    signal: 图信号 (N,)
    eigenvectors: 拉普拉斯矩阵的特征向量 (N, N)
    """
    return eigenvectors.T @ signal

def inverse_graph_fourier_transform(fourier_coeffs, eigenvectors):
    """
    逆图傅里叶变换
    """
    return eigenvectors @ fourier_coeffs

# 示例：在图上平滑信号
n = 100
# 创建路径图
adj = sp.diags([np.ones(n-1), np.ones(n-1)], offsets=[-1, 1])
L = compute_graph_laplacian(adj)

# 计算特征向量
eigenvalues, eigenvectors = eigsh(L, k=n-1, which='SM')
eigenvectors = np.column_stack([np.ones(n) / np.sqrt(n), eigenvectors])
eigenvalues = np.concatenate([[0], eigenvalues])

# 创建带噪声的信号
signal_clean = np.sin(2 * np.pi * np.arange(n) / n * 3)
signal_noisy = signal_clean + 0.5 * np.random.randn(n)

# 图傅里叶变换
fourier_coeffs = graph_fourier_transform(signal_noisy, eigenvectors)

# 低通滤波：保留低频成分
k = 10  # 保留前k个频率
fourier_filtered = fourier_coeffs.copy()
fourier_filtered[k:] = 0

# 逆变换
signal_filtered = inverse_graph_fourier_transform(fourier_filtered, eigenvectors)

print(f"去噪前MSE: {np.mean((signal_noisy - signal_clean)**2):.4f}")
print(f"去噪后MSE: {np.mean((signal_filtered - signal_clean)**2):.4f}")

from scipy import stats
import numpy as np

# 常见分布
# 1. 正态分布
normal_dist = stats.norm(loc=0, scale=1)
print(f"正态分布期望: {normal_dist.mean()}")
print(f"正态分布方差: {normal_dist.var()}")

# 2. 二项分布
binomial_dist = stats.binom(n=10, p=0.3)
print(f"\n二项分布期望: {binomial_dist.mean()}")
print(f"二项分布方差: {binomial_dist.var()}")

# 3. 泊松分布
poisson_dist = stats.poisson(mu=3)
print(f"\n泊松分布期望: {poisson_dist.mean()}")
print(f"泊松分布方差: {poisson_dist.var()}")

# 大数定律验证
sample_sizes = [10, 100, 1000, 10000, 100000]
for n in sample_sizes:
    samples = np.random.randn(n)
    sample_mean = np.mean(samples)
    print(f"样本量 {n:6d}: 样本均值 = {sample_mean:.6f}")

# 贝叶斯定理应用：垃圾邮件分类
class NaiveBayesClassifier:
    """
    多项式朴素贝叶斯（用于文本分类）
    """
    def __init__(self, alpha=1.0):  # __init__构造方法，创建对象时自动调用
        self.alpha = alpha  # 拉普拉斯平滑参数
        self.class_priors = {}
        self.word_probs = {}
        self.vocab = set()

    def fit(self, X, y):
        """
        X: 文档列表，每个文档是词列表
        y: 类别标签
        """
        n_samples = len(X)
        classes = np.unique(y)

        # 构建词汇表
        for doc in X:
            self.vocab.update(doc)
        self.vocab = list(self.vocab)
        word_to_idx = {word: i for i, word in enumerate(self.vocab)}

        for c in classes:
            # 计算先验概率 P(y=c)
            class_docs = [X[i] for i in range(n_samples) if y[i] == c]
            self.class_priors[c] = len(class_docs) / n_samples

            # 计算词频
            word_counts = np.zeros(len(self.vocab))
            for doc in class_docs:
                for word in doc:
                    word_counts[word_to_idx[word]] += 1

            # 计算条件概率 P(word|y=c)（带平滑）
            total_words = word_counts.sum()
            self.word_probs[c] = (word_counts + self.alpha) / (total_words + self.alpha * len(self.vocab))

        self.word_to_idx = word_to_idx
        return self

    def predict(self, X):
        predictions = []
        for doc in X:
            scores = {}
            for c in self.class_priors:
                # log P(y=c) + sum(log P(word|y=c))
                score = np.log(self.class_priors[c])
                for word in doc:
                    if word in self.word_to_idx:
                        score += np.log(self.word_probs[c][self.word_to_idx[word]])
                scores[c] = score
            predictions.append(max(scores, key=scores.get))
        return predictions

# 示例
docs = [
    ['free', 'win', 'money', 'now'],
    ['meeting', 'schedule', 'tomorrow'],
    ['win', 'prize', 'free'],
    ['project', 'deadline', 'meeting']
]
labels = ['spam', 'ham', 'spam', 'ham']

clf = NaiveBayesClassifier()
clf.fit(docs, labels)
test_doc = [['free', 'win']]
prediction = clf.predict(test_doc)
print(f"测试文档分类结果: {prediction[0]}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# 多元高斯分布可视化
def plot_gaussian_2d(mean, cov, title):
    x, y = np.mgrid[-3:3:.01, -3:3:.01]
    pos = np.dstack((x, y))
    rv = multivariate_normal(mean, cov)

    plt.figure(figsize=(8, 6))
    plt.contourf(x, y, rv.pdf(pos), levels=20, cmap='viridis')
    plt.colorbar(label='Probability Density')
    plt.title(title)
    plt.xlabel('x1')
    plt.ylabel('x2')

    # 绘制特征向量
    eigenvalues, eigenvectors = np.linalg.eig(cov)
    for i in range(2):
        vec = eigenvectors[:, i] * np.sqrt(eigenvalues[i])
        plt.arrow(mean[0], mean[1], vec[0], vec[1],
                 head_width=0.1, head_length=0.1, fc='red', ec='red')

    plt.show()

# 不同协方差矩阵的形状
plot_gaussian_2d([0, 0], [[1, 0], [0, 1]], '各向同性（球形）')
plot_gaussian_2d([0, 0], [[2, 0], [0, 0.5]], '轴对齐椭圆')
plot_gaussian_2d([0, 0], [[1, 0.8], [0.8, 1]], '旋转椭圆（正相关）')
plot_gaussian_2d([0, 0], [[1, -0.8], [-0.8, 1]], '旋转椭圆（负相关）')

def conditional_gaussian(mu, Sigma, x2_val, idx1, idx2):
    """
    计算多元高斯的条件分布

    mu: 均值向量
    Sigma: 协方差矩阵
    x2_val: 观测到的x2值
    idx1: 需要推断的变量索引
    idx2: 观测变量的索引
    """
    mu1 = mu[idx1]
    mu2 = mu[idx2]
    Sigma11 = Sigma[np.ix_(idx1, idx1)]
    Sigma12 = Sigma[np.ix_(idx1, idx2)]
    Sigma21 = Sigma[np.ix_(idx2, idx1)]
    Sigma22 = Sigma[np.ix_(idx2, idx2)]

    # 条件均值
    mu_cond = mu1 + Sigma12 @ np.linalg.inv(Sigma22) @ (x2_val - mu2)

    # 条件协方差
    Sigma_cond = Sigma11 - Sigma12 @ np.linalg.inv(Sigma22) @ Sigma21

    return mu_cond, Sigma_cond

# 示例
mu = np.array([0, 0])
Sigma = np.array([[1, 0.8], [0.8, 1]])

# 已知 x2 = 1.5，求 x1 的条件分布
x2_val = np.array([1.5])
mu_cond, Sigma_cond = conditional_gaussian(mu, Sigma, x2_val, [0], [1])

print(f"条件分布: x1 | x2=1.5 ~ N({mu_cond[0]:.3f}, {Sigma_cond[0,0]:.3f})")
print(f"先验分布: x1 ~ N({mu[0]:.3f}, {Sigma[0,0]:.3f})")
print(f"观测x2后，x1的不确定性从{Sigma[0,0]:.3f}降低到{Sigma_cond[0,0]:.3f}")

# 指数族分布的统一视角
import numpy as np
from scipy.special import logsumexp

class ExponentialFamily:
    """
    指数族分布的基类
    """
    def __init__(self, natural_params):
        self.eta = natural_params

    def log_partition(self):
        """对数配分函数 A(eta)"""
        raise NotImplementedError

    def sufficient_stats(self, x):
        """充分统计量 T(x)"""
        raise NotImplementedError

    def log_pdf(self, x):
        """对数概率密度"""
        return np.dot(self.eta, self.sufficient_stats(x)) - self.log_partition()  # np.dot计算点积/矩阵乘法

    def mean(self):
        """E[T(X)] = dA/deta"""
        raise NotImplementedError

    def variance(self):
        """Var(T(X)) = d^2A/deta^2"""
        raise NotImplementedError

class GaussianEF(ExponentialFamily):
    """正态分布的指数族形式"""
    def __init__(self, mu, sigma2):
        # 自然参数: [mu/sigma^2, -1/(2*sigma^2)]
        eta = np.array([mu / sigma2, -0.5 / sigma2])
        super().__init__(eta)  # super()调用父类方法
        self.mu = mu
        self.sigma2 = sigma2

    def log_partition(self):
        eta1, eta2 = self.eta
        return -eta1**2 / (4*eta2) - 0.5 * np.log(-2*eta2)

    def sufficient_stats(self, x):
        return np.array([x, x**2])

    def mean(self):
        return np.array([self.mu, self.sigma2 + self.mu**2])

    def variance(self):
        return np.array([[self.sigma2, 2*self.mu*self.sigma2],
                        [2*self.mu*self.sigma2, 2*self.sigma2**2 + 4*self.mu**2*self.sigma2]])

# 验证
gaussian = GaussianEF(mu=2, sigma2=4)
print(f"自然参数: {gaussian.eta}")
print(f"对数配分函数: {gaussian.log_partition():.4f}")
print(f"均值: {gaussian.mean()}")

import numpy as np

class MarkovChain:
    def __init__(self, transition_matrix):
        self.P = transition_matrix
        self.n_states = transition_matrix.shape[0]

    def simulate(self, initial_state, n_steps):
        """模拟马尔可夫链"""
        states = [initial_state]
        current = initial_state

        for _ in range(n_steps - 1):
            current = np.random.choice(self.n_states, p=self.P[current])
            states.append(current)

        return states

    def stationary_distribution(self):
        """计算平稳分布"""
        # 解方程: pi P = pi, sum(pi) = 1
        A = np.vstack([self.P.T - np.eye(self.n_states), np.ones(self.n_states)])
        b = np.zeros(self.n_states + 1)
        b[-1] = 1  # [-1]负索引取最后一个元素

        pi = np.linalg.lstsq(A, b, rcond=None)[0]
        return pi

    def n_step_transition(self, n):
        """n步转移矩阵"""
        return np.linalg.matrix_power(self.P, n)

# 示例：天气模型
# 状态: 0=晴天, 1=多云, 2=雨天
P = np.array([
    [0.7, 0.2, 0.1],  # 晴天 -> 晴天/多云/雨天
    [0.3, 0.4, 0.3],  # 多云 -> 晴天/多云/雨天
    [0.2, 0.3, 0.5]   # 雨天 -> 晴天/多云/雨天
])

mc = MarkovChain(P)
pi = mc.stationary_distribution()
print(f"平稳分布: 晴天={pi[0]:.3f}, 多云={pi[1]:.3f}, 雨天={pi[2]:.3f}")

# 长期预测
P_30 = mc.n_step_transition(30)
print(f"\n30天后（无论从什么天气开始）:")
print(f"晴天概率: {P_30[0, 0]:.3f}")
print(f"多云概率: {P_30[0, 1]:.3f}")
print(f"雨天概率: {P_30[0, 2]:.3f}")

import numpy as np
from scipy.spatial.distance import cdist

class GaussianProcess:
    """
    高斯过程回归
    """
    def __init__(self, kernel='rbf', noise=1e-5, **kernel_params):
        self.kernel_name = kernel
        self.noise = noise
        self.kernel_params = kernel_params
        self.X_train = None
        self.y_train = None
        self.K = None

    def kernel(self, X1, X2):
        """计算核矩阵"""
        if self.kernel_name == 'rbf':
            length_scale = self.kernel_params.get('length_scale', 1.0)
            sigma_f = self.kernel_params.get('sigma_f', 1.0)

            dist = cdist(X1, X2, 'sqeuclidean')
            return sigma_f**2 * np.exp(-0.5 * dist / length_scale**2)

        elif self.kernel_name == 'linear':
            return X1 @ X2.T

        else:
            raise ValueError(f"Unknown kernel: {self.kernel_name}")

    def fit(self, X, y):
        """训练高斯过程"""
        self.X_train = X
        self.y_train = y

        # 计算核矩阵
        K = self.kernel(X, X)
        self.K = K + self.noise * np.eye(len(X))
        self.L = np.linalg.cholesky(self.K)
        self.alpha = np.linalg.solve(self.L.T, np.linalg.solve(self.L, y))

        return self

    def predict(self, X_test, return_std=True):
        """预测"""
        # 计算测试点与训练点的核
        K_s = self.kernel(self.X_train, X_test)

        # 均值
        mu = K_s.T @ self.alpha

        if return_std:
            # 方差
            K_ss = self.kernel(X_test, X_test)
            v = np.linalg.solve(self.L, K_s)
            var = np.diag(K_ss) - np.sum(v**2, axis=0)
            std = np.sqrt(np.maximum(var, 0))
            return mu, std

        return mu

# 示例：用GP拟合函数
np.random.seed(42)

# 真实函数
def f(x):
    return np.sin(x) + 0.1 * x**2

# 训练数据
X_train = np.array([-3, -2, -1, 0, 1, 2, 3]).reshape(-1, 1)
y_train = f(X_train).ravel() + 0.1 * np.random.randn(len(X_train))

# 测试点
X_test = np.linspace(-4, 4, 100).reshape(-1, 1)

# 拟合GP
gp = GaussianProcess(kernel='rbf', length_scale=1.0, sigma_f=1.0)
gp.fit(X_train, y_train)
mu, std = gp.predict(X_test)

# 可视化
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.plot(X_test, f(X_test), 'r--', label='True function')
plt.plot(X_test, mu, 'b-', label='GP mean')
plt.fill_between(X_test.ravel(), mu - 2*std, mu + 2*std, alpha=0.2, label='95% confidence')
plt.scatter(X_train, y_train, c='red', s=100, zorder=5, label='Training data')
plt.legend()
plt.title('Gaussian Process Regression')
plt.show()

from scipy.optimize import minimize
import numpy as np
from scipy.stats import norm

def maximum_likelihood_estimation(data, distribution='normal'):
    """
    最大似然估计

    data: 观测数据
    distribution: 分布类型
    """
    if distribution == 'normal':
        # 正态分布的MLE有解析解
        mu_mle = np.mean(data)
        sigma2_mle = np.var(data, ddof=0)  # 注意：MLE用n而不是n-1
        return {'mu': mu_mle, 'sigma2': sigma2_mle}

    elif distribution == 'exponential':
        # 指数分布: lambda_MLE = 1/mean
        lambda_mle = 1 / np.mean(data)
        return {'lambda': lambda_mle}

    else:
        raise ValueError(f"Unknown distribution: {distribution}")

# 示例：估计正态分布参数
true_mu, true_sigma = 5, 2
data = np.random.normal(true_mu, true_sigma, size=1000)

params = maximum_likelihood_estimation(data, 'normal')
print(f"真实参数: mu={true_mu}, sigma^2={true_sigma**2}")
print(f"MLE估计: mu={params['mu']:.4f}, sigma^2={params['sigma2']:.4f}")

# 可视化似然函数
mu_range = np.linspace(3, 7, 100)
sigma_range = np.linspace(1, 4, 100)
Mu, Sigma = np.meshgrid(mu_range, sigma_range)

log_likelihood = np.zeros_like(Mu)
for i in range(len(sigma_range)):
    for j in range(len(mu_range)):
        log_likelihood[i, j] = np.sum(norm.logpdf(data, Mu[i, j], Sigma[i, j]))

plt.figure(figsize=(10, 6))
contour = plt.contour(Mu, Sigma, log_likelihood, levels=20)
plt.clabel(contour, inline=True, fontsize=8)
plt.plot(params['mu'], np.sqrt(params['sigma2']), 'r*', markersize=15, label='MLE')
plt.xlabel('mu')
plt.ylabel('sigma')
plt.title('Log-Likelihood Contours')
plt.legend()
plt.show()

def compute_fisher_information_normal(mu, sigma2, n):
    """
    计算正态分布的Fisher信息矩阵

    对于正态分布 N(mu, sigma^2):
    I(mu) = n / sigma^2
    I(sigma^2) = n / (2 * sigma^4)
    """
    I_mu = n / sigma2
    I_sigma2 = n / (2 * sigma2**2)

    return I_mu, I_sigma2

# 验证Cramér-Rao下界
n = 100
mu_true, sigma2_true = 5, 4

# 多次实验，计算样本均值的方差
n_experiments = 1000
mu_estimates = []

for _ in range(n_experiments):
    data = np.random.normal(mu_true, np.sqrt(sigma2_true), n)
    mu_estimates.append(np.mean(data))

empirical_variance = np.var(mu_estimates, ddof=1)

# 理论下界
I_mu, _ = compute_fisher_information_normal(mu_true, sigma2_true, n)
cramer_rao_bound = 1 / I_mu

print(f"样本均值的实证方差: {empirical_variance:.6f}")
print(f"Cramér-Rao下界: {cramer_rao_bound:.6f}")
print(f"MLE达到了下界: {np.isclose(empirical_variance, cramer_rao_bound, rtol=0.1)}")

def map_estimate_normal(data, prior_mu=0, prior_sigma2=1, prior_weight=1):
    """
    正态分布均值的MAP估计

    先验: mu ~ N(prior_mu, prior_sigma2)
    似然: data ~ N(mu, sigma^2)（假设sigma^2已知）
    """
    n = len(data)
    sample_mean = np.mean(data)

    # 假设sigma^2 = 1（简化）
    sigma2 = 1

    # MAP估计是高斯先验和高斯似然的加权平均
    # 后验: mu|data ~ N(mu_n, sigma_n^2)
    # 其中:
    # 1/sigma_n^2 = n/sigma^2 + 1/prior_sigma2
    # mu_n = sigma_n^2 * (n*sample_mean/sigma^2 + prior_mu/prior_sigma2)

    posterior_precision = n / sigma2 + 1 / prior_sigma2
    posterior_variance = 1 / posterior_precision

    posterior_mean = posterior_variance * (n * sample_mean / sigma2 + prior_mu / prior_sigma2)

    return posterior_mean, posterior_variance

# 示例：小样本情况下的MAP优势
np.random.seed(42)
true_mu = 5
n_samples = 5  # 小样本

data = np.random.normal(true_mu, 1, n_samples)

# MLE
mle = np.mean(data)

# MAP（强先验：认为mu在0附近）
map_est, map_var = map_estimate_normal(data, prior_mu=0, prior_sigma2=1)

# MAP（弱先验）
map_est_weak, _ = map_estimate_normal(data, prior_mu=0, prior_sigma2=10)

print(f"真实值: {true_mu}")
print(f"MLE: {mle:.4f}")
print(f"MAP（强先验）: {map_est:.4f}")
print(f"MAP（弱先验）: {map_est_weak:.4f}")
print(f"\n强先验将估计向0拉近，起到正则化作用")

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 可视化梯度
def f(x, y):
    return x**2 + 2*y**2

def grad_f(x, y):
    return np.array([2*x, 4*y])

# 绘制函数和梯度场
x = np.linspace(-3, 3, 20)
y = np.linspace(-3, 3, 20)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

U, V = grad_f(X, Y)

plt.figure(figsize=(10, 8))
contour = plt.contour(X, Y, Z, levels=15)
plt.clabel(contour, inline=True, fontsize=8)
plt.quiver(X, Y, U, V, alpha=0.6)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Function Contours and Gradient Field')
plt.colorbar(contour)
plt.show()

# 梯度指向函数增长最快的方向，垂直于等高线

import numpy as np

class ComputationalGraph:
    """
    计算图和反向传播的简单实现
    """
    def __init__(self):
        self.gradients = {}

    def forward(self, x, y, z):
        """
        计算: f = (x + y) * z
        """
        self.x = x
        self.y = y
        self.z = z

        self.q = x + y  # 中间变量
        self.f = self.q * z

        return self.f

    def backward(self, df=1):
        """
        反向传播计算梯度

        f = q * z
        q = x + y

        df/df = 1
        df/dq = z
        df/dz = q
        df/dx = df/dq * dq/dx = z * 1 = z
        df/dy = df/dq * dq/dy = z * 1 = z
        """
        # 从输出开始
        df_dq = self.z
        df_dz = self.q

        # 传播到输入
        df_dx = df_dq * 1  # dq/dx = 1
        df_dy = df_dq * 1  # dq/dy = 1

        self.gradients = {
            'x': df_dx,
            'y': df_dy,
            'z': df_dz
        }

        return self.gradients

# 示例
graph = ComputationalGraph()
f = graph.forward(x=2, y=3, z=4)
grads = graph.backward()  # 反向传播计算梯度

print(f"f = (2 + 3) * 4 = {f}")
print(f"df/dx = {grads['x']}")
print(f"df/dy = {grads['y']}")
print(f"df/dz = {grads['z']}")

# 数值验证
eps = 1e-5
f_x_plus = graph.forward(x=2+eps, y=3, z=4)
grad_x_numerical = (f_x_plus - f) / eps
print(f"\n数值验证 df/dx: {grad_x_numerical}")