建模资源中心

🔢

∫

🎲

考虑矩阵 A = [[2, 0], [0, 3]]

这是一个伸缩变换矩阵：

• x轴方向（向量[1,0]）被拉伸2倍 → 特征值λ₁=2

• y轴方向（向量[0,1]）被拉伸3倍 → 特征值λ₂=3

单位圆经过变换后变成椭圆，长轴沿y方向（因为3>2）

import numpy as np

# 创建矩阵
A = np.array([[4, 1],
              [2, 3]])

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

print("特征值:", eigenvalues)
# 输出: [5. 2.]

print("特征向量（每列是一个特征向量）:")
print(eigenvectors)
# 第一列是对应λ=5的特征向量
# 第二列是对应λ=2的特征向量

# 验证: Av = λv
v1 = eigenvectors[:, 0]  # 第一个特征向量
print("Av =", A @ v1)
print("λv =", eigenvalues[0] * v1)
# 两者应该相等！

Aₖ = σ₁u₁v₁ᵀ + σ₂u₂v₂ᵀ + ... + σₖuₖvₖᵀ

当k远小于原矩阵秩时，可以用很少的空间近似原矩阵，实现数据压缩和去噪。

import numpy as np
import matplotlib.pyplot as plt

# 创建示例矩阵（可以想象成一张灰度图像）
np.random.seed(42)
A = np.random.randn(100, 80)

# SVD分解
U, S, Vt = np.linalg.svd(A, full_matrices=False)

print("原始矩阵形状:", A.shape)
print("奇异值个数:", len(S))
print("前5个奇异值:", S[:5])

# 低秩近似
def low_rank_approx(A, k):
    U, S, Vt = np.linalg.svd(A, full_matrices=False)
    return U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]

# 用不同秩近似
for k in [1, 5, 10, 50]:
    A_k = low_rank_approx(A, k)
    error = np.linalg.norm(A - A_k)
    print(f"秩{k}近似: 误差 = {error:.4f}")

# 压缩率
original_size = 100 * 80
compressed_size = 100 * 10 + 10 + 80 * 10  # k=10
print(f"\n压缩率: {compressed_size/original_size*100:.1f}%")

AᵀAx̂ = Aᵀb

如果AᵀA可逆，则最小二乘解为：

x̂ = (AᵀA)⁻¹Aᵀb

设计矩阵 X = [[1, x₁], [1, x₂], ..., [1, xₙ]]ᵀ
响应向量 y = [y₁, y₂, ..., yₙ]ᵀ

最小二乘估计：β̂ = (XᵀX)⁻¹Xᵀy

import numpy as np
from sklearn.linear_model import LinearRegression

# 生成数据
np.random.seed(42)
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = 2 * x + 1 + np.random.normal(0, 1, 10)  # y = 2x + 1 + 噪声

# 方法1：手动计算（正规方程）
X = np.column_stack([np.ones(len(x)), x])  # 添加截距列
beta = np.linalg.inv(X.T @ X) @ X.T @ y
print(f"手动计算: β₀={beta[0]:.3f}, β₁={beta[1]:.3f}")

# 方法2：使用numpy最小二乘
beta_np, residuals, rank, s = np.linalg.lstsq(X, y, rcond=None)
print(f"numpy lstsq: β₀={beta_np[0]:.3f}, β₁={beta_np[1]:.3f}")

# 方法3：使用sklearn
model = LinearRegression()
model.fit(x.reshape(-1, 1), y)
print(f"sklearn: β₀={model.intercept_:.3f}, β₁={model.coef_[0]:.3f}")

# 计算R²（决定系数）
y_pred = model.predict(x.reshape(-1, 1))
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r_squared = 1 - ss_res / ss_tot
print(f"R² = {r_squared:.4f}")

∇f = (∂f/∂x₁, ∂f/∂x₂, ..., ∂f/∂xₙ)

想象你站在一座山上：

• 梯度方向就是你感觉"最陡的上坡方向"
• 梯度的大小就是坡度的陡峭程度
• 如果你想最快下山，应该沿梯度的反方向走

算法步骤：
1. 初始化起始点 x₀
2. 计算当前位置的梯度 ∇f(xₖ)
3. 沿梯度反方向更新：xₖ₊₁ = xₖ - η·∇f(xₖ)
4. 重复步骤2-3，直到梯度接近0

其中 η 是学习率（步长）：
• η 太大：可能震荡或发散
• η 太小：收敛速度慢

import numpy as np

# 目标函数：f(x,y) = x² + 2y² （椭圆抛物面）
def f(x, y):
    return x**2 + 2*y**2

