Jacky's blog
Jacky
2025-06-19
目录

scipy

# 基础知识

SciPy是基于NumPy构建的开源科学计算库，提供了大量用于科学计算、工程计算和数据分析的高级函数和算法。示例 (opens new window)

# 安装与导入

# 安装
pip install scipy

# 导入
import numpy as np
import scipy
from scipy import optimize, integrate, interpolate, linalg, stats, signal, ndimage
1
2
3
4
5
6
7

# SciPy与NumPy的关系

  • NumPy：提供基础数组操作和基本数学函数
  • SciPy：基于NumPy，提供高级科学计算算法
  • 依赖关系：SciPy依赖NumPy，但功能更专业化

# 主要子模块概览

# 核心子模块

  • scipy.optimize：优化算法与拟合
  • scipy.integrate：数值积分与微分方程
  • scipy.interpolate：插值与拟合
  • scipy.linalg：线性代数运算
  • scipy.stats：统计函数与概率分布
  • scipy.signal：信号处理
  • scipy.ndimage：图像处理
  • scipy.sparse：稀疏矩阵
  • scipy.spatial：空间数据结构与算法
  • scipy.constants：物理常数
  • scipy.special：特殊函数

# 优化与拟合 (scipy.optimize)

# 函数最小化

from scipy.optimize import minimize, minimize_scalar

# 单变量最小化
def f(x):
    return (x - 2)**2 + 1

result = minimize_scalar(f)
print(f"最小值: {result.x:.4f}, 函数值: {result.fun:.4f}")

# 多变量最小化
def rosen(x):
    return (1 - x[0])**2 + 100 * (x[1] - x[0]**2)**2

