📐 数学基础教程

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar

# 绘制函数及其导数，直观看单调区间
x = np.linspace(-3, 3, 400)

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# 一次函数
axes[0].plot(x, 2*x - 1, color='#b85c2a', linewidth=2, label='f(x)=2x-1')
axes[0].axhline(0, color='#3d352a', linewidth=0.8)
axes[0].axvline(0, color='#3d352a', linewidth=0.8)
axes[0].set_title('一次函数 (单调)'); axes[0].legend()

# 二次函数：顶点、对称轴
y_quad = x**2 - 2*x - 1
axes[1].plot(x, y_quad, color='#2c5282', linewidth=2, label='f(x)=x²-2x-1')
axes[1].axvline(1, color='gray', linestyle='--', linewidth=1, label='对称轴 x=1')
axes[1].scatter([1], [-2], color='red', zorder=5, label=f'顶点(1,-2)')
axes[1].set_title('二次函数 (先减后增)'); axes[1].legend(fontsize=8)

# Logistic 函数
def logistic(x, L=1, k=1, x0=0):
    return L / (1 + np.exp(-k*(x - x0)))

y_lg = logistic(x)
axes[2].plot(x, y_lg, color='#553c9a', linewidth=2, label='σ(x)=1/(1+e⁻ˣ)')
axes[2].axhline(0.5, color='gray', linestyle='--', linewidth=1, label='y=0.5 (拐点)')
axes[2].axvline(0, color='gray', linestyle=':', linewidth=1)
axes[2].set_title('Logistic 函数 (S型)'); axes[2].legend(fontsize=8)

plt.tight_layout()
plt.savefig('elementary_functions.png', dpi=120)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

n = 20  # 项数

# ① 等差数列  a_n = 2 + 3(n-1)
a1, d = 2, 3
arith = [a1 + (k-1)*d for k in range(1, n+1)]
S_arith = [k*a1 + k*(k-1)//2*d for k in range(1, n+1)]

# ② 等比数列  a_n = 1 * 1.5^(n-1)
b1, q = 1, 1.5
geo = [b1 * q**(k-1) for k in range(1, n+1)]

# ③ Fibonacci 数列
fib = [1, 1]
for k in range(2, n):
    fib.append(fib[-1] + fib[-2])

# ④ 递推 Logistic（种群模型）：N_{n+1} = 2.8 * N_n * (1 - N_n)
N = [0.1]
for _ in range(n-1):
    N.append(2.8 * N[-1] * (1 - N[-1]))

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

ns = list(range(1, n+1))
axes[0,0].stem(ns, arith, basefmt=" ", linefmt='#b85c2a', markerfmt='o')
axes[0,0].set_title(f'等差数列 (d={d})  S₂₀={S_arith[-1]}')

axes[0,1].stem(ns, geo, basefmt=" ", linefmt='#2c5282', markerfmt='o')
axes[0,1].set_title(f'等比数列 (q={q})')

axes[1,0].stem(ns, fib, basefmt=" ", linefmt='#1a6e6a', markerfmt='o')
axes[1,0].set_title('Fibonacci 数列')

axes[1,1].plot(ns, N, 'o-', color='#553c9a', linewidth=1.5, markersize=5)
axes[1,1].axhline(1 - 1/2.8, color='gray', linestyle='--', label=f'稳定点 K={(1-1/2.8):.3f}')
axes[1,1].set_title('Logistic 递推种群模型 (λ=2.8)'); axes[1,1].legend()

plt.tight_layout()
plt.savefig('sequences.png', dpi=120)
plt.show()

# 数列通项与求和（符号计算）
from sympy import symbols, summation, simplify, oo, factorial, Sum

n_sym, k = symbols('n k', positive=True, integer=True)

# 等差求和：Σk = n(n+1)/2
print(summation(k, (k, 1, n_sym)))               # n*(n+1)/2

# 等比求和：Σ2^k
print(summation(2**k, (k, 0, n_sym-1)))          # 2**n - 1

# 错位相减结果验证：Σk*2^k
print(simplify(summation(k*2**k, (k, 1, n_sym)))) # (n-1)*2^(n+1)+2

import numpy as np

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

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

print("矩阵 A:")
print(A)
print("\n特征值:")
print(eigenvalues)  # [5. 2.]
print("\n特征向量（每列是一个特征向量）:")
print(eigenvectors)

# 验证: Av = λv
v1 = eigenvectors[:, 0]
print("\n验证 Av = λv:")
print("Av₁ =", A @ v1)
print("λ₁v₁ =", eigenvalues[0] * v1)

# 特征值分解
V = eigenvectors
Lambda = np.diag(eigenvalues)
A_reconstructed = V @ Lambda @ np.linalg.inv(V)
print("\n重构矩阵（应该等于A）:")
print(A_reconstructed)
print("误差:", np.linalg.norm(A - A_reconstructed))

import numpy as np

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

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

print("原始矩阵 A:")
print(A)
print("\n奇异值:")
print(S)  # 奇异值从大到小排列

# 低秩近似
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, :]

