Python 算例实现Levenberg-Marquardt算法
2018-01-11 13:35
3379 查看
第一次写技术博文,有错误的地方欢迎指点。
本博文是通过一个算例对LM算法的学习进行总结,编程语言是python。
理论就不讲了,网上一大堆。
拟合函数 y(x) = exp(a*x^2 + b * x + c) 中的参数 a, b, c 。废话不多说,直接上代码:
拟合效果如下图所示:
拟合结果如下:
step = 1,abs(mse-lase_mse) = 6427.823881782
step = 2,abs(mse-lase_mse) = 1348.739146779
step = 3,abs(mse-lase_mse) = 7725.330011166
step = 4,abs(mse-lase_mse) = 51.226834643
step = 5,abs(mse-lase_mse) = 0.013664543
step = 6,abs(mse-lase_mse) = 0.000026830
step = 7,abs(mse-lase_mse) = 0.000000194
参数:[[ 0.97256842]
[ 3.03242402]
[ 1.99001528]]
可以看到,算法迭代了七次就收敛了。并且和真实值 [1, 3, 2] 很接近。
本博文参考了下篇博文:
http://blog.csdn.net/jinshengtao/article/details/53310804
本博文是通过一个算例对LM算法的学习进行总结,编程语言是python。
理论就不讲了,网上一大堆。
拟合函数 y(x) = exp(a*x^2 + b * x + c) 中的参数 a, b, c 。废话不多说,直接上代码:
# -*- coding:utf-8 -*- # autor:HuangYuliang import numpy as np from numpy import matrix as mat from matplotlib import pyplot as plt import random n = 100 a1,b1,c1 = 1,3,2 # 这个是需要拟合的函数y(x) 的真实参数 h = np.linspace(0,1,n) # 产生包含噪声的数据 y = [np.exp(a1*i**2+b1*i+c1)+random.gauss(0,4) for i in h] y = mat(y) # 转变为矩阵形式 def Func(abc,iput): # 需要拟合的函数,abc是包含三个参数的一个矩阵[[a],[b],[c]] a = abc[0,0] b = abc[1,0] c = abc[2,0] return np.exp(a*iput**2+b*iput+c) def Deriv(abc,iput,n): # 对函数求偏导 x1 = abc.copy() x2 = abc.copy() x1[n,0] -= 0.000001 x2[n,0] += 0.000001 p1 = Func(x1,iput) p2 = Func(x2,iput) d = (p2-p1)*1.0/(0.000002) return d J = mat(np.zeros((n,3))) #雅克比矩阵 fx = mat(np.zeros((n,1))) # f(x) 100*1 误差 fx_tmp = mat(np.zeros((n,1))) xk = mat([[0.8],[2.7],[1.5]]) # 参数初始化 lase_mse = 0 step = 0 u,v= 1,2 conve = 100 while (conve): mse,mse_tmp = 0,0 step += 1 for i in range(n): fx[i] = Func(xk,h[i]) - y[0,i] # 注意不能写成 y - Func ,否则发散 mse += fx[i,0]**2 for j in range(3): J[i,j] = Deriv(xk,h[i],j) # 数值求导 mse /= n # 范围约束 H = J.T*J + u*np.eye(3) # 3*3 dx = -H.I * J.T*fx # 注意这里有一个负号,和fx = Func - y的符号要对应 xk_tmp = xk.copy() xk_tmp += dx for j in range(n): fx_tmp[i] = Func(xk_tmp,h[i]) - y[0,i] mse_tmp += fx_tmp[i,0]**2 mse_tmp /= n q = (mse - mse_tmp)/((0.5*dx.T*(u*dx - J.T*fx))[0,0]) if q > 0: s = 1.0/3.0 v = 2 mse = mse_tmp xk = xk_tmp temp = 1 - pow(2*q-1,3) if s > temp: u = u*s else: u = u*temp else: u = u*v v = 2*v xk = xk_tmp print "step = %d,abs(mse-lase_mse) = %.8f" %(step,abs(mse-lase_mse)) if abs(mse-lase_mse)<0.000001: break lase_mse = mse # 记录上一个 mse 的位置 conve -= 1 print xk z = [Func(xk,i) for i in h] #用拟合好的参数画图 plt.figure(0) plt.scatter(h,y,s = 4) plt.plot(h,z,'r') plt.show()
拟合效果如下图所示:
拟合结果如下:
step = 1,abs(mse-lase_mse) = 6427.823881782
step = 2,abs(mse-lase_mse) = 1348.739146779
step = 3,abs(mse-lase_mse) = 7725.330011166
step = 4,abs(mse-lase_mse) = 51.226834643
step = 5,abs(mse-lase_mse) = 0.013664543
step = 6,abs(mse-lase_mse) = 0.000026830
step = 7,abs(mse-lase_mse) = 0.000000194
参数:[[ 0.97256842]
[ 3.03242402]
[ 1.99001528]]
可以看到,算法迭代了七次就收敛了。并且和真实值 [1, 3, 2] 很接近。
本博文参考了下篇博文:
http://blog.csdn.net/jinshengtao/article/details/53310804
相关文章推荐
- 基于Levenberg-Marquardt训练算法的BP网络Python实现
- 基于Levenberg-Marquardt训练算法的BP网络Python实现
- 基于Levenberg-Marquardt训练算法的BP网络Python实现
- 基于Levenberg-Marquardt训练算法的BP网络Python实现
- 基于Levenberg-Marquardt训练算法的BP网络Python实现
- LM(Levenberg-Marquard) c语言实现
- LM(Levenberg-Marquard)算法的实现
- Eigen中Levenberg-Marquardt算法的应用
- Levenberg–Marquardt算法
- 相机标定:关于用Levenberg-Marquardt算法在相机标定中应用
- 迭代求解最优化问题——Levenberg-Marquardt算法
- 利用Levenberg_Marquardt算法求解无约束的非线性最小二乘问题~
- 【转】LM(Levenberg-Marquard) Matlab及C语言实现
- LM(Levenberg-Marquard)算法的实现
- 【转】LM(Levenberg-Marquard)算法的实现
- Levenberg-Marquardt算法简介和C++实现
- python实现的json数据以HTTP GET,POST,PUT,DELETE方式页面请求
- python实现虎扑网站图片爬虫
- Python实现简单的udp打洞(P2P)
- python实现 全局变量的两种解决办法