科学计算的核心包——SciPy

SciPy是一个用于数学、科学及工程方面的常用软件包，Scipy包含科学计算中常见各种工具箱。它不同的子模块对应不同的应用，如插值、积分、优化、图像处理等。SciPy可以与其他标准科学计算程序库进行比较，如GSL（GNU C或C++科学计算库），或者Matlab工具箱。
SciPy是Python中科学计算程序的核心包，SciPy函数库在NumPy库的基础上增加了众多的数学、科学及工程计算中常用的库函数，如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等。它用于有效地计算NumPy矩阵，来让NumPy和SciPy协同工作。

最小二乘法

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。

实例：使用最小二乘法进行直线拟合。

import numpy as np
from scipy.optimize import leastsq
# 已知的训练数据
Xi = np.array([8.19, 2.72, 6.39, 8.71, 4.7, 2.66, 3.78])
Yi = np.array([7.01, 2.78, 6.47, 6.71, 4.1, 4.23, 4.05])
# 定义要拟合的函数形式和误差
def func(p, x):
    k, b = p
    return k * x + b

def error(p, x, y, s):
    print(s)
    return func(p, x) - y

# 随机给出初始值
p0 = [100, 2]
# 利用SciPy的leastsq函数估计参数
s = "Test the number of iteration"
Para = leastsq(error, p0, args=(Xi, Yi, s))
k, b = Para[0]
print("k=", k, '\n', "b=", b)

估计的结果如下：

计算结果

为了更加直观地观察拟合效果，可以使用图形显示，代码如下。

import matplotlib.pyplot as plt
plt.figure(figsize=(8, 6))
plt.scatter(Xi, Yi, color="red", label="Sample Point", linewidth=3)
x = np.linspace(0, 10, 1000)
y = k * x + b
plt.plot(x, y, color="orange", label="Fitting Line", linewidth=2)
plt.legend()
plt.show()

运行结果如下：

拟合结果

实例：使用最小二乘法进行正弦函数的拟合。

在这个例子中要拟合的函数是一个正弦波函数，他有3个参数A、k、theta，分别对应振幅、频率和相角。假设实验数据是一组包含噪声的数据x,y1，其中y1是在真实数据y0的基础上加入噪声得到的。

import numpy as np
from scipy.optimize import leastsq
import pylab as pl

def func(x, p):
    A, k, theta = p
    return A * np.sin(2 * np.pi * k * x + theta)

def residuals(p, y, x):
    return y - func(x, p)

# 训练数据的生成如下
x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6
y0 = func(x, [A, k, theta])
y1 = y0 + 2 * np.random.randn(len(x))
p0 = [7, 0.2, 0]  # 第一次猜测的函数拟合参数
plsq = leastsq(residuals, p0, args=(y1, x))
print("真实参数：", [A, k, theta])
print("拟合参数：", plsq[0])  # 实验数据拟合后的参数

pl.plot(x, y0, label='real data')
pl.plot(x, y1, label='data having noise')
pl.plot(x, func(x, plsq[0]), label="fitting data")
pl.legend()
pl.show()

最终训练的参数如下：

拟合结果

为了更好地显示结果，使用下图所示的图形显示。

绘图结果

非线性方程组求解

optimize库中的fsolve函数可以用来对非线性方程组进行求解。它的基本调用形式为：fsolve(func, x0)，其中func(x)是计算方程组的函数，它的参数x是一个矢量，表示方程组的各个未知数的一组可能解，func返回将x代入方程组之后得到的结果；x0为未知数矢量的初始值。

实例：解一个非线性方程组。

from scipy.optimize import fsolve
from math import sin, cos

def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
            5 * x1 + 3,
            4 * x0 * x0 - 2 * sin(x1 * x2),
            x1 * x2 - 1.5
            ]
result = fsolve(f, [1,1,1])

print(result)
print(f(result))  

最终计算的结果为：

[-0.70622057 -0.6        -2.5       ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]