import numpy as np

def polynomial_fit(x, t, M):
    """多项式曲线拟合的最小二乘解"""
    N = len(x)
    # 构造设计矩阵
    Phi = np.zeros((N, M + 1))
    for j in range(M + 1):
        Phi[:, j] = x ** j
    # 最小二乘解
    w = np.linalg.inv(Phi.T @ Phi) @ Phi.T @ t
    return w

def polynomial_predict(x, w):
    """使用拟合的多项式进行预测"""
    M = len(w) - 1
    y = np.zeros_like(x)
    for j in range(M + 1):
        y += w[j] * x ** j
    return y

# 生成示例数据
np.random.seed(42)
N = 10
x = np.linspace(0, 1, N)
t = np.sin(2 * np.pi * x) + np.random.normal(0, 0.3, N)

# 不同阶数的多项式拟合
for M in [1, 3, 9]:
    w = polynomial_fit(x, t, M)
    print(f"M={M}, 权重: {np.round(w, 3)}")

import matplotlib.pyplot as plt

def evaluate_model(x_train, t_train, x_test, t_test, M):
    """评估不同阶数多项式的训练误差和测试误差"""
    w = polynomial_fit(x_train, t_train, M)
    t_train_pred = polynomial_predict(x_train, w)
    t_test_pred = polynomial_predict(x_test, w)
    
    train_error = np.sqrt(np.mean((t_train_pred - t_train) ** 2))
    test_error = np.sqrt(np.mean((t_test_pred - t_test) ** 2))
    return train_error, test_error

# 生成训练数据和测试数据
np.random.seed(42)
N_train = 10
x_train = np.linspace(0, 1, N_train)
t_train = np.sin(2 * np.pi * x_train) + np.random.normal(0, 0.3, N_train)

N_test = 100
x_test = np.linspace(0, 1, N_test)
t_test = np.sin(2 * np.pi * x_test) + np.random.normal(0, 0.3, N_test)

print(f"{'M':>3} {'训练误差':>10} {'测试误差':>10}")
print("-" * 30)
for M in range(10):
    train_err, test_err = evaluate_model(x_train, t_train, x_test, t_test, M)
    print(f"{M:3d} {train_err:10.4f} {test_err:10.4f}")

def bayesian_update(prior, likelihood, evidence):
    """贝叶斯更新：后验 = 似然 * 先验 / 证据"""
    posterior = likelihood * prior / evidence
    return posterior

# 示例：硬币抛掷问题
# 假设我们有一枚硬币，不确定它是否公平
# 先验：认为硬币是公平的（正面概率 theta=0.5 的概率最高）
import numpy as np
from scipy import stats

theta = np.linspace(0, 1, 100)
# 先验：Beta(2, 2)，倾向于 theta=0.5
prior = stats.beta.pdf(theta, 2, 2)
# 观察数据：10次抛掷中7次正面
# 似然：Binomial(10, theta) 在 k=7 时的值
likelihood = stats.binom.pmf(7, 10, theta)
# 证据：归一化常数
evidence = np.trapz(likelihood * prior, theta)
# 后验
posterior = bayesian_update(prior, likelihood, evidence)

print(f"后验均值: {np.trapz(theta * posterior, theta):.3f}")
print(f"后验众数: {theta[np.argmax(posterior)]:.3f}")

def gaussian_pdf(x, mu, sigma):
    """一元高斯分布概率密度函数"""
    return (1 / np.sqrt(2 * np.pi * sigma**2)) * np.exp(-(x - mu)**2 / (2 * sigma**2))

def multivariate_gaussian_pdf(x, mu, Sigma):
    """多元高斯分布概率密度函数"""
    D = len(mu)
    diff = x - mu
    normalization = 1 / ((2 * np.pi)**(D/2) * np.linalg.det(Sigma)**0.5)
    exponent = -0.5 * diff.T @ np.linalg.inv(Sigma) @ diff
    return normalization * np.exp(exponent)

# 生成高斯分布数据
mu = 0.0
sigma = 1.0
x = np.linspace(-4, 4, 100)
y = gaussian_pdf(x, mu, sigma)
print(f"高斯分布 N(0,1) 在 x=0 处的密度: {gaussian_pdf(0, 0, 1):.4f}")
print(f"高斯分布 N(0,1) 在 x=1 处的密度: {gaussian_pdf(1, 0, 1):.4f}")

def entropy(p):
    """计算离散分布的熵"""
    return -np.sum(p * np.log2(p + 1e-10))

# 示例：比较不同分布的熵
p_uniform = np.array([0.25, 0.25, 0.25, 0.25])
p_biased = np.array([0.7, 0.1, 0.1, 0.1])
p_extreme = np.array([0.97, 0.01, 0.01, 0.01])

print(f"均匀分布的熵: {entropy(p_uniform):.4f} bits")
print(f"有偏分布的熵: {entropy(p_biased):.4f} bits")
print(f"极端分布的熵: {entropy(p_extreme):.4f} bits")

def kl_divergence(p, q):
    """计算KL散度 KL(p||q)"""
    return np.sum(p * np.log((p + 1e-10) / (q + 1e-10)))