print("\n秩1近似:")
print(low_rank_approx(A, 1))

print("\n秩2近似:")
print(low_rank_approx(A, 2))

# 压缩率
original_params = A.shape[0] * A.shape[1]
compressed_params = k * (A.shape[0] + A.shape[1] + 1)
print(f"\n原始参数: {original_params}")
print(f"秩2近似参数: {compressed_params}")
print(f"压缩率: {compressed_params/original_params*100:.1f}%")

import numpy as np

# 目标函数：f(x,y) = x² + 2y²
def f(x, y):
    return x**2 + 2*y**2

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

# 梯度下降算法
def gradient_descent(start, lr=0.1, epochs=100):
    x = np.array(start, dtype=float)
    path = [x.copy()]
    
    for i in range(epochs):
        grad = gradient(x[0], x[1])
        x = x - lr * grad
        path.append(x.copy())
        
        if np.linalg.norm(grad) < 1e-6:
            print(f"第{i+1}次迭代收敛")
            break
    
    return x, path

# 运行
start = [3.0, 3.0]
opt, path = gradient_descent(start, 0.1, 50)

print(f"起始点: {start}")
print(f"最优点: ({opt[0]:.4f}, {opt[1]:.4f})")
print(f"最小值: f = {f(opt[0], opt[1]):.6f}")
print(f"迭代次数: {len(path)-1}")

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

# 定义被积函数
def f(x):
    return x**2 * np.sin(x)

# 方法1：SciPy数值积分
result, error = integrate.quad(f, 0, np.pi)
print(f"∫₀^π x²sin(x)dx = {result:.6f}")
print(f"误差估计: {error:.2e}")

# 方法2：梯形法则
x = np.linspace(0, np.pi, 1000)
y = f(x)
trap_result = np.trapz(y, x)
print(f"梯形法则结果: {trap_result:.6f}")

# 方法3：Simpson法则
simp_result = integrate.simps(y, x)
print(f"Simpson法则结果: {simp_result:.6f}")

# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
ax.fill_between(x, y, alpha=0.3, label='积分区域')
ax.plot(x, y, linewidth=2)
ax.set_xlabel('x', fontsize=14)
ax.set_ylabel('f(x)', fontsize=14)
ax.set_title(f'∫₀^π x²sin(x)dx = {result:.4f}', fontsize=16)
ax.legend()
ax.grid(True, alpha=0.3)
plt.savefig('integral_demo.png', dpi=150)
plt.show()

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

# 示例1：人口增长模型 dP/dt = rP
def population_model(t, P, r=0.02):
    return r * P

# 求解
t_span = (0, 100)
t_eval = np.linspace(0, 100, 200)
P0 = [100]  # 初始人口

sol = solve_ivp(population_model, t_span, P0, t_eval=t_eval)

# 解析解对比
P_analytical = 100 * np.exp(0.02 * t_eval)

