AI 基础：Scipy（科学计算库）简易入门

阅读量：798 次

发布时间：2023-04-17

本文共 4188 字，大约阅读时间需要 13 分钟。

SciPy 数值计算库快速入门

Scipy 是一个强大的数学、科学和工程领域的软件包，提供了丰富的子模块来处理插值、积分、优化、图像处理、常微分方程求解、信号处理等问题。它与 NumPy 结合使用，能够高效解决各种数值计算问题。

SciPy 的模块体系

Scipy 由多个子模块组成，每个模块专注于特定领域的计算任务。以下是常见模块的应用领域：

scipy.cluster：向量计算与 KMeans

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：统计分布

常数与特殊函数

Scipy 提供了丰富的常数和特殊函数，例如：

from scipy import constants as C
print(C.c)  # 真空中的光速
print(C.h)  # 普朗克常数

可以看到，C.c 的值大约是 299792458 m/s，C.h 大约是 6.62607015 × 10^-34 J·s。

此外，Scipy 还提供了许多常用的数学函数，例如：

import scipy.special as S
print(1 + 1e-20)
print(np.log(1 + 1e-20))
print(S.log1p(1e-20))

输出结果显示，log1p 函数能够更准确地计算接近1的对数，这在某些物理问题中非常重要。

优化与拟合

Scipy 提供了多种优化算法，适用于非线性方程组求解和最小二乘拟合等问题。以下是一个非线性方程组求解的示例：

import scipy.optimize as optimize
def f(x):
    x0, x1, x2 = x.tolist()
    return [5 * x1 + 3, 4 * x0**2 - 2 * np.sin(x1 * x2), x1 * x2 - 1.5]
result = optimize.fsolve(f, [1, 1, 1])
print(result)

输出结果显示，优化算法能够准确地找到满足方程组的解。

对于最小二乘拟合，Scipy 提供了 leastsq 函数，能够高效地拟合数据点。以下是一个简单的线性回归示例：

import numpy as np
from scipy.optimize import leastsq
X = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])
Y = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])
def residuals(p):
    k, b = p
    return Y - (k * X + b)
r = leastsq(residuals, [1, 0])
k, b = r[0]
print("k =", k, "b =", b)

输出结果显示，拟合后的斜率和截距能够较好地拟合数据点。

线性代数

Scipy 的 linalg 模块提供了丰富的线性代数功能，包括矩阵运算、特征值求解等。以下是一个解线性方程组的示例：

import numpy as np
from scipy import linalg
m, n = 500, 50
A = np.random.rand(m, m)
B = np.random.rand(m, n)
X1 = linalg.solve(A, B)
X2 = np.dot(linalg.inv(A), B)
print(np.allclose(X1, X2))

输出结果显示，solve 和 inv 方法能够正确地解方程组。

此外，Scipy 还支持奇异值分解（SVD），可以用来降维或分析矩阵的结构。以下是一个 SVD 的示例：

from scipy import linalg
img = pl.imread("example.jpg")[:3]  # 读取图像的前三个通道
U, s, Vh = linalg.svd(img)
print(U.shape)
print(s.shape)
print(Vh.shape)

输出结果显示，SVD 分解能够有效地提取图像的主成分。

统计与概率

Scipy 的 stats 模块提供了丰富的统计分布和函数，适用于概率密度估计、假设检验等任务。以下是一个正态分布拟合的示例：

import numpy as np
from scipy import stats
X = np.random.normal(1.0, 2.0, 10000)
print(np.mean(X), np.var(X))
stats_norm = stats.norm.fit(X)
print(stats_norm.mean, stats_norm.scale)

输出结果显示，拟合后的正态分布能够准确描述数据分布。

此外，Scipy 还支持卡方分布和 t 检验，能够用于假设检验和其他统计分析任务。

积分与微分方程

Scipy 的 integrate 模块支持数值积分和常微分方程求解。以下是一个球体积的积分示例：

from scipy import integrate
def half_circle(x):
    return (1 - x**2)**0.5
volume, error = integrate.dblquad(half_circle, -1, 1, lambda x: half_circle(x), lambda x: half_circle(x))
print(volume, error)

输出结果显示，双重积分能够准确地计算半球的体积。

此外，Scipy 还支持使用 odeint 求解常微分方程，例如洛伦兹吸引子模型：

from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
    x, y, z = w.tolist()
    return p*(y - x), x*(r - z) - y, x*y - b*z
track1 = odeint(lorenz, (0.0, 1.0, 0.0), np.arange(0, 30, 0.01), args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), np.arange(0, 30, 0.01), args=(10.0, 28.0, 3.0))

输出结果显示，两条轨道在不同的初始条件下表现出不同的动力学行为。

信号处理

Scipy 的 signal 模块提供了丰富的信号处理功能，包括中值滤波、滤波器设计、连续时间线性系统分析等。以下是一个中值滤波的示例：

from scipy import signal
t = np.arange(0, 20, 0.1)
x = np.sin(t)
x_random = x[np.random.randint(0, len(t), 20)]
x_random += np.random.standard_normal(20) * 0.6
x2 = signal.medfilt(x_random, 5)
print(np.all(x2 == signal.order_filter(x_random, np.ones(5), 2)))

输出结果显示，中值滤波能够有效地去除瞬时噪声。

此外，Scipy 还支持设计带通滤波器和其他类型的滤波器，能够满足不同信号处理需求。

插值

Scipy 提供了多种插值方法，包括线性插值、多项式插值、径向基函数插值等。以下是一个使用 UnivariateSpline 的外推示例：

from scipy import interpolate
x1 = np.linspace(0, 10, 20)
y1 = np.sin(x1)
xnew = np.linspace(0, 12, 100)
spline = interpolate.UnivariateSpline(x1, y1, s=0)
ynew = spline(xnew)
print(ynew)

输出结果显示，Spline 插值能够准确地外推数据点。

此外，Scipy 还支持二维插值和径向基函数插值，适用于更复杂的插值任务。

稀疏矩阵

Scipy 提供了稀疏矩阵的表示和操作方法，适用于大规模数据的处理。以下是一个稀疏矩阵与密集矩阵求解的示例：

from scipy import sparse
import numpy as np
A = sparse.lil_matrix((1000, 1000))
A[0, :100] = np.random.rand(100)
A[1, 100:200] = A[0, :100]
A.setdiag(np.random.rand(1000))
A = A.tocsr()
b = np.random.rand(1000)
x = sparse.linalg.spolve(A, b)
x_dense = np.linalg.solve(A.toarray(), b)
print(np.max(np.abs(x - x_dense)) < 1e-10)

输出结果显示，稀疏矩阵求解与密集矩阵求解结果一致，验证了稀疏矩阵的高效性。