1. Scipy简介与安装
Scipy（Scientific Python）是一个开源的Python算法库和数学工具包，它建立在Numpy的基础上，提供了许多高级的数学函数、信号处理工具、优化算法、统计函数等。Scipy是Python科学计算生态系统中的核心组件之一，与Numpy、Matplotlib、Pandas等库共同构成了强大的科学计算环境。
1.1 Scipy的特点
高效的数值计算能力
丰富的科学计算功能模块
与Numpy数组无缝集成
开源免费，社区活跃
跨平台支持
1.2 安装Scipy
在开始使用Scipy之前，我们需要先安装它。安装Scipy最简单的方法是使用pip包管理器：
pip install scipy
如果你使用的是Anaconda发行版，可以使用conda进行安装：
conda install scipy
1.3 验证安装
安装完成后，可以通过以下代码验证Scipy是否安装成功：
import scipy
print(scipy.__version__)
如果成功安装，上述代码将输出你安装的Scipy版本号。
2. Scipy核心模块概述
Scipy由多个子模块组成，每个模块专注于特定的科学计算领域。以下是Scipy的主要模块及其功能：
scipy.cluster
: 向量量化/K-means
scipy.constants
: 物理和数学常数
scipy.fftpack
: 傅里叶变换
scipy.integrate
: 积分和常微分方程求解
scipy.interpolate
: 插值
scipy.io
: 数据输入和输出
scipy.linalg
: 线性代数
scipy.ndimage
: n维图像处理
scipy.odr
: 正交距离回归
scipy.optimize
: 优化和求根
scipy.signal
: 信号处理
scipy.sparse
: 稀疏矩阵
scipy.spatial
: 空间数据结构和算法
scipy.special
: 特殊函数
scipy.stats
: 统计分布和函数
在接下来的章节中，我们将详细探讨其中一些最重要的模块。
3. 数值优化
数值优化是Scipy中最强大的功能之一，
scipy.optimize
模块提供了多种优化算法，可以用于寻找函数的最小值或最大值、求解非线性方程等。
3.1 无约束优化
3.1.1 标量函数最小化
对于单变量函数，我们可以使用
minimize_scalar
函数来寻找最小值。
import numpy as np
from scipy.optimize import minimize_scalar
# 定义一个单变量函数
def f(x):
return x**2 + 10*np.sin(x)
# 寻找函数的最小值
result = minimize_scalar(f)
print(“最小值点:”, result.x)
print(“最小值:”, result.fun)
输出结果：
最小值点: -1.306440991491866
最小值: -7.945823375615215
3.1.2 多变量函数最小化
对于多变量函数，我们可以使用
minimize
函数。
import numpy as np
from scipy.optimize import minimize
# 定义一个多变量函数
def rosen(x):
“””Rosenbrock函数”””
return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
# 初始点
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
# 使用BFGS算法寻找最小值
result = minimize(rosen, x0, method=’BFGS’)
print(“最小值点:”, result.x)
print(“最小值:”, result.fun)
输出结果：
最小值点: [1. 1. 1. 1. 1.]
最小值: 2.156663455538717e-11
3.2 约束优化
在实际问题中，我们常常需要考虑约束条件。Scipy的
minimize
函数支持多种约束优化方法。
import numpy as np
from scipy.optimize import minimize
# 定义目标函数
def objective(x):
return x[0]**2 + x[1]**2
# 定义约束条件
def constraint1(x):
return x[0] + x[1] – 1
# 初始点
x0 = np.array([0.5, 0.5])
# 定义约束
con1 = {‘type’: ‘eq’, ‘fun’: constraint1}
# 求解约束优化问题
solution = minimize(objective, x0, method=’SLSQP’, constraints=con1)
print(“最优解:”, solution.x)
print(“最优值:”, solution.fun)
输出结果：
最优解: [0.5 0.5]
最优值: 0.5
3.3 非线性方程求解
Scipy提供了多种方法来求解非线性方程组。
import numpy as np
from scipy.optimize import fsolve
# 定义方程组
def equations(p):
x, y = p
return (x**2 + y**2 – 1, np.exp(x) + y – 2)
# 初始猜测
x0, y0 = 0, 1
# 求解方程组
x, y = fsolve(equations, (x0, y0))
print(“解为: x =”, x, “, y =”, y)
输出结果：
解为: x = 0.865471033792519 , y = 0.500955377235357
4. 统计分析
scipy.stats
模块提供了大量的统计函数和概率分布，可以用于统计分析、假设检验、概率分布拟合等。
4.1 概率分布
Scipy支持多种连续和离散概率分布。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# 正态分布
mu, sigma = 0, 1 # 均值和标准差
x = np.linspace(mu – 3*sigma, mu + 3*sigma, 100)
y = stats.norm.pdf(x, mu, sigma)
plt.plot(x, y)
plt.title(‘正态分布概率密度函数’)
plt.xlabel(‘x’)
plt.ylabel(‘PDF’)
plt.grid(True)
plt.show()
4.2 统计检验
4.2.1 t检验
t检验用于比较两个样本的均值是否有显著差异。
import numpy as np
from scipy import stats
# 生成两组随机数据
np.random.seed(42)
group1 = np.random.normal(0, 1, 100)
group2 = np.random.normal(0.5, 1, 100)
# 执行独立样本t检验
t_stat, p_value = stats.ttest_ind(group1, group2)
print(“t统计量:”, t_stat)
print(“p值:”, p_value)
# 判断是否显著
alpha = 0.05
if p_value < alpha:
print("两组样本均值有显著差异")
else:
print("两组样本均值无显著差异")
输出结果：
t统计量: -3.799316502273649
p值: 0.00019426361648821436
两组样本均值有显著差异
4.2.2 卡方检验
卡方检验用于检验分类变量之间的独立性。
import numpy as np
from scipy import stats
# 创建一个列联表
observed = np.array([[10, 20, 30], [20, 15, 25]])
# 执行卡方检验
chi2, p_value, dof, expected = stats.chi2_contingency(observed)
print("卡方统计量:", chi2)
print("p值:", p_value)
print("自由度:", dof)
print("期望频数:")
print(expected)
# 判断是否独立
alpha = 0.05
if p_value < alpha:
print("两个分类变量不独立")
else:
print("两个分类变量独立")
输出结果：
卡方统计量: 2.9166666666666665
p值: 0.2326447958569924
自由度: 2
期望频数:
[[15. 17.5 27.5]
[15. 17.5 27.5]]
两个分类变量独立
4.3 线性回归
Scipy提供了多种线性回归方法，包括简单线性回归和多元线性回归。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# 生成一些随机数据
np.random.seed(42)
x = np.random.random(20) * 10
y = 2 * x + 1 + np.random.normal(0, 1, 20)
# 执行线性回归
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
print("斜率:", slope)
print("截距:", intercept)
print("R平方值:", r_value**2)
print("p值:", p_value)
print("标准误差:", std_err)
# 绘制数据和回归线
plt.scatter(x, y, label='数据点')
plt.plot(x, intercept + slope*x, 'r', label='回归线')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('线性回归')
plt.grid(True)
plt.show()
输出结果：
斜率: 1.9889265303147984
截距: 1.0996168158555173
R平方值: 0.966965415620522
p值: 1.323721436545421e-13
标准误差: 0.0896125229197226
5. 信号处理
scipy.signal
模块提供了多种信号处理工具，包括滤波、频谱分析、小波变换等。
5.1 滤波
滤波是信号处理中的基本操作，用于去除噪声或提取特定频率成分。
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# 创建一个信号
t = np.linspace(0, 1, 1000, False) # 1秒
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) # 10Hz和20Hz的正弦波
sig_noisy = sig + np.random.normal(0, 0.5, len(t)) # 添加噪声
# 设计一个Butterworth低通滤波器
b, a = signal.butter(4, 15, 'low', fs=1000) # 截止频率15Hz
filtered = signal.filtfilt(b, a, sig_noisy)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t, sig_noisy, label='带噪声信号')
plt.plot(t, filtered, label='滤波后信号')
plt.plot(t, sig, 'k--', label='原始信号')
plt.xlabel('时间 [s]')
plt.ylabel('幅值')
plt.legend()
plt.grid(True)
plt.title('Butterworth低通滤波')
plt.show()
5.2 频谱分析
频谱分析可以帮助我们了解信号的频率成分。
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# 创建一个信号
fs = 1000 # 采样频率
t = np.arange(0, 1, 1/fs) # 1秒
f1, f2 = 50, 120 # 频率
sig = np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t)
# 添加噪声
noise = 2.5 * np.random.randn(len(t))
sig_noisy = sig + noise
# 计算功率谱密度
f, Pxx = signal.welch(sig_noisy, fs, nperseg=1024)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.semilogy(f, Pxx)
plt.xlabel('频率 [Hz]')
plt.ylabel('功率谱密度')
plt.grid(True)
plt.title('功率谱密度')
plt.show()
5.3 小波变换
小波变换是一种时频分析方法，适用于非平稳信号的分析。
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# 创建一个信号
t = np.linspace(-1, 1, 200, endpoint=False)
sig = np.cos(2 * np.pi * 7 * t) * signal.gausspulse(t, fc=2)
# 计算连续小波变换
widths = np.arange(1, 31)
cwtmatr = signal.cwt(sig, signal.ricker, widths)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
plt.colorbar()
plt.title('连续小波变换')
plt.ylabel('尺度')
plt.xlabel('时间')
plt.show()
6. 线性代数
scipy.linalg
模块提供了线性代数的各种功能，包括矩阵分解、线性方程求解、特征值计算等。
6.1 矩阵分解
6.1.1 LU分解
LU分解将矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积。
import numpy as np
from scipy import linalg
# 创建一个矩阵
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
# 执行LU分解
P, L, U = linalg.lu(A)
print("原始矩阵 A:")
print(A)
print("\n置换矩阵 P:")
print(P)
print("\n下三角矩阵 L:")
print(L)
print("\n上三角矩阵 U:")
print(U)
# 验证分解是否正确
print("\n验证 P @ L @ U == A:")
print(np.allclose(P @ L @ U, A))
输出结果：
原始矩阵 A:
[[ 1 2 3]
[ 4 5 6]
[ 7 8 10]]
置换矩阵 P:
[[0. 1. 0.]
[0. 0. 1.]
[1. 0. 0.]]
下三角矩阵 L:
[[1. 0. 0. ]
[0.14285714 1. 0. ]
[0.57142857 0.5 1. ]]
上三角矩阵 U:
[[7. 8. 10. ]
[0. 0.85714286 1.57142857]
[0. 0. 0.5 ]]
验证 P @ L @ U == A:
True
6.1.2 奇异值分解（SVD）
奇异值分解将矩阵分解为三个矩阵的乘积：A = UΣV*。
import numpy as np
from scipy import linalg
# 创建一个矩阵
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])
# 执行SVD分解
U, s, Vh = linalg.svd(A)
print("原始矩阵 A:")
print(A)
print("\n矩阵 U:")
print(U)
print("\n奇异值 s:")
print(s)
print("\n矩阵 Vh:")
print(Vh)
# 验证分解是否正确
# 需要将奇异值转换为对角矩阵
S = np.zeros(A.shape)
S[:A.shape[1], :A.shape[1]] = np.diag(s)
print("\n验证 U @ S @ Vh == A:")
print(np.allclose(U @ S @ Vh, A))
输出结果：
原始矩阵 A:
[[ 1 2 3]
[ 4 5 6]
[ 7 8 9]
[10 11 12]]
矩阵 U:
[[-0.14087668 0.82471435 0.54729401]
[-0.34394628 0.42626394 -0.79484036]
[-0.54701588 0.02781353 0.24754635]
[-0.75008547 -0.37063688 0. ]]
奇异值 s:
[2.54624074e+01 1.29066168e+00 7.40748856e-16]
矩阵 Vh:
[[-0.50453315 -0.5745157 -0.64449826]
[ 0.76081555 0.05714048 -0.6465346 ]
[-0.40824829 0.81649658 -0.40824829]]
验证 U @ S @ Vh == A:
True
6.2 线性方程组求解
Scipy提供了多种方法来求解线性方程组。
import numpy as np
from scipy import linalg
# 创建一个线性方程组 Ax = b
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
b = np.array([3, 6, 9])
# 求解线性方程组
x = linalg.solve(A, b)
print("系数矩阵 A:")
print(A)
print("\n常数向量 b:")
print(b)
print("\n解向量 x:")
print(x)
# 验证解是否正确
print("\n验证 Ax == b:")
print(np.allclose(A @ x, b))
输出结果：
系数矩阵 A:
[[ 1 2 3]
[ 4 5 6]
[ 7 8 10]]
常数向量 b:
[3 6 9]
解向量 x:
[-1. 1. 1.]
验证 Ax == b:
True
6.3 特征值和特征向量
特征值和特征向量是线性代数中的重要概念，在许多应用中都有使用。
import numpy as np
from scipy import linalg
# 创建一个矩阵
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# 计算特征值和特征向量
eigenvalues, eigenvectors = linalg.eig(A)
print("矩阵 A:")
print(A)
print("\n特征值:")
print(eigenvalues)
print("\n特征向量:")
print(eigenvectors)
# 验证特征值和特征向量是否正确
for i in range(len(eigenvalues)):
print(f"\n验证 A @ v_{i} == λ_{i} @ v_{i}:")
print(np.allclose(A @ eigenvectors[:, i], eigenvalues[i] * eigenvectors[:, i]))
输出结果：
矩阵 A:
[[1 2 3]
[4 5 6]
[7 8 9]]
特征值:
[ 1.61168440e+01 -1.11684397e+00 -1.30367773e-15]
特征向量:
[[-0.23197069 -0.78583024 0.40824829]
[-0.52532209 -0.08675134 -0.81649658]
[-0.8186735 0.61232756 0.40824829]]
验证 A @ v_0 == λ_0 @ v_0:
True
验证 A @ v_1 == λ_1 @ v_1:
True
验证 A @ v_2 == λ_2 @ v_2:
True
7. 插值
scipy.interpolate
模块提供了多种插值方法，可以用于在已知数据点之间估计值。
7.1 一维插值
7.1.1 线性插值
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# 创建一些数据点
x = np.linspace(0, 10, 10)
y = np.sin(x)
# 创建插值函数
f_linear = interpolate.interp1d(x, y, kind='linear')
# 生成更多的点用于绘制平滑曲线
x_new = np.linspace(0, 10, 100)
y_new = f_linear(x_new)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'o', label='数据点')
plt.plot(x_new, y_new, '-', label='线性插值')
plt.plot(x_new, np.sin(x_new), '--', label='真实函数')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('一维线性插值')
plt.grid(True)
plt.show()
7.1.2 三次样条插值
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# 创建一些数据点
x = np.linspace(0, 10, 10)
y = np.sin(x)
# 创建插值函数
f_cubic = interpolate.interp1d(x, y, kind='cubic')
# 生成更多的点用于绘制平滑曲线
x_new = np.linspace(0, 10, 100)
y_new = f_cubic(x_new)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'o', label='数据点')
plt.plot(x_new, y_new, '-', label='三次样条插值')
plt.plot(x_new, np.sin(x_new), '--', label='真实函数')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('一维三次样条插值')
plt.grid(True)
plt.show()
7.2 二维插值
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# 创建一个二维函数
def func(x, y):
return np.sin(x) * np.cos(y)
# 创建网格点
x = np.linspace(0, 2*np.pi, 10)
y = np.linspace(0, 2*np.pi, 10)
X, Y = np.meshgrid(x, y)
Z = func(X, Y)
# 创建插值函数
f = interpolate.interp2d(x, y, Z, kind='cubic')
# 生成更密集的网格点
x_new = np.linspace(0, 2*np.pi, 50)
y_new = np.linspace(0, 2*np.pi, 50)
Z_new = f(x_new, y_new)
# 绘制结果
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.pcolormesh(X, Y, Z, shading='auto')
plt.colorbar()
plt.title('原始数据')
plt.xlabel('x')
plt.ylabel('y')
plt.subplot(1, 2, 2)
plt.pcolormesh(x_new, y_new, Z_new, shading='auto')
plt.colorbar()
plt.title('二维插值结果')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
7.3 径向基函数插值
径向基函数插值适用于散乱数据点的插值。
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
# 创建一些随机散乱数据点
np.random.seed(42)
x = np.random.random(50) * 10
y = np.random.random(50) * 10
z = np.sin(x) * np.cos(y)
# 创建径向基函数插值
rbf = interpolate.Rbf(x, y, z, function='multiquadric')
# 创建网格点
xi = np.linspace(0, 10, 100)
yi = np.linspace(0, 10, 100)
XI, YI = np.meshgrid(xi, yi)
ZI = rbf(XI, YI)
# 绘制结果
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(x, y, c=z, s=50, cmap='viridis')
plt.colorbar()
plt.title('原始散乱数据点')
plt.xlabel('x')
plt.ylabel('y')
plt.subplot(1, 2, 2)
plt.pcolormesh(XI, YI, ZI, shading='auto', cmap='viridis')
plt.colorbar()
plt.title('径向基函数插值')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
8. 图像处理
scipy.ndimage
模块提供了n维图像处理功能，包括滤波、形态学操作、图像测量等。
8.1 基本图像操作
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from scipy import misc
# 加载示例图像
face = misc.face(gray=True)
# 显示原始图像
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.imshow(face, cmap='gray')
plt.title('原始图像')
plt.axis('off')
# 旋转图像
rotated = ndimage.rotate(face, 45, reshape=False)
plt.subplot(2, 2, 2)
plt.imshow(rotated, cmap='gray')
plt.title('旋转45度')
plt.axis('off')
# 缩放图像
zoomed = ndimage.zoom(face, 0.5)
plt.subplot(2, 2, 3)
plt.imshow(zoomed, cmap='gray')
plt.title('缩小50%')
plt.axis('off')
# 裁剪图像
cropped = face[100:500, 200:600]
plt.subplot(2, 2, 4)
plt.imshow(cropped, cmap='gray')
plt.title('裁剪')
plt.axis('off')
plt.tight_layout()
plt.show()
8.2 图像滤波
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from scipy import misc
# 加载示例图像并添加噪声
face = misc.face(gray=True)
noisy_face = face + 0.5 * face.std() * np.random.random(face.shape)
# 显示原始图像和带噪声图像
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.imshow(face, cmap='gray')
plt.title('原始图像')
plt.axis('off')
plt.subplot(2, 2, 2)
plt.imshow(noisy_face, cmap='gray')
plt.title('带噪声图像')
plt.axis('off')
# 高斯滤波
gaussian_filtered = ndimage.gaussian_filter(noisy_face, sigma=3)
plt.subplot(2, 2, 3)
plt.imshow(gaussian_filtered, cmap='gray')
plt.title('高斯滤波')
plt.axis('off')
# 中值滤波
median_filtered = ndimage.median_filter(noisy_face, size=5)
plt.subplot(2, 2, 4)
plt.imshow(median_filtered, cmap='gray')
plt.title('中值滤波')
plt.axis('off')
plt.tight_layout()
plt.show()
8.3 图像分割
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from scipy import misc
# 加载示例图像
face = misc.face(gray=True)
# 应用阈值分割
threshold = 100
binary = face > threshold
# 标记连通区域
label_array, num_features = ndimage.label(binary)
# 显示结果
plt.figure(figsize=(12, 8))
plt.subplot(1, 3, 1)
plt.imshow(face, cmap=’gray’)
plt.title(‘原始图像’)
plt.axis(‘off’)
plt.subplot(1, 3, 2)
plt.imshow(binary, cmap=’gray’)
plt.title(f’阈值分割 (阈值={threshold})’)
plt.axis(‘off’)
plt.subplot(1, 3, 3)
plt.imshow(label_array, cmap=’jet’)
plt.title(f’连通区域标记 ({num_features}个区域)’)
plt.axis(‘off’)
plt.tight_layout()
plt.show()
8.4 边缘检测
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from scipy import misc
# 加载示例图像
face = misc.face(gray=True)
# Sobel边缘检测
dx = ndimage.sobel(face, 0) # 水平导数
dy = ndimage.sobel(face, 1) # 垂直导数
mag = np.hypot(dx, dy) # 幅度
# 显示结果
plt.figure(figsize=(12, 8))
plt.subplot(1, 2, 1)
plt.imshow(face, cmap=’gray’)
plt.title(‘原始图像’)
plt.axis(‘off’)
plt.subplot(1, 2, 2)
plt.imshow(mag, cmap=’gray’)
plt.title(‘Sobel边缘检测’)
plt.axis(‘off’)
plt.tight_layout()
plt.show()
9. 空间数据结构与算法
scipy.spatial
模块提供了空间数据结构和算法，如KD树、凸包、Voronoi图等。
9.1 KD树
KD树是一种用于快速最近邻搜索的数据结构。
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 创建一些随机点
np.random.seed(42)
points = np.random.rand(100, 2)
# 构建KD树
tree = spatial.KDTree(points)
# 查询一个点的最近邻
point = np.array([0.5, 0.5])
distance, index = tree.query(point)
nearest_point = points[index]
print(f”查询点: {point}”)
print(f”最近邻点: {nearest_point}”)
print(f”距离: {distance}”)
# 查询k个最近邻
k = 5
distances, indices = tree.query(point, k=k)
nearest_points = points[indices]
print(f”\n查询点: {point}”)
print(f”{k}个最近邻点:”)
for i in range(k):
print(f” 点{i+1}: {nearest_points[i]}, 距离: {distances[i]}”)
# 绘制结果
plt.figure(figsize=(10, 8))
plt.scatter(points[:, 0], points[:, 1], c=’blue’, label=’数据点’)
plt.scatter(point[0], point[1], c=’red’, s=100, label=’查询点’)
plt.scatter(nearest_points[:, 0], nearest_points[:, 1], c=’green’, s=50, label=’最近邻点’)
plt.legend()
plt.grid(True)
plt.title(‘KD树最近邻搜索’)
plt.show()
9.2 凸包
凸包是包含所有点的最小凸多边形。
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 创建一些随机点
np.random.seed(42)
points = np.random.rand(30, 2)
# 计算凸包
hull = spatial.ConvexHull(points)
# 绘制结果
plt.figure(figsize=(10, 8))
plt.scatter(points[:, 0], points[:, 1], c=’blue’, label=’数据点’)
# 绘制凸包
for simplex in hull.simplices:
plt.plot(points[simplex, 0], points[simplex, 1], ‘r-‘)
# 填充凸包
plt.fill(points[hull.vertices, 0], points[hull.vertices, 1], ‘red’, alpha=0.2)
plt.legend()
plt.grid(True)
plt.title(‘凸包’)
plt.show()
9.3 Voronoi图
Voronoi图将平面划分为多个区域，每个区域包含一个生成点，且区域内任意位置到该生成点的距离小于到其他生成点的距离。
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 创建一些随机点
np.random.seed(42)
points = np.random.rand(20, 2)
# 计算Voronoi图
vor = spatial.Voronoi(points)
# 绘制结果
plt.figure(figsize=(10, 8))
spatial.voronoi_plot_2d(vor, show_vertices=False)
# 绘制生成点
plt.scatter(points[:, 0], points[:, 1], c=’red’)
plt.grid(True)
plt.title(‘Voronoi图’)
plt.show()
9.4 Delaunay三角剖分
Delaunay三角剖分是一种将点集连接成三角形网格的方法，满足所有三角形的最小角最大的特性。
import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt
# 创建一些随机点
np.random.seed(42)
points = np.random.rand(30, 2)
# 计算Delaunay三角剖分
tri = spatial.Delaunay(points)
# 绘制结果
plt.figure(figsize=(10, 8))
plt.triplot(points[:, 0], points[:, 1], tri.simplices)
plt.plot(points[:, 0], points[:, 1], ‘o’)
plt.grid(True)
plt.title(‘Delaunay三角剖分’)
plt.show()
10. 高级应用与实例
10.1 信号去噪
import numpy as np
from scipy import signal, ndimage
import matplotlib.pyplot as plt
# 创建一个信号
t = np.linspace(0, 1, 1000, False)
sig = np.sin(2*np.pi*10*t) * (t >= 0.2) * (t <= 0.8) # 10Hz正弦波，只在0.2-0.8秒内存在
sig_noisy = sig + 0.5 * np.random.randn(len(t)) # 添加噪声
# 小波去噪
wavelet = 'db4' # Daubechies小波
level = 4 # 分解层数
# 使用PyWavelets进行小波去噪
try:
import pywt
coeffs = pywt.wavedec(sig_noisy, wavelet, level=level)
threshold = np.sqrt(2*np.log(len(sig_noisy))) * np.median(np.abs(coeffs[-level])) / 0.6745
coeffs_thresh = [pywt.threshold(c, threshold, 'soft') for c in coeffs]
sig_denoised = pywt.waverec(coeffs_thresh, wavelet)
except ImportError:
print("PyWavelets未安装，使用中值滤波替代")
sig_denoised = ndimage.median_filter(sig_noisy, size=5)
# 绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(t, sig)
plt.title('原始信号')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(t, sig_noisy)
plt.title('带噪声信号')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(t, sig_denoised)
plt.title('去噪后信号')
plt.grid(True)
plt.tight_layout()
plt.show()
10.2 图像配准
import numpy as np
from scipy import ndimage, misc
import matplotlib.pyplot as plt
# 加载示例图像
image = misc.face(gray=True)
# 创建一个变换后的图像（旋转+平移）
transformed = ndimage.rotate(image, 15, reshape=False)
transformed = ndimage.shift(transformed, (20, -10))
# 显示原始图像和变换后的图像
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(image, cmap='gray')
plt.title('原始图像')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(transformed, cmap='gray')
plt.title('变换后的图像')
plt.axis('off')
plt.tight_layout()
plt.show()
# 使用相位相关进行图像配准
def phase_correlation(a, b):
# 计算二维傅里叶变换
G_a = np.fft.fft2(a)
G_b = np.fft.fft2(b)
# 计算交叉功率谱
cross_power_spectrum = (G_a * np.conj(G_b)) / (np.abs(G_a) * np.abs(G_b))
# 计算逆傅里叶变换
correlation = np.fft.ifft2(cross_power_spectrum).real
# 找到峰值位置
peak = np.unravel_index(np.argmax(correlation), correlation.shape)
# 计算位移
shift = np.array(peak) - np.array(a.shape) // 2
return shift, correlation
# 计算位移
shift, correlation = phase_correlation(image, transformed)
print(f"检测到的位移: {shift}")
# 使用检测到的位移对图像进行校正
corrected = ndimage.shift(transformed, -shift)
# 显示配准结果
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(transformed, cmap='gray')
plt.title('变换后的图像')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(corrected, cmap='gray')
plt.title('配准后的图像')
plt.axis('off')
plt.tight_layout()
plt.show()
10.3 机器学习中的Scipy应用
import numpy as np
from scipy import optimize, linalg, spatial
import matplotlib.pyplot as plt
# 生成一些示例数据
np.random.seed(42)
n_samples = 100
X = np.random.rand(n_samples, 2)
y = np.zeros(n_samples)
y[X[:, 0] + X[:, 1] > 1] = 1 # 简单的线性决策边界
# 添加一些噪声
y[np.random.choice(n_samples, 10, replace=False)] = 1 – y[np.random.choice(n_samples, 10, replace=False)]
# 定义逻辑回归模型
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def cost_function(theta, X, y, reg_param=0):
m = len(y)
h = sigmoid(X @ theta)
J = -1/m * (y.T @ np.log(h) + (1-y).T @ np.log(1-h)) + reg_param/(2*m) * np.sum(theta[1:]**2)
return J
def gradient(theta, X, y, reg_param=0):
m = len(y)
h = sigmoid(X @ theta)
grad = 1/m * X.T @ (h – y)
grad[1:] += reg_param/m * theta[1:]
return grad
# 添加偏置项
X_with_bias = np.hstack([np.ones((n_samples, 1)), X])
# 初始化参数
initial_theta = np.zeros(X_with_bias.shape[1])
# 使用优化算法找到最佳参数
result = optimize.minimize(
fun=cost_function,
x0=initial_theta,
args=(X_with_bias, y, 0.1),
method=’BFGS’,
jac=gradient
)
theta = result.x
print(“优化结果:”, result.message)
print(“参数 theta:”, theta)
# 绘制决策边界
plt.figure(figsize=(10, 8))
plt.scatter(X[y == 0, 0], X[y == 0, 1], c=’blue’, label=’类别 0′)
plt.scatter(X[y == 1, 0], X[y == 1, 1], c=’red’, label=’类别 1′)
# 计算决策边界
x_values = [np.min(X[:, 0]), np.max(X[:, 0])]
y_values = -(theta[0] + np.array(x_values) * theta[1]) / theta[2]
plt.plot(x_values, y_values, ‘g-‘, label=’决策边界’)
plt.legend()
plt.xlabel(‘特征 1’)
plt.ylabel(‘特征 2’)
plt.title(‘逻辑回归分类’)
plt.grid(True)
plt.show()
# 计算准确率
predictions = sigmoid(X_with_bias @ theta) >= 0.5
accuracy = np.mean(predictions == y)
print(f”模型准确率: {accuracy*100:.2f}%”)
10.4 数据聚类
import numpy as np
from scipy import cluster, spatial
import matplotlib.pyplot as plt
# 生成一些示例数据
np.random.seed(42)
n_samples = 300
# 生成三个不同的聚类
X1 = np.random.randn(n_samples//3, 2) + np.array([0, 0])
X2 = np.random.randn(n_samples//3, 2) + np.array([5, 5])
X3 = np.random.randn(n_samples//3, 2) + np.array([0, 5])
X = np.vstack([X1, X2, X3])
# 使用K-means聚类
k = 3 # 聚类数量
centroids, distortion = cluster.vq.kmeans(X, k)
# 分配每个点到最近的聚类中心
idx, _ = cluster.vq.vq(X, centroids)
# 绘制结果
plt.figure(figsize=(10, 8))
plt.scatter(X[:, 0], X[:, 1], c=idx, cmap=’viridis’)
plt.scatter(centroids[:, 0], centroids[:, 1], c=’red’, s=100, marker=’x’)
plt.title(‘K-means聚类’)
plt.xlabel(‘特征 1’)
plt.ylabel(‘特征 2’)
plt.grid(True)
plt.show()
# 使用层次聚类
# 计算距离矩阵
dist_matrix = spatial.distance.pdist(X, ‘euclidean’)
# 执行层次聚类
Z = cluster.hierarchy.linkage(dist_matrix, method=’ward’)
# 绘制树状图
plt.figure(figsize=(12, 6))
cluster.hierarchy.dendrogram(Z)
plt.title(‘层次聚类树状图’)
plt.xlabel(‘样本索引’)
plt.ylabel(‘距离’)
plt.show()
# 从树状图获取聚类
max_d = 10 # 最大距离阈值
clusters = cluster.hierarchy.fcluster(Z, max_d, criterion=’distance’)
# 绘制层次聚类结果
plt.figure(figsize=(10, 8))
plt.scatter(X[:, 0], X[:, 1], c=clusters, cmap=’viridis’)
plt.title(‘层次聚类’)
plt.xlabel(‘特征 1’)
plt.ylabel(‘特征 2’)
plt.grid(True)
plt.show()
总结
本教程全面介绍了Scipy库的核心功能和应用，从基础入门到高级应用，涵盖了数值优化、统计分析、信号处理、线性代数、插值、图像处理、空间数据结构等多个方面。通过大量的实例演示，我们展示了如何使用Scipy解决实际问题。
Scipy作为Python科学计算生态系统的重要组成部分，为科学研究和工程应用提供了强大的工具。无论是进行数据分析、信号处理、图像处理还是机器学习，Scipy都能提供必要的支持。
要真正掌握Scipy，除了学习本教程中的内容外，还需要在实际项目中不断实践和探索。Scipy的官方文档（
https://docs.scipy.org/doc/scipy/）是获取更多信息和深入学习的宝贵资源。
希望本教程能够帮助你快速掌握Scipy的核心技术，并在你的科学计算和数据分析工作中发挥重要作用。