# 示例
p = np.array([0.4, 0.3, 0.2, 0.1])
q1 = np.array([0.25, 0.25, 0.25, 0.25])  # 均匀分布
q2 = np.array([0.35, 0.30, 0.25, 0.10])  # 接近p的分布

print(f"KL(p||q1) = {kl_divergence(p, q1):.4f}")
print(f"KL(p||q2) = {kl_divergence(p, q2):.4f}")

def polynomial_fit_regularized(x, t, M, lam):
    """带L2正则化的多项式拟合（岭回归）"""
    N = len(x)
    Phi = np.zeros((N, M + 1))
    for j in range(M + 1):
        Phi[:, j] = x ** j
    # 岭回归的解析解
    I = np.eye(M + 1)
    I[0, 0] = 0  # 不正则化偏置项
    w = np.linalg.inv(Phi.T @ Phi + lam * I) @ Phi.T @ t
    return w

# 比较不同正则化系数的效果
M = 9
for lam in [0, 0.001, 0.01, 0.1, 1.0]:
    w = polynomial_fit_regularized(x_train, t_train, M, lam)
    train_err = np.sqrt(np.mean((polynomial_predict(x_train, w) - t_train) ** 2))
    test_err = np.sqrt(np.mean((polynomial_predict(x_test, w) - t_test) ** 2))
    print(f"λ={lam:.3f}, 训练误差={train_err:.4f}, 测试误差={test_err:.4f}")

def k_fold_cross_validation(x, t, M, lam, K=5):
    """K折交叉验证"""
    N = len(x)
    indices = np.arange(N)
    np.random.shuffle(indices)
    fold_size = N // K
    
    errors = []
    for k in range(K):
        # 划分训练集和验证集
        val_start = k * fold_size
        val_end = val_start + fold_size if k < K - 1 else N
        val_idx = indices[val_start:val_end]
        train_idx = np.concatenate([indices[:val_start], indices[val_end:]])
        
        x_tr, t_tr = x[train_idx], t[train_idx]
        x_val, t_val = x[val_idx], t[val_idx]
        
        # 训练并评估
        w = polynomial_fit_regularized(x_tr, t_tr, M, lam)
        t_pred = polynomial_predict(x_val, w)
        error = np.sqrt(np.mean((t_pred - t_val) ** 2))
        errors.append(error)
    
    return np.mean(errors)

# 使用交叉验证选择最优正则化系数
lambdas = [0.0001, 0.001, 0.01, 0.1, 1.0]
best_lambda = None
best_error = float('inf')

for lam in lambdas:
    cv_error = k_fold_cross_validation(x_train, t_train, M=9, lam=lam, K=5)
    print(f"λ={lam:.4f}, CV误差={cv_error:.4f}")
    if cv_error < best_error:
        best_error = cv_error
        best_lambda = lam

print(f"\n最优正则化系数: λ={best_lambda:.4f}")

def bayesian_linear_regression(x_train, t_train, x_test, alpha, beta):
    """贝叶斯线性回归"""
    N = len(x_train)
    M = 9  # 多项式阶数
    
    # 构造设计矩阵
    Phi_train = np.column_stack([x_train ** j for j in range(M + 1)])
    Phi_test = np.column_stack([x_test ** j for j in range(M + 1)])
    
    # 后验参数
    S_N_inv = alpha * np.eye(M + 1) + beta * Phi_train.T @ Phi_train
    S_N = np.linalg.inv(S_N_inv)
    m_N = beta * S_N @ Phi_train.T @ t_train
    
    # 预测分布
    pred_mean = Phi_test @ m_N
    pred_var = np.array([1/beta + Phi_test[i] @ S_N @ Phi_test[i] for i in range(len(x_test))])
    
    return pred_mean, pred_var

# 示例
x_test_fine = np.linspace(0, 1, 200)
pred_mean, pred_var = bayesian_linear_regression(x_train, t_train, x_test_fine, alpha=2.0, beta=25.0)
pred_std = np.sqrt(pred_var)

print(f"预测均值范围: [{pred_mean.min():.3f}, {pred_mean.max():.3f}]")
print(f"预测标准差范围: [{pred_std.min():.3f}, {pred_std.max():.3f}]")

import numpy as np
import matplotlib.pyplot as plt