plt.figure(figsize=(10, 6))
plt.plot(t_eval, sol.y[0], 'b-', linewidth=2, label='数值解')
plt.plot(t_eval, P_analytical, 'r--', linewidth=2, label='解析解')
plt.xlabel('时间 t', fontsize=12)
plt.ylabel('人口 P(t)', fontsize=12)
plt.title('Malthus人口模型', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('population_model.png', dpi=150)
plt.show()

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

# 求解（两个非常接近的初始条件）
t_span = (0, 50)
t_eval = np.linspace(0, 50, 5000)
y0_1 = [1.0, 1.0, 1.0]
y0_2 = [1.0001, 1.0, 1.0]  # 微小差异

sol1 = solve_ivp(lorenz, t_span, y0_1, t_eval=t_eval)
sol2 = solve_ivp(lorenz, t_span, y0_2, t_eval=t_eval)

# 绘制混沌吸引子和蝴蝶效应
fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot(sol1.y[0], sol1.y[1], sol1.y[2], 'b-', linewidth=0.5)
ax1.set_title('Lorenz吸引子', fontsize=14)
ax1.set_xlabel('x', fontsize=10)
ax1.set_ylabel('y', fontsize=10)
ax1.set_zlabel('z', fontsize=10)

ax2 = fig.add_subplot(122)
ax2.plot(t_eval, sol1.y[0], 'b-', linewidth=1, label='初始 x₀=1.0')
ax2.plot(t_eval, sol2.y[0], 'r--', linewidth=1, label='初始 x₀=1.0001')
ax2.set_title('蝴蝶效应（敏感依赖于初始条件）', fontsize=14)
ax2.set_xlabel('时间 t', fontsize=12)
ax2.set_ylabel('x(t)', fontsize=12)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lorenz_chaos.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def taylor_sin(x, n):
    """sin(x) 的 n 阶泰勒展开"""
    result = 0
    for k in range(n+1):
        term = ((-1)**k) * x**(2*k+1) / np.math.factorial(2*k+1)
        result += term
    return result

def taylor_exp(x, n):
    """eˣ 的 n 阶泰勒展开"""
    return sum(x**k / np.math.factorial(k) for k in range(n+1))

x = np.linspace(-2*np.pi, 2*np.pi, 500)
plt.figure(figsize=(10, 6))
plt.plot(x, np.sin(x), 'k-', linewidth=2, label='sin(x) 精确值')

for n in [1, 3, 5, 9]:
    y_approx = [taylor_sin(xi, n) for xi in x]
    plt.plot(x, y_approx, '--', linewidth=1.5, label=f'n={n}')

plt.axhline(0, color='gray', linewidth=0.5)
plt.axvline(0, color='gray', linewidth=0.5)
plt.title('sin(x) 的泰勒展开', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-3, 3)
plt.savefig('taylor_sin.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 方波的傅里叶级数: f(x) = 4/π * Σ sin((2k-1)x)/(2k-1)
def fourier_square(x, n_terms):
    result = 0
    for k in range(1, n_terms+1):
        result += np.sin((2*k-1)*x) / (2*k-1)
    return (4/np.pi) * result

# 锯齿波: f(x) = 2/π * Σ (-1)^(k+1) * sin(kx)/k
def fourier_sawtooth(x, n_terms):
    result = 0
    for k in range(1, n_terms+1):
        result += ((-1)**(k+1)) * np.sin(k*x) / k
    return (2/np.pi) * result

x = np.linspace(-2*np.pi, 2*np.pi, 1000)

plt.figure(figsize=(12, 5))

# 方波
plt.subplot(1, 2, 1)
plt.plot(x, np.sign(np.sin(x)), 'k--', linewidth=1, alpha=0.5, label='理想方波')
for n in [1, 3, 7, 20]:
    plt.plot(x, fourier_square(x, n), linewidth=1.5, label=f'N={n}')
plt.title('方波的傅里叶级数展开', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-1.5, 1.5)

# 锯齿波
plt.subplot(1, 2, 2)
target = (2/np.pi) * np.arctan2(np.sin(x), np.cos(x))
plt.plot(x, target, 'k--', linewidth=1, alpha=0.5, label='理想锯齿波')
for n in [1, 3, 7, 20]:
    plt.plot(x, fourier_sawtooth(x, n), linewidth=1.5, label=f'N={n}')
plt.title('锯齿波的傅里叶级数展开', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-1.5, 1.5)

plt.tight_layout()
plt.savefig('fourier_series.png', dpi=150)
plt.show()

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

# 创建子图
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. 正态分布（不同参数）
x = np.linspace(-5, 5, 500)
for mu, sigma in [(0, 1), (0, 2), (2, 1)]:
    y = stats.norm.pdf(x, mu, sigma)
    label = f'μ={mu}, σ={sigma}'
    axes[0,0].plot(x, y, linewidth=2, label=label)
axes[0,0].set_title('正态分布', fontsize=14)
axes[0,0].legend()
axes[0,0].grid(True, alpha=0.3)

# 2. 二项分布 vs Poisson近似
n, p = 50, 0.1
k = np.arange(0, 20)
binom_pmf = stats.binom.pmf(k, n, p)
poisson_pmf = stats.poisson.pmf(k, n*p)
axes[0,1].bar(k, binom_pmf, alpha=0.6, label=f'二项 B({n},{p})')
axes[0,1].plot(k, poisson_pmf, 'r-o', linewidth=2, label=f'Poisson(λ={n*p})')
axes[0,1].set_title('二项分布与Poisson近似', fontsize=14)
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)

# 3. 中心极限定理演示
sample_means = []
for _ in range(10000):
    samples = np.random.exponential(scale=2, size=100)
    sample_means.append(np.mean(samples))
axes[1,0].hist(sample_means, bins=50, density=True, alpha=0.7, label='样本均值分布')
# 理论正态
x_norm = np.linspace(0, 4, 200)
axes[1,0].plot(x_norm, stats.norm.pdf(x_norm, 2, 2/np.sqrt(100)),
               'r-', linewidth=2, label='理论正态')
axes[1,0].set_title('中心极限定理（指数分布样本均值）', fontsize=14)
axes[1,0].legend()
axes[1,0].grid(True, alpha=0.3)

# 4. 卡方分布
x = np.linspace(0, 30, 500)
for df in [1, 2, 5, 10]:
    axes[1,1].plot(x, stats.chi2.pdf(x, df), linewidth=2, label=f'df={df}')
axes[1,1].set_title('卡方分布', fontsize=14)
axes[1,1].legend()
axes[1,1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('probability_distributions.png', dpi=150)
plt.show()

import numpy as np
from scipy import stats

# 示例数据：某批零件的测量值（单位：mm）
np.random.seed(42)
data = np.random.normal(loc=50.0, scale=2.5, size=30)

# ===== 1. 矩估计 =====
mu_hat = np.mean(data)           # 均值估计
sigma_hat = np.std(data, ddof=1) # 标准差估计（无偏）

print("=== 矩估计 ===")
print(f"样本均值（μ的估计）: {mu_hat:.4f}")
print(f"样本标准差（σ的估计）: {sigma_hat:.4f}")
print(f"样本量: {len(data)}")

# ===== 2. 正态分布均值的置信区间 =====
confidence = 0.95
alpha = 1 - confidence
n = len(data)
se = sigma_hat / np.sqrt(n)  # 标准误差

# t分布临界值
t_critical = stats.t.ppf(1 - alpha/2, df=n-1)
ci_lower = mu_hat - t_critical * se
ci_upper = mu_hat + t_critical * se

print(f"\n=== {confidence*100:.0f}% 置信区间 ===")
print(f"t临界值 t_{{α/2}}(n-1) = t_{{0.025}}({n-1}) = {t_critical:.4f}")
print(f"标准误差 SE = {se:.4f}")
print(f"均值的{confidence*100:.0f}%置信区间: [{ci_lower:.4f}, {ci_upper:.4f}]")

# 验证：使用scipy直接计算
ci = stats.t.interval(confidence, df=n-1, loc=mu_hat, scale=se)
print(f"scipy验证: [{ci[0]:.4f}, {ci[1]:.4f}]")

# ===== 3. 正态分布的MLE参数估计 =====
mu_mle, sigma_mle = stats.norm.fit(data)
print(f"\n=== 极大似然估计（MLE）===")
print(f"μ_MLE = {mu_mle:.4f}（等于样本均值）")
print(f"σ_MLE = {sigma_mle:.4f}（有偏，分母为n而非n-1）")

# ===== 4. 样本量对置信区间宽度的影响 =====
print("\n=== 样本量对置信区间宽度的影响 ===")
for n_size in [10, 30, 100, 1000]:
    sample = np.random.normal(50, 2.5, n_size)
    ci_s = stats.t.interval(0.95, df=n_size-1,
                            loc=np.mean(sample),
                            scale=stats.sem(sample))
    width = ci_s[1] - ci_s[0]
    print(f"n={n_size:4d}: 区间宽度 = {width:.4f}")

import numpy as np
from scipy import stats

# 示例1：单样本t检验
# 问题：某厂生产的零件平均长度是否为10cm？
sample_data = np.array([10.1, 9.8, 10.3, 9.9, 10.0, 10.2, 9.7, 10.1, 9.8, 10.0])

# 执行t检验
t_stat, p_value = stats.ttest_1samp(sample_data, popmean=10.0)

print("=== 单样本t检验 ===")
print(f"样本均值: {np.mean(sample_data):.3f}")
print(f"样本标准差: {np.std(sample_data, ddof=1):.3f}")
print(f"t统计量: {t_stat:.4f}")
print(f"p值: {p_value:.4f}")
if p_value > 0.05:
    print("结论: p > 0.05, 不拒绝H₀, 没有足够证据说明均值不同于10cm")
else:
    print("结论: p ≤ 0.05, 拒绝H₀, 均值显著不同于10cm")

# 示例2：两样本t检验
# 问题：两种教学方法的效果是否有差异？
group_a = np.array([85, 88, 90, 82, 87, 91, 86, 89])
group_b = np.array([78, 80, 85, 76, 82, 79, 83, 81])

t_stat2, p_value2 = stats.ttest_ind(group_a, group_b)

print("\n=== 两样本t检验 ===")
print(f"A组均值: {np.mean(group_a):.2f}")
print(f"B组均值: {np.mean(group_b):.2f}")
print(f"t统计量: {t_stat2:.4f}")
print(f"p值: {p_value2:.4f}")
if p_value2 < 0.05:
    print("结论: p ≤ 0.05, 拒绝H₀, 两种教学方法效果有显著差异")
else:
    print("结论: p > 0.05, 不拒绝H₀")

# 示例3：卡方独立性检验
# 问题：吸烟与患肺癌是否有关？
observed = np.array([[40, 60],   # 吸烟: 患癌40, 未患癌60
                     [10, 90]])  # 不吸烟: 患癌10, 未患癌90

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

print("\n=== 卡方独立性检验 ===")
print(f"卡方统计量: {chi2:.4f}")
print(f"自由度: {dof}")
print(f"p值: {p_chi:.6f}")
print(f"期望频数:\n{expected}")
if p_chi < 0.05:
    print("结论: p ≤ 0.05, 拒绝H₀, 吸烟与患肺癌有关联")
else:
    print("结论: p > 0.05, 不拒绝H₀")

# 示例4：线性回归与相关性检验
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = np.array([2.3, 4.1, 6.2, 8.0, 10.1, 11.8, 14.2, 16.0, 17.9, 20.1])

slope, intercept, r_value, p_value_reg, std_err = stats.linregress(x, y)

print("\n=== 线性回归 ===")
print(f"回归方程: y = {slope:.3f}x + {intercept:.3f}")
print(f"相关系数 r = {r_value:.4f}")
print(f"R² = {r_value**2:.4f}")
print(f"p值 = {p_value_reg:.6f}")
if p_value_reg < 0.05:
    print("结论: x与y之间存在显著的线性关系")

import numpy as np
from scipy import stats

# 示例数据：学习时间与考试成绩
study_hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
exam_scores = np.array([52, 55, 58, 62, 68, 72, 78, 82, 88, 95])

# 方法1：手动计算
x_mean = np.mean(study_hours)
y_mean = np.mean(exam_scores)
numerator = np.sum((study_hours - x_mean) * (exam_scores - y_mean))
denominator = np.sqrt(np.sum((study_hours - x_mean)**2) * np.sum((exam_scores - y_mean)**2))
r_manual = numerator / denominator
print(f"手动计算 r = {r_manual:.4f}")

# 方法2：使用 numpy
r_numpy = np.corrcoef(study_hours, exam_scores)[0, 1]
print(f"numpy 计算 r = {r_numpy:.4f}")

# 方法3：使用 scipy（同时获得p值）
r_scipy, p_value = stats.pearsonr(study_hours, exam_scores)
print(f"scipy 计算 r = {r_scipy:.4f}, p = {p_value:.6f}")

# 相关强度判断
def interpret_correlation(r):
    abs_r = abs(r)
    if abs_r >= 0.8: return "极强相关"
    elif abs_r >= 0.6: return "强相关"
    elif abs_r >= 0.4: return "中等相关"
    elif abs_r >= 0.2: return "弱相关"
    else: return "极弱或无相关"

print(f"\n相关强度: {interpret_correlation(r_scipy)}")
print(f"R² = {r_scipy**2:.4f}（{r_scipy**2*100:.1f}%的方差被解释）")

# 显著性检验
t_stat = r_scipy * np.sqrt((len(study_hours) - 2) / (1 - r_scipy**2))
print(f"t统计量 = {t_stat:.4f}")
print(f"在α=0.05水平下，相关性{'显著' if p_value < 0.05 else '不显著'}")

import numpy as np
from scipy import stats

# 示例1：单调非线性关系
x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
y = np.exp(x * 0.3)  # 指数增长

# 皮尔逊 vs 斯皮尔曼
r_pearson, p_pearson = stats.pearsonr(x, y)
r_spearman, p_spearman = stats.spearmanr(x, y)

print("=== 单调非线性关系：y = e^(0.3x) ===")
print(f"Pearson  r = {r_pearson:.4f}, p = {p_pearson:.6f}")
print(f"Spearman ρ = {r_spearman:.4f}, p = {p_spearman:.6f}")
print("→ Spearman = 1.0（完美单调关系），Pearson < 1")

# 示例2：含异常值的数据
np.random.seed(42)
x2 = np.random.randn(100)
y2 = x2 + np.random.randn(100) * 0.5
# 加入异常值
x2[99] = 10
y2[99] = -10

r2_p, p2_p = stats.pearsonr(x2, y2)
r2_s, p2_s = stats.spearmanr(x2, y2)

print("\n=== 含异常值的数据 ===")
print(f"Pearson  r = {r2_p:.4f}, p = {p2_p:.6f}")
print(f"Spearman ρ = {r2_s:.4f}, p = {p2_s:.6f}")
print("→ Spearman 受异常值影响小，更稳健")

# 示例3：手动计算（演示秩次转换）
x3 = [85, 92, 78, 95, 88]
y3 = [80, 90, 75, 92, 85]

rank_x = [sorted(x3).index(v)+1 for v in x3]
rank_y = [sorted(y3).index(v)+1 for v in y3]
d_squared = sum((rx - ry)**2 for rx, ry in zip(rank_x, rank_y))
n = len(x3)
rho_manual = 1 - 6*d_squared / (n*(n**2 - 1))

print("\n=== 手动计算斯皮尔曼 ===")
print(f"秩次X: {rank_x}")
print(f"秩次Y: {rank_y}")
print(f"Σd² = {d_squared}")
print(f"ρ = 1 - 6×{d_squared}/({n}×{n**2-1}) = {rho_manual:.4f}")

import numpy as np
from scipy import stats

np.random.seed(42)
n = 100

# 场景1：完美线性
x1 = np.linspace(1, 10, n)
y1 = 2 * x1 + np.random.randn(n) * 0.5

# 场景2：非线性单调（指数）
x2 = np.linspace(0.5, 3, n)
y2 = np.exp(x2)

# 场景3：非单调（二次）
x3 = np.linspace(-3, 3, n)
y3 = x3**2 + np.random.randn(n) * 0.3

# 场景4：含异常值
x4 = np.random.randn(n)
y4 = x4 + np.random.randn(n) * 0.5
x4[-1], y4[-1] = 8, -8  # 异常值

def compare_all(x, y, label):
    rp, pp = stats.pearsonr(x, y)
    rs, ps = stats.spearmanr(x, y)
    rk, pk = stats.kendalltau(x, y)
    print(f"\n=== {label} ===")
    print(f"Pearson:  r = {rp:+.4f} (p={pp:.4f})")
    print(f"Spearman: ρ = {rs:+.4f} (p={ps:.4f})")
    print(f"Kendall:  τ = {rk:+.4f} (p={pk:.4f})")

compare_all(x1, y1, "线性+噪声")
compare_all(x2, y2, "指数关系")
compare_all(x3, y3, "二次关系（非单调）")
compare_all(x4, y4, "含异常值")

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d, CubicSpline
from scipy.optimize import curve_fit

# 数据点
x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([0, 0.8, 0.9, 0.1, -0.8, -1.0])

# 1. 拉格朗日插值（高次多项式）
from scipy.interpolate import BarycentricInterpolator
lagrange = BarycentricInterpolator(x_data, y_data)

# 2. 三次样条插值
cs = CubicSpline(x_data, y_data)

# 3. 最小二乘拟合（二次多项式）
coeffs = np.polyfit(x_data, y_data, 2)
poly2 = np.poly1d(coeffs)

# 绘图比较
x_fine = np.linspace(-0.5, 5.5, 500)
plt.figure(figsize=(10, 6))
plt.plot(x_data, y_data, 'ko', markersize=8, label='数据点')
plt.plot(x_fine, lagrange(x_fine), 'b--', linewidth=2, label='拉格朗日插值')
plt.plot(x_fine, cs(x_fine), 'r-', linewidth=2, label='三次样条')
plt.plot(x_fine, poly2(x_fine), 'g:', linewidth=2.5, label='最小二乘(二次)')
plt.legend()
plt.title('插值 vs 拟合')
plt.grid(True, alpha=0.3)
plt.savefig('interp_vs_fit.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 目标函数：f(x,y) = x² + 2y²
def f(x, y): return x**2 + 2*y**2
def grad_f(x, y): return np.array([2*x, 4*y])
def hessian_f(x, y): return np.array([[2, 0], [0, 4]])

# 梯度下降
def gradient_descent(start, lr=0.1, epochs=50):
    x = np.array(start, dtype=float)
    path = [x.copy()]
    for i in range(epochs):
        g = grad_f(x[0], x[1])
        x = x - lr * g
        path.append(x.copy())
        if np.linalg.norm(g) < 1e-6: break
    return x, path

# 牛顿法
def newton_method(start, epochs=20):
    x = np.array(start, dtype=float)
    path = [x.copy()]
    for i in range(epochs):
        g = grad_f(x[0], x[1])
        H = hessian_f(x[0], x[1])
        x = x - np.linalg.solve(H, g)
        path.append(x.copy())
        if np.linalg.norm(g) < 1e-6: break
    return x, path

# 运行比较
start = [3.0, 3.0]
opt_gd, path_gd = gradient_descent(start, 0.1)
opt_newton, path_newton = newton_method(start)

print(f"梯度下降: 最优点 ({opt_gd[0]:.6f}, {opt_gd[1]:.6f}), f={f(*opt_gd):.8f}")
print(f"牛顿法:   最优点 ({opt_newton[0]:.6f}, {opt_newton[1]:.6f}), f={f(*opt_newton):.8f}")
print(f"梯度下降迭代: {len(path_gd)-1}, 牛顿法迭代: {len(path_newton)-1}")

import numpy as np

# 1. 估算 π
def estimate_pi(n=100000):
    x = np.random.uniform(-1, 1, n)
    y = np.random.uniform(-1, 1, n)
    inside = (x**2 + y**2 <= 1)
    return 4 * np.sum(inside) / n

# 2. 生日问题
def birthday_problem(n=23, simulations=10000):
    count = 0
    for _ in range(simulations):
        birthdays = np.random.randint(1, 366, n)
        if len(birthdays) != len(set(birthdays)):
            count += 1
    return count / simulations

# 3. 蒙特卡洛积分：∫₀¹ x²dx = 1/3
def monte_carlo_integral(n=100000):
    x = np.random.uniform(0, 1, n)
    return np.mean(x**2)

# 4. 布丰投针估算 π
def buffon_needle(n=10000, needle_length=0.5, line_spacing=1.0):
    x = np.random.uniform(0, line_spacing/2, n)
    theta = np.random.uniform(0, np.pi/2, n)
    crosses = (x <= needle_length/2 * np.cos(theta))
    p = np.mean(crosses)
    return (2 * needle_length) / (p * line_spacing)

print(f"π 估算: {estimate_pi():.6f}")
print(f"生日问题 (n=23): {birthday_problem():.4f}")
print(f"∫₀¹ x²dx ≈ {monte_carlo_integral():.6f} (理论: 0.3333)")
print(f"布丰投针 π ≈ {buffon_needle():.6f}")