# 梯度：∇f = (2x, 4y)
def gradient(x, y):
    return np.array([2*x, 4*y])

# 梯度下降算法
def gradient_descent(start_x, start_y, learning_rate=0.1, epochs=100):
    x, y = start_x, start_y
    path = [(x, y)]  # 记录路径
    
    for i in range(epochs):
        grad = gradient(x, y)
        x = x - learning_rate * grad[0]
        y = y - learning_rate * grad[1]
        path.append((x, y))
        
        # 检查收敛
        if np.linalg.norm(grad) < 1e-6:
            print(f"第{i+1}次迭代收敛")
            break
    
    return x, y, path

# 运行
start_x, start_y = 3.0, 3.0
opt_x, opt_y, path = gradient_descent(start_x, start_y, 0.1, 50)

print(f"起始点: ({start_x}, {start_y})")
print(f"最优点: ({opt_x:.4f}, {opt_y:.4f})")
print(f"最小值: f({opt_x:.4f}, {opt_y:.4f}) = {f(opt_x, opt_y):.6f}")
print(f"迭代次数: {len(path)-1}")

# 查看前5次迭代
print("\n前5次迭代过程:")
for i, (x, y) in enumerate(path[:5]):
    print(f"  迭代{i}: x={x:.3f}, y={y:.3f}, f(x,y)={f(x,y):.3f}")

min f(x)，满足 g(x) = 0

L(x, λ) = f(x) + λ·g(x)

∂L/∂x = 0 → ∇f(x) + λ·∇g(x) = 0
∂L/∂λ = 0 → g(x) = 0

问题：在圆 x² + y² = 1 上找一点，使 f(x,y) = x + 2y 最大

解：

1. 构造拉格朗日函数：L = x + 2y + λ(x² + y² - 1)
2. 求偏导：
   • ∂L/∂x = 1 + 2λx = 0
   • ∂L/∂y = 2 + 2λy = 0
   • ∂L/∂λ = x² + y² - 1 = 0
3. 解得：x = 1/√5, y = 2/√5, λ = -√5/2
4. 最大值：f = x + 2y = √5 ≈ 2.236

from scipy.optimize import minimize
import numpy as np

# 目标函数：f(x,y) = x + 2y
def objective(x):
    return -(x[0] + 2*x[1])  # 负号因为scipy求最小值

# 约束条件：x² + y² = 1
def constraint(x):
    return x[0]**2 + x[1]**2 - 1

# 定义约束
con = {'type': 'eq', 'fun': constraint}

# 初始猜测
x0 = [0.5, 0.5]

# 求解
result = minimize(objective, x0, constraints=con)

print(f"最优解: x={result.x[0]:.4f}, y={result.x[1]:.4f}")
print(f"最大值: f = {result.x[0] + 2*result.x[1]:.4f}")
print(f"理论值: √5 = {np.sqrt(5):.4f}")
print(f"约束满足: x²+y² = {result.x[0]**2 + result.x[1]**2:.6f}")

dy/dt = f(t, y)，初始条件 y(t₀) = y₀

迭代公式：
yₙ₊₁ = yₙ + h·f(tₙ, yₙ)

其中 h 是步长
• 优点：简单易懂
• 缺点：精度低，需要很小步长
• 误差：O(h)

k₁ = f(tₙ, yₙ)
k₂ = f(tₙ + h/2, yₙ + h·k₁/2)
k₃ = f(tₙ + h/2, yₙ + h·k₂/2)
k₄ = f(tₙ + h, yₙ + h·k₃)

yₙ₊₁ = yₙ + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)

误差：O(h⁴)，精度高

(1) 指数增长模型：

dy/dt = ky，解：y = y₀eᵏᵗ

应用：人口增长、放射性衰变、复利计算

(2) Logistic增长模型：

dP/dt = rP(1 - P/K)

应用：有上限的人口增长、传染病传播、产品推广

(3) Lorenz系统（混沌）：

dx/dt = σ(y-x)
dy/dt = x(ρ-z) - y
dz/dt = xy - βz

应用：气象预报、混沌理论、蝴蝶效应

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# === 例1：Logistic增长 ===
def logistic(P, t, r, K):
    return r * P * (1 - P / K)

# 参数
r = 0.5    # 增长率
K = 100    # 环境容纳量
P0 = 5     # 初始种群
t = np.linspace(0, 20, 100)

# 求解
solution = odeint(logistic, P0, t, args=(r, K))

print("Logistic增长模型结果:")
print(f"初始种群: {P0}")
print(f"10时刻种群: {solution[50,0]:.2f}")
print(f"20时刻种群: {solution[-1,0]:.2f} (接近K={K})")

