Python学习-Scipy库优化与拟合optimize(最小二乘法拟合、B-样条拟合)

Python学习-Scipy库优化与拟合optimize

目录

1、最小二乘法拟合least_squares()

2、B-样条拟合interpolate.BSpline()

导入库

import scipy.optimize as otm
import scipy.interpolate as ipl
import numpy as np
import matplotlib.pyplot as plt

plt.rc('font', family='simhei', size=15)  # 设置中文显示，字体大小
plt.rc('axes', unicode_minus=False)  # 该参数解决负号显示的问题

1、最小二乘法拟合

通过最小化误差的平方和寻找实际的最佳逼近值，把这些点连接起来形成拟合曲线。
最小二乘法线性回归拟合，线性回归的目的就是实现预测值与实际值的残差平方和最小。

求解带变量边界的非线性最小回归问题 least_squares()

重要参数说明：
fun: 计算残差向量的函数，形式为fun(x0,args,**kwargs), 传递给x0，为多维数组
x0: 多维数组、浮点数，对自变量的初步猜测值，如果是浮点数，它将被视为带有一个元素的一维数组
jac: ‘2-point’、‘3-point’、‘cs’、‘callable’，计算雅可比矩阵的方法（mn矩阵，其中元素(i,j)是f[i]相当于x[j]的偏导数）
bounds: 带2元组的数组对象，设置自变量的下限和上限，默认无界限。每个数组必须匹配x0参数的大小或标量
method: ‘trf’，Trust Region Reflective 算法，适用于带边界的大型稀疏问题；
‘dogbox’，具有矩形信任区域的狗腿算法，适用于边界小问题
‘lm’，在MINPACK中实现的Levenberg-Marquardt 算法，适用于小的无约束问题
ftol: 浮点数，通过改变成本函数来容忍终止
xtol: 浮点数，通过改变自变量来容忍终止
gtol: 浮点数，通过梯度的范数容忍终止
args: 元组、字典，传递给fun函数和jac的其他函数

返回值说明：
x: 多维数组，最小二乘法拟合求解方案
cost: 浮点数，解决方案中成本函数的值
fun: 多维数组，解决方案中的残差
jac: 多维数组，稀疏矩阵或线性
grad: 多维数组，解决方案中成本函数的梯度
optimality: 浮点数，一阶最优度量：在无约束问题中，是梯度的统一规范；在约束问题中，是在迭代期间与gtol进行比较的数量
active_mask: 基于整型的多维数组，显示method方法相应的约束是否处于活动状态，0未激活、-1下限有效、1上限有效
nfev: 整数，完成的功能评估数
njev: 整数或None，雅可比评估的数量
status: 算法终止的原因
message: 提供终止原因的描述
success: 如果执行结果中满足其中一个收敛条件，为True

代码：

Y = np.array([7, 16, 10, 19, 12, 20.5, 16, 25, 17, 27, 18, 27, 20, 31])
X = np.array([0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6])


def p(x):
    a, b = x
    return Y-(a*X+b)


line = otm.least_squares(p, [1, 0])
print(line)

# 绘图
y = (Y+line.fun)
plt.figure(1)
plt.plot(X, Y, 'o', label='温度实际测量值')
plt.plot(X, y, 'g-', label='带残差值温度曲线')
plt.plot([0, 6], line.x+11.5, 'b--', lw=2, label='线性拟合线')
plt.legend()
plt.show()

输出：

2、B-样条拟合

B-样条曲线曲面具有几何不变性、凸包性、保凸性、变差减小性、局部支撑性等许多优良性质

interpolate.BSpline() 函数

参数说明：
t: 多维数组，形状为（n+k+1，），指定原始节点数
c: 多维数组，形状为（>=n, …），样条系数
k: 整数，指定样条阶数
extrapolate: 布尔值或’periodic’, 可选推断是否超出基本区间t[k]…t[n];
默认值为True，则推断基本区间上活动的第一个和最后一个多项式的B-样条函数；
‘periodic’，则使用周期性外推法
axis: 默认值为0，指定插值轴方向

代码：

def B(x, k, i, t):  # 递归定义函数
    if k == 0:
        return 1.0 if t[i] <= x < t[i + 1] else 0.0
    if t[i+k] == t[i]:
        c1 = 0.0
    else:
        c1 = (x - t[i])/(t[i+k]-t[i]) * B(x, k-1, i, t)
    if t[i+k+1] == t[i+1]:
        c2 = 0.0
    else:
        c2 = (t[i+k+1] - x) / (t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t)
    return c1 + c2


def bspline(x, t, c, k):  # 建立一个B-样条评估函数
    n = len(t) - k - 1
    assert (n >= k+1) and (len(c) >= n)
    return sum(c[i] * B(x, k, i, t) for i in range(n))


# 在基础区间2<= x<=4 上构造二次样条函数，并与评估样条的简单方法进行比较
k = 2  # 指定2阶
x = np.linspace(0, np.pi, 7)
t = x * 1.8  # 原始采样数据
c = [-1, 2, 0, -1]  # 样条系数
spl = ipl.BSpline(t, c, k)  # B-样条拟合计算
spl(2.5)
bspline(2.5, t, c, k)  # 自定义评估B-样条曲线
fig, ax = plt.subplots()
xx = np.linspace(1.5, 4.5, 50)
ax.plot(xx, [bspline(x, t, c, k) for x in xx], 'r--', lw=3, label='自定义B-样条评估曲线')
ax.plot(xx, spl(xx), 'b-', lw=4, label='自定义B-样条拟合曲线')
ax.grid(True)
ax.legend(loc='best')
plt.show()

输出：