赞
踩
为了方便研究最小二乘法问题,现提供如下车辆时间和行驶距离的观测数据用于讨论和分析。
令:车辆的初速度
v
0
v_0
v0 、加速度
a
a
a、时间
t
t
t、距离
y
^
\hat{y}
y^ ,第
i
i
i组观测值为
(
t
i
,
y
^
i
)
(t_i,\hat{y}_i )
(ti,y^i),时间
t
t
t与距离
y
^
\hat{y}
y^满足的数学模型如下:
y
^
=
f
(
v
0
,
a
,
t
)
=
1
2
a
t
2
+
v
0
t
(
公式
59
)
\hat{y}=f(v_0,a,t)=\frac{1}{2}at^2+v_0t \qquad (公式59)
y^=f(v0,a,t)=21at2+v0t(公式59)
min
F
(
a
,
v
0
)
=
1
2
∑
i
=
1
10
(
f
(
v
0
,
a
,
t
i
)
−
y
i
^
)
2
(
公式
60
)
\text{min}F(a,v_0)=\frac{1}{2}\sum_{i=1}^{10} (f(v_0,a,t_i)-\hat{y_i})^2 \qquad (公式60)
minF(a,v0)=21i=1∑10(f(v0,a,ti)−yi^)2(公式60)
观察公式 60 发现,
y
i
^
\hat{y_i}
yi^ 和
t
i
t_i
ti 是 常 数 ,
a
a
a 和
v
0
v_0
v0 是变量,可知该最小二乘问题是函数
F
F
F关于
a
a
a 和
v
v
v 的最小极值问题。
import numpy as np import matplotlib.pyplot as plt # 定义数学模型 def f(args, t): return 0.5 * args[0] * np.power(t, 2) + args[1] * t # 定义残差函数 def l(args, t, y): return f(args, t) - y # 定义目标函数 def o(lv): return 0.5 * np.power(lv, 2) # 定义关于目标函数第一个参数的一阶偏导函数 def oj1(args, t, y): return (o(l(args + [delta, 0], t, y)) - o(l(args, t, y))) / delta # 定义关于目标函数第二个参数的一阶偏导函数 def oj2(args, t, y): return (o(l(args + [0, delta], t, y)) - o(l(args, t, y))) / delta # 初始数据 t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) # 加速度实验的观测值,时间(单位:秒) y = np.array([0, 11.52, 26.2, 43.5, 64.12, 87.57, 114.12, 143.5, 176.3, 211.5, 250.12]) # 加速度实验的观测值,距离(单位:米) x = np.array([1, 1]) # 初始化参数[args(1);args(2)] goal = 0.0001 # 目标误差 epochs = 10000 # 训练次数 delta = 0.0001 # 偏导计算的增量 ov = 0 # 初始化目标值 for k in range(epochs): # 误差计算 tl = l(x, t, y) # 计算残差值 tmp_ov = np.sum(o(tl)) # 计算最新目标值 tmp_error = np.abs(tmp_ov - ov) # 目标值的变化 print(f'{k}->目标值: [{tmp_ov}], 误差: [{tmp_error}]') ov = tmp_ov # 更新目标值 if tmp_error < goal: print(f'训练次数: {k+1}次...') break # 执行训练 J = np.array([np.sum(oj1(x, t, y)), np.sum(oj2(x, t, y))]) # 雅克比矩阵 # 更新参数 x = x - 0.0001 * J # 打印训练结果 print(f'训练后的参数: {x}') plt.plot(t, y, 'o', np.arange(min(t), max(t), 0.001), f(x, np.arange(min(t), max(t), 0.001)), '-') plt.show()
2301->目标值: [0.06649798204087369], 误差: [0.00010499914693314072]
2302->目标值: [0.06639346785247356], 误差: [0.00010451418840012883]
2303->目标值: [0.06628943641358499], 误差: [0.0001040314388885688]
2304->目标值: [0.06618588552519271], 误差: [0.00010355088839228421]
2305->目标值: [0.06608281299820142], 误差: [0.00010307252699129354]
2306->目标值: [0.06598021665344106], 误差: [0.00010259634476035562]
2307->目标值: [0.06587809432158329], 误差: [0.00010212233185777353]
2308->目标值: [0.0657764438431518], 误差: [0.00010165047843148367]
2309->目标值: [0.06567526306842077], 误差: [0.0001011807747310356]
2310->目标值: [0.0655745498573545], 误差: [0.00010071321106626396]
2311->目标值: [0.06547430207969993], 误差: [0.00010024777765457737]
2312->目标值: [0.06537451761472354], 误差: [9.978446497638238e-05]
训练次数: 2313次...
训练后的参数: [3.00520621 9.99137523]
import numpy as np import matplotlib.pyplot as plt # 定义数学模型 def f(args, t): return 0.5 * args[0] * np.power(t, 2) + args[1] * t # 定义残差函数 def l(args, t, y): return f(args, t) - y # 定义目标函数 def o(lv): return 0.5 * np.power(lv, 2) # 定义关于目标函数第一个参数的一阶偏导函数 def oj1(args, t, y): return (o(l(args + [delta, 0], t, y)) - o(l(args, t, y))) / delta # 定义关于目标函数第二个参数的一阶偏导函数 def oj2(args, t, y): return (o(l(args + [0, delta], t, y)) - o(l(args, t, y))) / delta # 定义关于目标函数的(第一个参数,第一个参数)的二阶偏导函数 def oh00(args, t, y): return (oj1(args + [delta, 0], t, y) - oj1(args, t, y)) / delta # 定义关于目标函数的(第一个参数,第二个参数)的二阶偏导函数 def oh01(args, t, y): return (oj1(args + [0, delta], t, y) - oj1(args, t, y)) / delta # 定义关于目标函数的(第二个参数,第一个参数)的二阶偏导函数 def oh10(args, t, y): return (oj2(args + [delta, 0], t, y) - oj2(args, t, y)) / delta # 定义关于目标函数的(第二个参数,第二个参数)的二阶偏导函数 def oh11(args, t, y): return (oj2(args + [0, delta], t, y) - oj2(args, t, y)) / delta # 初始数据 t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) # 加速度实验的观测值,时间(单位:秒) y = np.array([0, 11.52, 26.2, 43.5, 64.12, 87.57, 114.12, 143.5, 176.3, 211.5, 250.12]) # 加速度实验的观测值,距离(单位:米) x = np.array([1, 1]) # 初始化参数[args(1);args(2)] goal = 0.0001 # 目标误差 epochs = 10000 # 训练次数 delta = 0.0001 # 偏导计算的增量 ov = 0 # 初始化目标值 for k in range(epochs): # 误差计算 tl = l(x, t, y) # 计算残差值 tmp_ov = np.sum(o(tl)) # 计算最新目标值 tmp_error = np.abs(tmp_ov - ov) # 目标值的变化 print(f'{k}->目标值: [{tmp_ov}], 误差: [{tmp_error}]') ov = tmp_ov # 更新目标值 if tmp_error < goal: print(f'训练次数: {k+1}次...') break # 执行训练 H = np.array([[np.sum(oh00(x, t, y)), np.sum(oh01(x, t, y))], [np.sum(oh10(x, t, y)), np.sum(oh11(x, t, y))]]) # 黑塞矩阵 J = np.array([np.sum(oj1(x, t, y)), np.sum(oj2(x, t, y))]) # 雅克比矩阵 # 更新参数 x = x - np.linalg.inv(H) @ J # 打印训练结果 print(f'训练后的参数: {x}') plt.plot(t, y, 'o', np.arange(min(t), max(t), 0.001), f(x, np.arange(min(t), max(t), 0.001)), '-') plt.show()
0->目标值: [55574.229250000004], 误差: [55574.229250000004]
1->目标值: [0.04454421541289134], 误差: [55574.18470578459]
2->目标值: [0.04452996893422017], 误差: [1.4246478671167684e-05]
训练次数: 3次...
训练后的参数: [ 2.99458661 10.0356844 ]
import numpy as np import matplotlib.pyplot as plt # 定义数学模型 def f(args, t): return 0.5 * args[0] * np.power(t, 2) + args[1] * t # 定义残差函数 def l(args, t, y): return f(args, t) - y # 定义目标函数 def o(lv): return 0.5 * np.power(lv, 2) # 定义关于残差函数第一个参数的一阶偏导函数 def j1(args, t, y): return (l(args + [delta, 0], t, y) - l(args, t, y)) / delta # 定义关于残差函数第二个参数的一阶偏导函数 def j2(args, t, y): return (l(args + [0, delta], t, y) - l(args, t, y)) / delta # 初始数据 t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) # 加速度实验的观测值,时间(单位:秒) y = np.array([0, 11.52, 26.2, 43.5, 64.12, 87.57, 114.12, 143.5, 176.3, 211.5, 250.12]) # 加速度实验的观测值,距离(单位:米) x = np.array([1, 1]) # 初始化参数[args(1);args(2)] goal = 0.0001 # 目标误差 epochs = 10000 # 训练次数 delta = 0.0001 # 偏导计算的增量 ov = 0 # 初始化目标值 for k in range(epochs): # 误差计算 tl = l(x, t, y) # 计算残差值 tmp_ov = np.sum(o(tl)) # 计算最新目标值 tmp_error = np.abs(tmp_ov - ov) # 目标值的变化 print(f'{k}->目标值: [{tmp_ov}], 误差: [{tmp_error}]') ov = tmp_ov # 更新目标值 if tmp_error < goal: print(f'训练次数: {k+1}次...') break # 执行训练 tmp_j1 = j1(x, t, y) # 关于第一个参数的一阶偏导函数的值 tmp_j2 = j2(x, t, y) # 关于第二个参数的一阶偏导函数的值 H = np.array([[np.sum(tmp_j1 * tmp_j1), np.sum(tmp_j1 * tmp_j2)], [np.sum(tmp_j2 * tmp_j1), np.sum(tmp_j2 * tmp_j2)]]) # 近似黑塞矩阵 J = np.array([np.sum(tl * tmp_j1), np.sum(tl * tmp_j2)]) # 近似雅克比矩阵 # 更新参数 x = x - np.linalg.inv(H) @ J # 打印训练结果 print(f'训练后的参数: {x}') plt.plot(t, y, 'o', np.arange(min(t), max(t), 0.001), f(x, np.arange(min(t), max(t), 0.001)), '-') plt.show()
0->目标值: [55574.229250000004], 误差: [55574.229250000004]
1->目标值: [0.04445524644031092], 误差: [55574.184794753564]
2->目标值: [0.044455246440310965], 误差: [4.85722573273506e-17]
训练次数: 3次...
训练后的参数: [ 2.99520263 10.03331435]
import numpy as np import matplotlib.pyplot as plt # 定义数学模型 def f(args, t): return 0.5 * args[0] * np.power(t, 2) + args[1] * t # 定义残差函数 def l(args, t, y): return f(args, t) - y # 定义目标函数 def o(lv): return 0.5 * np.power(lv, 2) # 偏导函数 delta = 0.0001 def j1(args, t, y): return (l(args + [delta, 0], t, y) - l(args, t, y)) / delta def j2(args, t, y): return (l(args + [0, delta], t, y) - l(args, t, y)) / delta # 初始数据 t = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) # 加速度实验的观测值,时间(单位:秒) y = np.array([0, 11.52, 26.2, 43.5, 64.12, 87.57, 114.12, 143.5, 176.3, 211.5, 250.12]) # 加速度实验的观测值,距离(单位:米) x = np.array([1, 1]) # 初始化参数[args(1);args(2)] goal = 0.0001 # 目标误差 epochs = 10000 # 训练次数 delta = 0.0001 # 偏导计算的增量 ov = 0 # 初始化目标值 for k in range(epochs): # 误差计算 tl = l(x, t, y) # 计算残差值 tmp_ov = np.sum(o(tl)) # 计算最新目标值 tmp_error = np.abs(tmp_ov - ov) # 目标值的变化 print(f'{k}->目标值: [{tmp_ov}], 误差: [{tmp_error}]') ov = tmp_ov # 更新目标值 if tmp_error < goal: print(f'训练次数: {k+1}次...') break # 执行训练 tmp_j1 = j1(x, t, y) # 关于第一个参数的一阶偏导函数的值 tmp_j2 = j2(x, t, y) # 关于第二个参数的一阶偏导函数的值 H = np.array([[np.sum(tmp_j1 * tmp_j1), np.sum(tmp_j1 * tmp_j2)], [np.sum(tmp_j2 * tmp_j1), np.sum(tmp_j2 * tmp_j2)]]) # 近似黑塞矩阵 J = np.array([np.sum(tl * tmp_j1), np.sum(tl * tmp_j2)]) # 近似雅克比矩阵 # 计算值λ if k == 0: lambda_val = 0.000001 * np.max(np.diag(H)) v = 2 else: tmp_delta = delta * np.ones(J.shape) beta = (tmp_ov - np.sum(o(l(x + tmp_delta, t, y)))) / (tmp_ov - (tmp_ov + J.T.dot(tmp_delta) + 0.5 * tmp_delta.T.dot(H).dot(tmp_delta))) if beta > 0: lambda_val = lambda_val * max([1/3, 1 - np.power(2 * beta - 1, 3)]) v = 2 else: lambda_val = lambda_val * v v = 2 * v print('-->值λ: ' + str(lambda_val)) # 更新参数 x = x - np.linalg.inv(H + lambda_val * np.eye(H.shape[0])).dot(J) # 打印训练结果 print('训练后的参数: ' + str(x)) plt.plot(t, y, 'o', np.arange(min(t), max(t), 0.001), f(x, np.arange(min(t), max(t), 0.001)), '-') plt.show()
0->目标值: [55574.229250000004], 误差: [55574.229250000004]
-->值λ: 0.006333249999994566
1->目标值: [0.04451695772021884], 误差: [55574.18473304228]
-->值λ: 0.002111083333331522
2->目标值: [0.04445524644085119], 误差: [6.171127936764609e-05]
训练次数: 3次...
训练后的参数: [ 2.99520268 10.03331413]
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。