拟合方法——leastsq
1. 概念:
scipy官网对该方法介绍是: 最小化一组方程的平方和
x
=
arg
min
y
(
∑
(
(
f
u
n
c
(
y
)
)
2
，
a
x
i
s
=
0
)
)
x =\arg \min\limits_{y}(\sum((func(y))^2，axis=0))
x=argymin(∑((func(y))2，axis=0))
简单介绍一下leastsq的参数:scipy.optimize.leastsq（func，x0，args = （），Dfun = None，full_output = 0，col_deriv = 0，ftol = 1.49012e-08，xtol = 1.49012e-08，gtol = 0.0，maxfev = 0，epsfcn = None，factor = 100，diag = None）
# 参数:
func: 所求平方和最小值的函数，一般都用于拟合时传入残差函数，当然也有其他用途，比如求函数的最小等等
x0: ndarray 一般为参数的初始值
args: 元组，可选，可以传入(x, y)其中x,y都是数组，作为求解传入的参数
# 返回值
leastsq的返回值是一个元组，里面包含了很多信息，如果是用作拟合的话，返回值的索引值为0的元素即为我们所求函数参数
有兴趣的读者可以自己看看官网介绍: leastsq官网
2. 使用:
1. 求最小值
比如求 2 ( x − 3 ) 2 + 1 2(x-3)^2+1 2(x−3)2+1的最小值:
from scipy.optimize import leastsq
def func(x):
return 2*(x-3)**2+1
print(leastsq(func, 0)) # 0为x的初始值，初始值可以由自己定,改成任意一个随机值都可以
结果为:
(array([2.99999999]), 1)
当
x
=
2.99999999
≈
3
x=2.99999999\approx3
x=2.99999999≈3时，函数最小值为1
2. 拟合
为了和上一篇文章数学建模方法—【04】拟合方法之np.polyfit、np.poly1d形成一个比较，我们这里也进行1-5阶的多项式拟合。
- 定义要拟合的函数关系式 (因为我们这里函数关系式是多项式，所以就用np.ploy1d获取了多项式了)
如果想要拟合的是其他函数，那只需要将from scipy.optimize import leastsq # 拟合数据集 x = [4, 8, 10, 12, 25, 32, 43, 58, 63, 69, 79] y = [20, 33, 50, 56, 42, 31, 33, 46, 65, 75, 78] x_arr = np.array(x) y_arr = np.array(y) def fitting_func(p, x): """ 获得拟合的目标数据数组 :param p: array[int] 多项式各项从高到低的项的系数数组 :param x: array[int] 自变量数组 :return: array[int] 拟合得到的结果数组 """ f = np.poly1d(p) # 获得拟合后得到的多项式 return f(x) # 将自变量数组带入多项式计算得到拟合所得的目标数组f = np.poly1d(p)改为你想要拟合的函数关系式即可。 - 定义残差计算函数
残差有很多计算方法，这里因为leastsq方法本身就会将我们的残差函数平方，所以这里我们直接相减，在经过leastsq方法的平方后就变为了最小二乘法了。def error_func(p, x, y): """ 计算残差 :param p: array[int] 多项式各项从高到低的项的系数数组 :param x: array[int] 自变量数组 :param y: array[int] 原始目标数组(因变量) :return: 拟合得到的结果和原始目标的差值 """ err = fitting_func(p, x) - y return err - 获取参数
def n_poly(n, x, y): """ n 次多项式拟合函数 :param n: 多项式的项数(包括常数项)，比如n=3的话最高次项为2次项 :return: 最终得到的系数数组 """ p_init = np.random.randn(n) # 生成 n个随机数，作为各项系数的初始值，用于迭代修正 parameters = leastsq(error_func, p_init, args=(np.array(x), np.array(y))) # 三个参数：误差函数、函数参数列表、数据点 return parameters[0] # leastsq返回的是一个元组，元组的第一个元素时多项式系数数组[wn、...、w2、w1、w0] - 函数定义好后，就可以进行拟合了(这里计算拟合优度的方法goodness_of_fit在我前面的文章有: 数学建模方法—【03】拟合优度的计算(python计算))
一阶
二、三、四、五阶:para_2 = n_poly(2, x_arr, y_arr) # 一阶 print("系数为:", para_2) print("多项式为:", np.poly1d(para_2)) print("拟合值为:", fitting_func(para_2, x_arr)) print("拟合优度为:", goodness_of_fit(fitting_func(para_2, x_arr), y_arr))
这里解释一下为什么一阶传入n=2，二阶传入n=3…。因为这里的n表示的不是多项式的阶数，而是表示多项式的项数，比如二阶的多项式有 x 2 、 x x^2、x x2、x还有常数项三项，所以传入的n=3。para_3 = n_poly(3, x_arr, y_arr) # 二阶 para_4 = n_poly(4, x_arr, y_arr) # 三阶 para_5 = n_poly(5, x_arr, y_arr) # 四阶 para_6 = n_poly(6, x_arr, y_arr) # 五阶 - 得到拟合的数据后就可以作图了
这里我只输出了一阶的拟合结果，结果为:figure = plt.figure(figsize=(8,6)) plt.plot(x_arr, fitting_func(para_2, x_arr), color="#72CD28", label='一阶拟合曲线') plt.plot(x_arr, fitting_func(para_3, x_arr), color="#EBBD43", label='二阶拟合曲线') plt.plot(x_arr, fitting_func(para_4, x_arr), color="#50BFFB", label='三阶拟合曲线') plt.plot(x_arr, fitting_func(para_5, x_arr), color="gold", label='四阶拟合曲线') plt.plot(x_arr, fitting_func(para_6, x_arr), color="#13EAC9", label='五阶拟合曲线') plt.scatter(x, y, color='#50BFFB', marker="p", label='原始数据') plt.title("leastsq拟合", fontsize=13) plt.xlabel('x', fontsize=13) plt.ylabel('y', fontsize=13) plt.legend(loc=4, fontsize=11) # 指定legend的位置右下角 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 plt.show()
结果图为:系数为: [ 0.5000123 29.77227653] 多项式为: 0.5 x + 29.77 拟合值为: [31.77232574 33.77237495 34.77239955 35.77242415 42.27258408 45.77267019 51.27280552 58.77299004 61.27305155 64.27312537 69.27324839] 拟合优度为: 0.5207874477394239
全部代码为:
from scipy.optimize import leastsq
# 拟合数据集
x = [4, 8, 10, 12, 25, 32, 43, 58, 63, 69, 79]
y = [20, 33, 50, 56, 42, 31, 33, 46, 65, 75, 78]
x_arr = np.array(x)
y_arr = np.array(y)
def fitting_func(p, x):
"""
获得拟合的目标数据数组
:param p: array[int] 多项式各项从高到低的项的系数数组
:param x: array[int] 自变量数组
:return: array[int] 拟合得到的结果数组
"""
f = np.poly1d(p) # 获得拟合后得到的多项式
return f(x) # 将自变量数组带入多项式计算得到拟合所得的目标数组
def error_func(p, x, y):
"""
计算残差
:param p: array[int] 多项式各项从高到低的项的系数数组
:param x: array[int] 自变量数组
:param y: array[int] 原始目标数组(因变量)
:return: 拟合得到的结果和原始目标的差值
"""
err = fitting_func(p, x) - y
return err
def n_poly(n, x, y):
"""
n 次多项式拟合函数
:param n: 多项式的项数(包括常数项)，比如n=3的话最高次项为2次项
:return: 最终得到的系数数组
"""
p_init = np.random.randn(n) # 生成 n个随机数，作为各项系数的初始值，用于迭代修正
parameters = leastsq(error_func, p_init, args=(np.array(x), np.array(y))) # 三个参数：误差函数、函数参数列表、数据点
return parameters[0] # leastsq返回的是一个元组，元组的第一个元素时多项式系数数组[wn、...、w2、w1、w0]
para_2 = n_poly(2, x_arr, y_arr) # 一阶
print("系数为:", para_2)
print("多项式为:", np.poly1d(para_2))
print("拟合值为:", fitting_func(para_2, x_arr))
print("拟合优度为:", goodness_of_fit(fitting_func(para_2, x_arr), y_arr))
para_3 = n_poly(3, x_arr, y_arr) # 二阶
para_4 = n_poly(4, x_arr, y_arr) # 三阶
para_5 = n_poly(5, x_arr, y_arr) # 四阶
para_6 = n_poly(6, x_arr, y_arr) # 五阶
figure2 = plt.figure(figsize=(8,6))
plt.plot(x_arr, fitting_func(para_2, x_arr), color="#72CD28", label='一阶拟合曲线')
plt.plot(x_arr, fitting_func(para_3, x_arr), color="#EBBD43", label='二阶拟合曲线')
plt.plot(x_arr, fitting_func(para_4, x_arr), color="#50BFFB", label='三阶拟合曲线')
plt.plot(x_arr, fitting_func(para_5, x_arr), color="gold", label='四阶拟合曲线')
plt.plot(x_arr, fitting_func(para_6, x_arr), color="#13EAC9", label='五阶拟合曲线')
plt.scatter(x_arr, y_arr, color='#50BFFB', marker="p", label='原始数据')
plt.title("leastsq拟合", fontsize=13)
plt.xlabel('x', fontsize=13)
plt.ylabel('y', fontsize=13)
plt.legend(loc=4, fontsize=11) # 指定legend的位置右下角
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
plt.show()