def run_curve_fitting_experiment():
    """完整的曲线拟合实验：展示过拟合、正则化和贝叶斯方法"""
    
    # 1. 数据生成
    np.random.seed(42)
    N = 10
    x = np.linspace(0, 1, N)
    true_function = np.sin(2 * np.pi * x)
    noise = np.random.normal(0, 0.3, N)
    t = true_function + noise
    
    # 2. 不同阶数的多项式拟合
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    x_fine = np.linspace(0, 1, 100)
    
    for idx, M in enumerate([0, 1, 3, 9]):
        ax = axes[idx // 2, idx % 2]
        w = polynomial_fit(x, t, M)
        y_fine = polynomial_predict(x_fine, w)
        
        ax.plot(x_fine, np.sin(2 * np.pi * x_fine), 'g--', label='True function', linewidth=2)
        ax.plot(x_fine, y_fine, 'b-', label=f'M={M} polynomial', linewidth=2)
        ax.scatter(x, t, c='red', s=50, zorder=5, label='Data points')
        ax.set_title(f'M={M}')
        ax.set_ylim(-1.5, 1.5)
        ax.legend(fontsize=8)
        ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('curve_fitting.png', dpi=150)
    plt.close()
    
    # 3. 偏差-方差分析
    n_trials = 100
    M_values = range(10)
    bias_squared = []
    variance = []
    
    for M in M_values:
        predictions = []
        for trial in range(n_trials):
            x_trial = np.linspace(0, 1, N)
            t_trial = np.sin(2 * np.pi * x_trial) + np.random.normal(0, 0.3, N)
            w = polynomial_fit(x_trial, t_trial, M)
            predictions.append(polynomial_predict(x_fine, w))
        
        predictions = np.array(predictions)
        mean_prediction = predictions.mean(axis=0)
        true_values = np.sin(2 * np.pi * x_fine)
        
        bias_sq = np.mean((mean_prediction - true_values) ** 2)
        var = np.mean(np.var(predictions, axis=0))
        
        bias_squared.append(bias_sq)
        variance.append(var)
    
    print("偏差-方差分析:")
    print(f"{'M':>3} {'Bias²':>10} {'Variance':>10} {'Total':>10}")
    for M, b, v in zip(M_values, bias_squared, variance):
        print(f"{M:3d} {b:10.4f} {v:10.4f} {b+v:10.4f}")
    
    return bias_squared, variance

bias_sq, var = run_curve_fitting_experiment()

def bayesian_model_comparison(x, t):
    """使用贝叶斯因子比较不同阶数的多项式模型"""
    
    M_values = range(6)
    log_evidences = []
    
    for M in M_values:
        N = len(x)
        alpha = 0.005
        beta = 25.0
        
        Phi = np.column_stack([x ** j for j in range(M + 1)])
        
        # 后验协方差
        S_N_inv = alpha * np.eye(M + 1) + beta * Phi.T @ Phi
        S_N = np.linalg.inv(S_N_inv)
        m_N = beta * S_N @ Phi.T @ t
        
        # 对数边缘似然（证据）
        # ln p(t|alpha,beta) = (M+1)/2 * ln(alpha) + N/2 * ln(beta) - E(m_N) 
        #                      - 1/2 * ln|S_N_inv| - N/2 * ln(2*pi)
        E_mN = beta / 2 * np.sum((t - Phi @ m_N) ** 2) + alpha / 2 * m_N @ m_N
        log_evidence = (M + 1) / 2 * np.log(alpha) + N / 2 * np.log(beta) \
                       - E_mN - 0.5 * np.log(np.linalg.det(S_N_inv)) \
                       - N / 2 * np.log(2 * np.pi)
        
        log_evidences.append(log_evidence)
        print(f"M={M}, ln(evidence)={log_evidence:.2f}")
    
    # 找到最优模型
    best_M = M_values[np.argmax(log_evidences)]
    print(f"\n贝叶斯模型选择结果: 最优 M = {best_M}")
    
    return log_evidences

# 运行模型比较
np.random.seed(42)
x = np.linspace(0, 1, 10)
t = np.sin(2 * np.pi * x) + np.random.normal(0, 0.3, 10)
log_evidences = bayesian_model_comparison(x, t)

def mutual_information_feature_selection(X, y, feature_names):
    """使用互信息进行特征选择"""
    n_features = X.shape[1]
    mi_scores = []
    
    for i in range(n_features):
        # 离散化连续特征
        x_binned = np.digitize(X[:, i], bins=np.linspace(X[:, i].min(), X[:, i].max(), 20))
        
        # 计算互信息
        # I(X;Y) = H(X) + H(Y) - H(X,Y)
        px = np.bincount(x_binned) / len(x_binned)
        py = np.bincount(y) / len(y)
        pxy = np.zeros((len(np.unique(x_binned)), len(np.unique(y))))
        
        for xi, yi in zip(x_binned, y):
            pxy[xi, yi] += 1
        pxy /= pxy.sum()
        
        hx = entropy(px)
        hy = entropy(py)
        hxy = entropy(pxy.flatten()[pxy.flatten() > 0])
        
        mi = hx + hy - hxy
        mi_scores.append(mi)
        print(f"特征 '{feature_names[i]}': MI = {mi:.4f} bits")
    
    # 按互信息排序
    sorted_indices = np.argsort(mi_scores)[::-1]
    print("\n特征重要性排序:")
    for rank, idx in enumerate(sorted_indices):
        print(f"  {rank+1}. {feature_names[idx]} (MI={mi_scores[idx]:.4f})")
    
    return mi_scores

# 示例：使用鸢尾花数据集
from sklearn.datasets import load_iris
iris = load_iris()
X = iris.data
y = iris.target
feature_names = iris.feature_names

mi_scores = mutual_information_feature_selection(X, y, feature_names)