x0 = [-1, 1]
result = minimize(rosen, x0, method='BFGS')
print(f"最优解: {result.x}")
print(f"最优值: {result.fun}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 曲线拟合

from scipy.optimize import curve_fit

# 定义拟合函数
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

# 生成数据
xdata = np.linspace(0, 4, 50)
ydata = func(xdata, 2.5, 1.3, 0.5) + 0.2 * np.random.normal(size=len(xdata))

# 拟合
popt, pcov = curve_fit(func, xdata, ydata)
print(f"拟合参数: a={popt[0]:.3f}, b={popt[1]:.3f}, c={popt[2]:.3f}")

# 计算拟合优度
residuals = ydata - func(xdata, *popt)
ss_res = np.sum(residuals**2)
ss_tot = np.sum((ydata - np.mean(ydata))**2)
r_squared = 1 - (ss_res / ss_tot)
print(f"R² = {r_squared:.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# 线性规划

线性规划是一种数学优化方法，用于在满足一组线性约束条件的情况下，找到线性目标函数的最大值或最小值。

from scipy.optimize import linprog

# 最小化: -x1 - 2x2
# 约束: x1 + x2 <= 3, x1 >= 0, x2 >= 0
c = [-1, -2]  # 目标函数系数
A = [[1, 1]]  # 约束矩阵
b = [3]       # 约束右端
bounds = [(0, None), (0, None)]  # 变量边界

result = linprog(c, A_ub=A, b_ub=b, bounds=bounds)
print(f"最优解: {result.x}")
print(f"最优值: {result.fun}")
1
2
3
4
5
6
7
8
9
10
11
12

# 数值积分与微分方程 (scipy.integrate)

# 数值积分

from scipy.integrate import quad, dblquad, nquad

# 一维定积分
def integrand(x):
    return x**2

result, error = quad(integrand, 0, 1)
print(f"积分结果: {result:.4f}, 误差: {error:.2e}")

# 二重积分
def integrand_2d(x, y):
    return x**2 + y**2

result_2d, error_2d = dblquad(integrand_2d, 0, 1, lambda x: 0, lambda x: 1)
print(f"二重积分: {result_2d:.4f}")

# 带参数的积分
def integrand_param(x, a, b):
    return a * x**2 + b * x

result_param, _ = quad(integrand_param, 0, 1, args=(2, 3))
print(f"参数积分: {result_param:.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

# 常微分方程求解

from scipy.integrate import odeint, solve_ivp

# 简单ODE: dy/dt = -2y
def model(y, t):
    return -2 * y

t = np.linspace(0, 5, 100)
y0 = 1.0
solution = odeint(model, y0, t)

# 使用solve_ivp（更现代的方法）
def model_ivp(t, y):
    return -2 * y

sol = solve_ivp(model_ivp, [0, 5], [1.0], t_eval=t)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 偏微分方程

# 热传导方程示例
def heat_equation_solver(nx=50, nt=100, dx=0.1, dt=0.001, alpha=1.0):
    """一维热传导方程求解"""
    x = np.linspace(0, nx*dx, nx)
    t = np.linspace(0, nt*dt, nt)
    
    # 初始条件
    u = np.zeros((nt, nx))
    u[0, :] = np.sin(np.pi * x / (nx * dx))
    
    # 边界条件
    u[:, 0] = 0
    u[:, -1] = 0
    
    # 时间步进
    for n in range(0, nt-1):
        for i in range(1, nx-1):
            u[n+1, i] = u[n, i] + alpha * dt / dx**2 * (u[n, i+1] - 2*u[n, i] + u[n, i-1])
    
    return x, t, u

x, t, u = heat_equation_solver()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

# 插值 (scipy.interpolate)

# 一维插值

from scipy.interpolate import interp1d, CubicSpline, lagrange

# 数据点
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

# 线性插值
linear_interp = interp1d(x, y, kind='linear')
x_new = np.linspace(0, 4, 100)
y_linear = linear_interp(x_new)

# 三次样条插值
cs = CubicSpline(x, y)
y_spline = cs(x_new)

# 拉格朗日插值
poly = lagrange(x, y)
y_lagrange = poly(x_new)

# 比较不同插值方法
plt.figure(figsize=(12, 8))
plt.plot(x, y, 'o', label='数据点', markersize=8)
plt.plot(x_new, y_linear, label='线性插值', alpha=0.7)
plt.plot(x_new, y_spline, label='三次样条', alpha=0.7)
plt.plot(x_new, y_lagrange, label='拉格朗日', alpha=0.7)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

# 多维插值

from scipy.interpolate import griddata, RBFInterpolator

# 二维散点数据
np.random.seed(42)
x = np.random.rand(100)
y = np.random.rand(100)
z = np.sin(10*x) + np.cos(10*y)

# 规则网格
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:100j]

# 网格插值
grid_z = griddata((x, y), z, (grid_x, grid_y), method='cubic')

# RBF插值
rbf = RBFInterpolator(np.column_stack([x, y]), z, kernel='thin_plate_spline')
grid_points = np.column_stack([grid_x.ravel(), grid_y.ravel()])
grid_z_rbf = rbf(grid_points).reshape(grid_x.shape)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 线性代数 (scipy.linalg)

# 基本矩阵运算

from scipy import linalg

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

# 矩阵求逆
A_inv = linalg.inv(A)
print(f"逆矩阵:\n{A_inv}")

# 行列式
det_A = linalg.det(A)
print(f"行列式: {det_A}")

# 矩阵乘法
C = linalg.matmul(A, B)
print(f"矩阵乘法:\n{C}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 特征值与特征向量

# 特征值分解
eigenvals, eigenvecs = linalg.eig(A)
print(f"特征值: {eigenvals}")
print(f"特征向量:\n{eigenvecs}")

# 实对称矩阵的特征值分解
A_sym = np.array([[2, 1], [1, 2]])
eigenvals_sym, eigenvecs_sym = linalg.eigh(A_sym)
print(f"对称矩阵特征值: {eigenvals_sym}")
1
2
3
4
5
6
7
8
9

# 奇异值分解

# SVD分解
U, s, Vh = linalg.svd(A)
print(f"U矩阵:\n{U}")
print(f"奇异值: {s}")
print(f"V^T矩阵:\n{Vh}")

# 重构矩阵
A_reconstructed = U @ np.diag(s) @ Vh
print(f"重构误差: {np.linalg.norm(A - A_reconstructed)}")
1
2
3
4
5
6
7
8
9

# 线性方程组求解

# 线性方程组 Ax = b
b = np.array([5, 6])
x = linalg.solve(A, b)
print(f"解: {x}")

# 最小二乘解
A_over = np.array([[1, 2], [3, 4], [5, 6]])
b_over = np.array([1, 2, 3])
x_lstsq = linalg.lstsq(A_over, b_over)[0]
print(f"最小二乘解: {x_lstsq}")
1
2
3
4
5
6
7
8
9
10

# 统计与概率分布 (scipy.stats)

# 描述性统计

from scipy import stats

# 生成数据
data = np.random.normal(0, 1, 1000)

# 基本统计量
print(f"样本数: {stats.describe(data).nobs}")
print(f"均值: {stats.describe(data).mean:.4f}")
print(f"方差: {stats.describe(data).variance:.4f}")
print(f"偏度: {stats.describe(data).skewness:.4f}")
print(f"峰度: {stats.describe(data).kurtosis:.4f}")

# 分位数
print(f"中位数: {np.median(data):.4f}")
print(f"25%分位数: {np.percentile(data, 25):.4f}")
print(f"75%分位数: {np.percentile(data, 75):.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 概率分布

# 正态分布
normal_dist = stats.norm(loc=0, scale=1)
x = np.linspace(-4, 4, 100)
pdf = normal_dist.pdf(x)
cdf = normal_dist.cdf(x)

# 其他常用分布
exponential_dist = stats.expon(scale=1)
gamma_dist = stats.gamma(a=2, scale=1)
beta_dist = stats.beta(a=2, b=5)

# 分布拟合
params = stats.norm.fit(data)
print(f"拟合参数: 均值={params[0]:.4f}, 标准差={params[1]:.4f}")

# 分布检验
statistic, p_value = stats.kstest(data, 'norm', params)
print(f"KS检验: 统计量={statistic:.4f}, p值={p_value:.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 假设检验

# t检验
t_stat, p_value = stats.ttest_1samp(data, 0)
print(f"单样本t检验: t={t_stat:.4f}, p={p_value:.4f}")

# 双样本t检验
data2 = np.random.normal(0.5, 1, 1000)
t_stat, p_value = stats.ttest_ind(data, data2)
print(f"双样本t检验: t={t_stat:.4f}, p={p_value:.4f}")

# 方差分析
group1 = np.random.normal(0, 1, 100)
group2 = np.random.normal(1, 1, 100)
group3 = np.random.normal(2, 1, 100)
f_stat, p_value = stats.f_oneway(group1, group2, group3)
print(f"ANOVA: F={f_stat:.4f}, p={p_value:.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 相关性分析

# 皮尔逊相关系数
x = np.random.randn(100)
y = 0.5 * x + np.random.randn(100) * 0.5
correlation, p_value = stats.pearsonr(x, y)
print(f"皮尔逊相关系数: r={correlation:.4f}, p={p_value:.4f}")

# 斯皮尔曼相关系数
correlation, p_value = stats.spearmanr(x, y)
print(f"斯皮尔曼相关系数: r={correlation:.4f}, p={p_value:.4f}")

# 肯德尔相关系数
correlation, p_value = stats.kendalltau(x, y)
print(f"肯德尔相关系数: τ={correlation:.4f}, p={p_value:.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13

# 信号处理 (scipy.signal)

# 滤波

from scipy import signal

# 生成信号
t = np.linspace(0, 1, 1000)
signal_data = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.1 * np.random.randn(1000)

# 低通滤波
b, a = signal.butter(4, 0.1, 'low')
filtered_low = signal.filtfilt(b, a, signal_data)

# 高通滤波
b, a = signal.butter(4, 0.1, 'high')
filtered_high = signal.filtfilt(b, a, signal_data)

# 带通滤波
b, a = signal.butter(4, [0.05, 0.15], 'band')
filtered_band = signal.filtfilt(b, a, signal_data)

# 可视化
plt.figure(figsize=(15, 10))
plt.subplot(4, 1, 1)
plt.plot(t, signal_data)
plt.title('原始信号')
plt.subplot(4, 1, 2)
plt.plot(t, filtered_low)
plt.title('低通滤波')
plt.subplot(4, 1, 3)
plt.plot(t, filtered_high)
plt.title('高通滤波')
plt.subplot(4, 1, 4)
plt.plot(t, filtered_band)
plt.title('带通滤波')
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

# 频谱分析

# 功率谱密度
freqs, psd = signal.welch(signal_data, fs=1000)
plt.figure(figsize=(12, 4))
plt.semilogy(freqs, psd)
plt.xlabel('频率 (Hz)')
plt.ylabel('功率谱密度')
plt.title('功率谱密度')
plt.grid(True, alpha=0.3)
plt.show()

# 短时傅里叶变换
f, t, Sxx = signal.spectrogram(signal_data, fs=1000, nperseg=256)
plt.figure(figsize=(12, 4))
plt.pcolormesh(t, f, 10 * np.log10(Sxx))
plt.ylabel('频率 (Hz)')
plt.xlabel('时间 (s)')
plt.title('频谱图')
plt.colorbar(label='功率 (dB)')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# 窗函数

# 常用窗函数
N = 100
hann_window = signal.hanning(N)
hamming_window = signal.hamming(N)
blackman_window = signal.blackman(N)
kaiser_window = signal.kaiser(N, beta=14)

plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.plot(hann_window)
plt.title('Hanning窗')
plt.subplot(2, 2, 2)
plt.plot(hamming_window)
plt.title('Hamming窗')
plt.subplot(2, 2, 3)
plt.plot(blackman_window)
plt.title('Blackman窗')
plt.subplot(2, 2, 4)
plt.plot(kaiser_window)
plt.title('Kaiser窗')
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

# 图像处理 (scipy.ndimage)

# 图像滤波

from scipy import ndimage

# 创建测试图像
image = np.random.rand(100, 100)

# 高斯滤波
blurred = ndimage.gaussian_filter(image, sigma=2)

# 中值滤波
median_filtered = ndimage.median_filter(image, size=5)

# 维纳滤波
wiener_filtered = ndimage.wiener(image, (5, 5))

# 可视化
plt.figure(figsize=(15, 10))
plt.subplot(2, 2, 1)
plt.imshow(image, cmap='gray')
plt.title('原始图像')
plt.subplot(2, 2, 2)
plt.imshow(blurred, cmap='gray')
plt.title('高斯滤波')
plt.subplot(2, 2, 3)
plt.imshow(median_filtered, cmap='gray')
plt.title('中值滤波')
plt.subplot(2, 2, 4)
plt.imshow(wiener_filtered, cmap='gray')
plt.title('维纳滤波')
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

# 图像变换

# 旋转
rotated = ndimage.rotate(image, angle=45)

# 缩放
scaled = ndimage.zoom(image, zoom=2.0)

# 平移
shifted = ndimage.shift(image, shift=[10, 10])

# 仿射变换
from scipy.ndimage import affine_transform
matrix = np.array([[np.cos(np.pi/4), -np.sin(np.pi/4)],
                   [np.sin(np.pi/4), np.cos(np.pi/4)]])
affine_transformed = affine_transform(image, matrix)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 形态学操作

from scipy.ndimage import binary_erosion, binary_dilation, binary_opening, binary_closing

# 二值化
binary_image = image > 0.5

# 腐蚀
eroded = binary_erosion(binary_image)

# 膨胀
dilated = binary_dilation(binary_image)

# 开运算（先腐蚀后膨胀）
opened = binary_opening(binary_image)

# 闭运算（先膨胀后腐蚀）
closed = binary_closing(binary_image)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 稀疏矩阵 (scipy.sparse)

# 稀疏矩阵类型

from scipy import sparse

# 创建稀疏矩阵
data = np.array([1, 2, 3, 4, 5, 6])
row = np.array([0, 0, 1, 1, 2, 2])
col = np.array([0, 1, 1, 2, 0, 2])

# COO格式
coo_matrix = sparse.coo_matrix((data, (row, col)), shape=(3, 3))

# CSR格式（行压缩）
csr_matrix = sparse.csr_matrix((data, (row, col)), shape=(3, 3))

# CSC格式（列压缩）
csc_matrix = sparse.csc_matrix((data, (row, col)), shape=(3, 3))

# 转换为密集矩阵
dense_matrix = csr_matrix.toarray()
print(f"稀疏矩阵:\n{dense_matrix}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# 稀疏矩阵运算

# 矩阵乘法
result = csr_matrix.dot(csr_matrix)
print(f"矩阵乘法结果:\n{result.toarray()}")

# 特征值
eigenvals = sparse.linalg.eigs(csr_matrix, k=2)
print(f"特征值: {eigenvals[0]}")

# 线性方程组求解
b = np.array([1, 2, 3])
x = sparse.linalg.spsolve(csr_matrix, b)
print(f"解: {x}")
1
2
3
4
5
6
7
8
9
10
11
12

# 空间数据结构 (scipy.spatial)

# 距离计算

from scipy.spatial import distance, distance_matrix

# 点集
points = np.random.rand(10, 2)

# 欧几里得距离
dist_matrix = distance_matrix(points, points)
print(f"距离矩阵:\n{dist_matrix}")

# 其他距离度量
manhattan_dist = distance.cityblock(points[0], points[1])
cosine_dist = distance.cosine(points[0], points[1])
print(f"曼哈顿距离: {manhattan_dist:.4f}")
print(f"余弦距离: {cosine_dist:.4f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# KD树最近邻搜索

from scipy.spatial import KDTree

# 创建KD树
tree = KDTree(points)

# 最近邻搜索
query_point = np.array([0.5, 0.5])
distance, index = tree.query(query_point)
print(f"最近点索引: {index}, 距离: {distance:.4f}")

# k近邻搜索
distances, indices = tree.query(query_point, k=3)
print(f"3个最近点: {indices}")
print(f"对应距离: {distances}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 凸包与三角剖分

from scipy.spatial import ConvexHull, Delaunay

# 凸包
hull = ConvexHull(points)
print(f"凸包顶点: {hull.vertices}")

# 三角剖分
tri = Delaunay(points)
print(f"三角形数量: {len(tri.simplices)}")
1
2
3
4
5
6
7
8
9

# 物理常数与特殊函数

# 物理常数

from scipy import constants

# 基本物理常数
print(f"光速: {constants.c} m/s")
print(f"普朗克常数: {constants.h} J·s")
print(f"玻尔兹曼常数: {constants.k} J/K")
print(f"电子质量: {constants.m_e} kg")
print(f"质子质量: {constants.m_p} kg")
print(f"电子电荷: {constants.e} C")
print(f"真空介电常数: {constants.epsilon_0} F/m")
print(f"真空磁导率: {constants.mu_0} H/m")

# 单位转换
print(f"1 eV = {constants.eV} J")
print(f"1 Å = {constants.angstrom} m")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 特殊函数

from scipy import special

# 贝塞尔函数
x = np.linspace(0, 10, 100)
bessel_j0 = special.j0(x)
bessel_j1 = special.j1(x)
bessel_y0 = special.y0(x)
bessel_y1 = special.y1(x)

# 伽马函数
gamma_values = special.gamma(x)
log_gamma = special.loggamma(x)

# 误差函数
error_function = special.erf(x)
complementary_error = special.erfc(x)

# 椭圆函数
elliptic_k = special.ellipk(0.5)
elliptic_e = special.ellipe(0.5)

# 可视化
plt.figure(figsize=(15, 10))
plt.subplot(2, 3, 1)
plt.plot(x, bessel_j0, label='J0')
plt.plot(x, bessel_j1, label='J1')
plt.title('贝塞尔函数')
plt.legend()
plt.subplot(2, 3, 2)
plt.plot(x, gamma_values)
plt.title('伽马函数')
plt.subplot(2, 3, 3)
plt.plot(x, error_function, label='erf')
plt.plot(x, complementary_error, label='erfc')
plt.title('误差函数')
plt.legend()
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

# 实际应用案例

# 1. 物理模拟：弹簧-质量系统

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

def spring_mass_system(state, t, k, m, c):
    """弹簧-质量-阻尼系统"""
    x, v = state
    dx_dt = v
    dv_dt = (-k/m) * x - (c/m) * v
    return [dx_dt, dv_dt]

# 参数
k = 1.0  # 弹簧常数
m = 1.0  # 质量
c = 0.1  # 阻尼系数

# 时间序列
t = np.linspace(0, 20, 1000)
initial_state = [1.0, 0.0]  # 初始位置和速度

# 求解ODE
solution = odeint(spring_mass_system, initial_state, t, args=(k, m, c))

# 可视化
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(t, solution[:, 0], 'b-', linewidth=2)
plt.ylabel('位置')
plt.title('弹簧-质量-阻尼系统')
plt.grid(True, alpha=0.3)
plt.subplot(2, 1, 2)
plt.plot(t, solution[:, 1], 'r-', linewidth=2)
plt.xlabel('时间')
plt.ylabel('速度')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

# 2. 数据分析：正态分布拟合与检验

from scipy import stats
import matplotlib.pyplot as plt

# 生成混合正态分布数据
np.random.seed(42)
data1 = np.random.normal(0, 1, 500)
data2 = np.random.normal(3, 1.5, 300)
data = np.concatenate([data1, data2])

# 描述性统计
print("=== 描述性统计 ===")
print(f"样本数: {len(data)}")
print(f"均值: {np.mean(data):.4f}")
print(f"标准差: {np.std(data):.4f}")
print(f"偏度: {stats.skew(data):.4f}")
print(f"峰度: {stats.kurtosis(data):.4f}")

# 分布拟合
params = stats.norm.fit(data)
print(f"\n=== 正态分布拟合 ===")
print(f"拟合均值: {params[0]:.4f}")
print(f"拟合标准差: {params[1]:.4f}")

# 分布检验
statistic, p_value = stats.kstest(data, 'norm', params)
print(f"\n=== KS检验 ===")
print(f"统计量: {statistic:.4f}")
print(f"p值: {p_value:.4f}")

# 可视化
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.hist(data, bins=30, density=True, alpha=0.7, color='skyblue', edgecolor='black')
x = np.linspace(data.min(), data.max(), 100)
plt.plot(x, stats.norm.pdf(x, *params), 'r-', linewidth=2, label='拟合正态分布')
plt.title('数据分布与拟合')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 2)
stats.probplot(data, dist="norm", plot=plt)
plt.title('Q-Q图')

plt.subplot(2, 2, 3)
plt.hist(data, bins=30, cumulative=True, density=True, alpha=0.7, color='lightgreen')
plt.plot(x, stats.norm.cdf(x, *params), 'r-', linewidth=2)
plt.title('累积分布函数')
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 4)
plt.boxplot(data)
plt.title('箱线图')
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

# 3. 信号处理：音频信号分析

from scipy import signal
import matplotlib.pyplot as plt

# 生成模拟音频信号
fs = 44100  # 采样频率
duration = 1.0  # 持续时间
t = np.linspace(0, duration, int(fs * duration))

# 复合信号：基频 + 谐波 + 噪声
f0 = 440  # 基频 (A4音符)
signal_clean = np.sin(2 * np.pi * f0 * t) + 0.3 * np.sin(2 * np.pi * 2 * f0 * t)
signal_noisy = signal_clean + 0.1 * np.random.randn(len(t))

# 频谱分析
freqs, psd = signal.welch(signal_noisy, fs, nperseg=1024)

# 滤波器设计
b, a = signal.butter(4, 0.1, 'low')
signal_filtered = signal.filtfilt(b, a, signal_noisy)

# 短时傅里叶变换
f, t_stft, Sxx = signal.spectrogram(signal_noisy, fs, nperseg=512, noverlap=256)

# 可视化
plt.figure(figsize=(15, 10))
plt.subplot(3, 2, 1)
plt.plot(t[:1000], signal_clean[:1000])
plt.title('原始信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')

plt.subplot(3, 2, 2)
plt.plot(t[:1000], signal_noisy[:1000])
plt.title('含噪声信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')

plt.subplot(3, 2, 3)
plt.plot(t[:1000], signal_filtered[:1000])
plt.title('滤波后信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')

plt.subplot(3, 2, 4)
plt.semilogy(freqs, psd)
plt.title('功率谱密度')
plt.xlabel('频率 (Hz)')
plt.ylabel('功率谱密度')
plt.grid(True, alpha=0.3)

plt.subplot(3, 2, 5)
plt.pcolormesh(t_stft, f, 10 * np.log10(Sxx))
plt.title('频谱图')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.colorbar(label='功率 (dB)')

plt.subplot(3, 2, 6)
plt.plot(freqs, psd)
plt.title('频谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.xlim(0, 2000)
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

# 4. 图像处理：医学图像分析

from scipy import ndimage
import matplotlib.pyplot as plt

# 创建模拟医学图像
np.random.seed(42)
image = np.random.rand(200, 200)

# 添加一些结构（模拟器官）
x, y = np.meshgrid(np.linspace(0, 1, 200), np.linspace(0, 1, 200))
center1 = np.exp(-((x-0.3)**2 + (y-0.3)**2) / 0.02)
center2 = np.exp(-((x-0.7)**2 + (y-0.7)**2) / 0.03)
image = image + 0.5 * center1 + 0.3 * center2

# 添加噪声
image_noisy = image + 0.1 * np.random.randn(200, 200)

# 图像处理
# 1. 高斯滤波去噪
image_smoothed = ndimage.gaussian_filter(image_noisy, sigma=1.0)

# 2. 边缘检测
edges_sobel = ndimage.sobel(image_smoothed)
edges_laplace = ndimage.laplace(image_smoothed)

# 3. 形态学操作
binary = image_smoothed > np.mean(image_smoothed)
eroded = ndimage.binary_erosion(binary)
dilated = ndimage.binary_dilation(binary)
opened = ndimage.binary_opening(binary)
closed = ndimage.binary_closing(binary)

# 4. 连通组件分析
labeled, num_features = ndimage.label(binary)
print(f"检测到 {num_features} 个连通区域")

# 可视化
plt.figure(figsize=(15, 12))
plt.subplot(3, 3, 1)
plt.imshow(image, cmap='gray')
plt.title('原始图像')
plt.axis('off')

plt.subplot(3, 3, 2)
plt.imshow(image_noisy, cmap='gray')
plt.title('含噪声图像')
plt.axis('off')

plt.subplot(3, 3, 3)
plt.imshow(image_smoothed, cmap='gray')
plt.title('平滑后图像')
plt.axis('off')

plt.subplot(3, 3, 4)
plt.imshow(edges_sobel, cmap='gray')
plt.title('Sobel边缘')
plt.axis('off')

plt.subplot(3, 3, 5)
plt.imshow(edges_laplace, cmap='gray')
plt.title('Laplace边缘')
plt.axis('off')

plt.subplot(3, 3, 6)
plt.imshow(binary, cmap='gray')
plt.title('二值化')
plt.axis('off')

plt.subplot(3, 3, 7)
plt.imshow(eroded, cmap='gray')
plt.title('腐蚀')
plt.axis('off')

plt.subplot(3, 3, 8)
plt.imshow(dilated, cmap='gray')
plt.title('膨胀')
plt.axis('off')

plt.subplot(3, 3, 9)
plt.imshow(labeled, cmap='tab10')
plt.title('连通组件')
plt.axis('off')

plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

# 性能优化建议

# 1. 算法选择

# 选择合适的算法
from scipy import optimize

# 对于凸优化问题，使用更高效的算法
def convex_function(x):
    return x**2 + 2*x + 1

# 使用BFGS算法（适合凸函数）
result_bfgs = optimize.minimize(convex_function, x0=0, method='BFGS')

# 使用L-BFGS-B算法（适合大规模问题）
result_lbfgs = optimize.minimize(convex_function, x0=0, method='L-BFGS-B')
1
2
3
4
5
6
7
8
9
10
11
12

# 2. 内存优化

# 使用稀疏矩阵处理大规模数据
from scipy import sparse

# 创建大型稀疏矩阵
n = 10000
sparse_matrix = sparse.random(n, n, density=0.01, format='csr')

# 稀疏矩阵运算比密集矩阵更高效
result = sparse_matrix.dot(sparse_matrix)
1
2
3
4
5
6
7
8
9

# 3. 并行计算

# 使用多进程进行并行计算
from multiprocessing import Pool
from scipy import integrate

def parallel_integration(func, limits_list):
    """并行积分计算"""
    with Pool() as pool:
        results = pool.starmap(integrate.quad, [(func, a, b) for a, b in limits_list])
    return results
1
2
3
4
5
6
7
8
9

# 常见问题与解决方案

# 1. 安装问题

# 推荐使用conda安装
conda install scipy

# 或者使用pip
pip install scipy

# 如果遇到编译问题，可以安装预编译版本
pip install scipy --only-binary=scipy
1
2
3
4
5
6
7
8

# 2. 性能问题

# 避免不必要的循环
# 错误的方式
def slow_function(data):
    result = []
    for i in range(len(data)):
        result.append(data[i] * 2)
    return result

# 正确的方式
def fast_function(data):
    return data * 2
1
2
3
4
5
6
7
8
9
10
11

# 3. 精度问题

# 处理浮点数精度问题
from scipy import constants

# 使用相对误差比较
def is_close(a, b, rtol=1e-5, atol=1e-8):
    return abs(a - b) <= (atol + rtol * abs(b))

# 在数值积分中指定精度
result, error = integrate.quad(lambda x: x**2, 0, 1, epsabs=1e-10, epsrel=1e-10)
1
2
3
4
5
6
7
8
9

# 扩展资源

# 1. 官方文档

# 2. 学习资源

# 3. 相关库

  • NumPy: 基础数组操作
  • Pandas: 数据分析
  • Matplotlib: 数据可视化
  • Scikit-learn: 机器学习
  • SymPy: 符号计算

# 4. 社区资源

上次更新: 2025/07/21, 21:50:24
matplotlib
python typing 模块

← python typing 模块→

最近更新
01
npx 使用指南
10-12
02
cursor
09-28
03
inspect
07-20
更多文章>
Theme by Vdoing | Copyright © 2019-2025 Jacky | MIT License