# === 例2：Lorenz混沌系统 ===
def lorenz(state, t, sigma, rho, beta):
    x, y, z = state
    return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]

# 参数和初始条件
sigma, rho, beta = 10, 28, 8/3
state0 = [0.1, 0, 0]
t = np.linspace(0, 50, 1000)

# 求解
solution = odeint(lorenz, state0, t, args=(sigma, rho, beta))

print("\nLorenz系统（混沌）:")
print(f"参数: σ={sigma}, ρ={rho}, β={beta:.2f}")
print(f"最终状态: x={solution[-1,0]:.3f}, y={solution[-1,1]:.3f}, z={solution[-1,2]:.3f}")
print("系统呈现混沌行为：初值敏感性（蝴蝶效应）")

设X₁, X₂, ..., Xₙ是独立同分布的随机变量，具有有限的期望E[Xᵢ] = μ和方差Var(Xᵢ) = σ²。

当n足够大时，样本均值的标准化变量近似服从标准正态分布：

Z = (X̄ - μ) / (σ/√n) ~ N(0, 1)

假设掷骰子服从均匀分布U(1,6)：

• 掷1次：均匀分布，每个面概率1/6
• 掷30次取平均：近似正态分布N(3.5, σ²/30)
• 掷1000次取平均：更接近正态分布，方差更小

结论：单个骰子是均匀分布，但大量骰子的平均值是正态分布！

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

# 均匀分布（掷骰子）
np.random.seed(42)
n_samples = 10000  # 实验次数
n_dice = 30        # 每次掷30个骰子

# 生成数据
dice_rolls = np.random.randint(1, 7, size=(n_samples, n_dice))
sample_means = np.mean(dice_rolls, axis=1)

# 理论值
mu = 3.5  # E[X] = (1+2+3+4+5+6)/6 = 3.5
sigma = np.sqrt(35/12) / np.sqrt(n_dice)  # σ/√n

print("中心极限定理验证:")
print(f"骰子原始分布: 均匀分布U(1,6)")
print(f"理论均值: {mu}")
print(f"样本均值的均值: {np.mean(sample_means):.4f}")
print(f"理论标准差: {sigma:.4f}")
print(f"样本均值的标准差: {np.std(sample_means):.4f}")

# 正态性检验
stat, p_value = stats.normaltest(sample_means)
print(f"\n正态性检验 p值: {p_value:.6f}")
print(f"结论: {'服从正态分布' if p_value > 0.05 else '不服从正态分布'}")

# 3σ原则验证
within_1sigma = np.sum(np.abs(sample_means - mu) < sigma) / n_samples
within_2sigma = np.sum(np.abs(sample_means - mu) < 2*sigma) / n_samples
within_3sigma = np.sum(np.abs(sample_means - mu) < 3*sigma) / n_samples

print(f"\n3σ原则验证:")
print(f"在μ±σ内: {within_1sigma:.2%} (理论: 68.27%)")
print(f"在μ±2σ内: {within_2sigma:.2%} (理论: 95.45%)")
print(f"在μ±3σ内: {within_3sigma:.2%} (理论: 99.73%)")

1. 建立假设：
   原假设 H₀：μ = μ₀（通常是无效果、无差异）
   备择假设 H₁：μ ≠ μ₀（双侧检验）
            或 μ > μ₀ / μ < μ₀（单侧检验）

2. 选择检验统计量：
   • 方差已知或大样本：Z统计量
   • 方差未知且小样本：t统计量

3. 计算p值：
   p值 = 在H₀成立的前提下，观察到当前或更极端结果的概率

4. 做出决策：
   • p < α（通常α=0.05）→ 拒绝H₀
   • p ≥ α → 不拒绝H₀

p值不是"原假设为真的概率"！

p值正确的理解：

• "如果原假设成立，观察到当前样本（或更极端样本）的概率"
• p值小 → 当前样本与原假设不符 → 拒绝原假设
• p值大 → 当前样本与原假设一致 → 不拒绝原假设

	不拒绝H₀	拒绝H₀
H₀为真	✓ 正确	Ⅰ类错误（概率α）
H₀为假	Ⅱ类错误（概率β）	✓ 正确

import numpy as np
from scipy import stats

# === 例1：单样本t检验 ===
# 问题：某药是否有效？（血压均值是否不等于120）
np.random.seed(42)
blood_pressure = np.random.normal(118, 12, 50)  # 50个样本

# t检验
t_stat, p_value = stats.ttest_1samp(blood_pressure, 120)

print("=== 单样本t检验 ===")
print(f"样本均值: {blood_pressure.mean():.2f}")
print(f"假设的总体均值: 120")
print(f"t统计量: {t_stat:.4f}")
print(f"p值: {p_value:.4f}")
print(f"结论: {'拒绝H₀，血压不等于120' if p_value < 0.05 else '不拒绝H₀，证据不足以说明血压偏离120'}")

# === 例2：独立样本t检验 ===
# 问题：新药 vs 旧药，效果是否有差异？
control = np.random.normal(100, 15, 50)    # 对照组
treatment = np.random.normal(108, 15, 50)  # 实验组

t_stat2, p_value2 = stats.ttest_ind(control, treatment)

print("\n=== 独立样本t检验 ===")
print(f"对照组均值: {control.mean():.2f}")
print(f"实验组均值: {treatment.mean():.2f}")
print(f"t统计量: {t_stat2:.4f}")
print(f"p值: {p_value2:.6f}")
print(f"结论: {'拒绝H₀，两组有显著差异，新药有效' if p_value2 < 0.05 else '不拒绝H₀，两组无显著差异'}")

# === 例3：卡方检验（分类变量独立性） ===
# 问题：吸烟与肺癌是否相关？
observed = np.array([[60, 40],   # 吸烟: 患癌60, 未患癌40
                     [30, 70]])  # 不吸烟: 患癌30, 未患癌70

chi2, p_value3, dof, expected = stats.chi2_contingency(observed)

print("\n=== 卡方检验 ===")
print(f"卡方统计量: {chi2:.4f}")
print(f"自由度: {dof}")
print(f"p值: {p_value3:.6f}")
print(f"结论: {'拒绝H₀，吸烟与肺癌相关' if p_value3 < 0.05 else '不拒绝H₀，吸烟与肺癌无关'}")

95%置信区间的意思是：如果重复抽样100次，构造100个置信区间，大约有95个区间会包含真实的总体参数。

注意：不是"参数有95%概率落在这个区间内"（参数是固定的，不是随机的）

方差已知或大样本（Z区间）：
   X̄ ± Z_α/2 · (σ/√n)

方差未知且小样本（t区间）：
   X̄ ± t_{α/2, n-1} · (S/√n)

其中：
   • Z_0.025 = 1.96（95%置信水平）
   • t值取决于样本量（自由度n-1）

import numpy as np
from scipy import stats
from sklearn.linear_model import LinearRegression

# === 例1：总体均值的置信区间 ===
np.random.seed(42)
data = np.random.normal(100, 15, 30)  # 30个样本

x_bar = np.mean(data)
s = np.std(data, ddof=1)  # 样本标准差
n = len(data)

# t临界值 (95%置信水平)
t_critical = stats.t.ppf(0.975, df=n-1)

# 置信区间
margin_of_error = t_critical * s / np.sqrt(n)
ci_lower = x_bar - margin_of_error
ci_upper = x_bar + margin_of_error

print("=== 总体均值的95%置信区间 ===")
print(f"样本均值: {x_bar:.2f}")
print(f"样本标准差: {s:.2f}")
print(f"样本量: {n}")
print(f"t临界值: {t_critical:.4f}")
print(f"误差范围: ±{margin_of_error:.2f}")
print(f"95%置信区间: [{ci_lower:.2f}, {ci_upper:.2f}]")
print(f"区间是否包含真实均值100? {ci_lower < 100 < ci_upper}")

# === 例2：回归系数的置信区间 ===
# 生成数据
X = np.linspace(0, 10, 50).reshape(-1, 1)
y = 2 * X.ravel() + 1 + np.random.normal(0, 2, 50)

# 拟合回归
model = LinearRegression()
model.fit(X, y)

# 计算回归系数的标准误
y_pred = model.predict(X)
residuals = y - y_pred
MSE = np.sum(residuals**2) / (n - 2)
SE_beta1 = np.sqrt(MSE / np.sum((X - np.mean(X))**2))

# 置信区间
t_crit = stats.t.ppf(0.975, df=n-2)
beta1_ci = [model.coef_[0] - t_crit*SE_beta1, 
            model.coef_[0] + t_crit*SE_beta1]

print("\n=== 回归系数的95%置信区间 ===")
print(f"斜率估计: {model.coef_[0]:.3f}")
print(f"95%置信区间: [{beta1_ci[0]:.3f}, {beta1_ci[1]:.3f}]")
print(f"区间是否包含0? {beta1_ci[0] < 0 < beta1_ci[1]}")
print(f"结论: {'斜率不显著' if beta1_ci[0] < 0 < beta1_ci[1] else '斜率显著不为0'}")