资源简介
Levenberg-Marquardt非线性最小二乘方法在一个二维位置参数拟合问题上的应用
代码片段和文件信息
#include “levernberg.h“
int main（）
{
//拟合数据
float data[9]={0.250.511.523468};
float obs[9]={19.2118.1515.3614.1012.899.327.455.243.01};
//拟合参数的初值
float a0=10;
float b0=0.5;
//y_init = a0*exp（-b0*data_1）;
float a_estb_est;
//a_est=1;b_est=1;
lm（dataobsa0b0&a_est&b_est）;
printf（“%f %f\n“a_estb_est）;
return 0;
}
void lm（float* datafloat* obsfloat a0float b0float* a_estfloat* b_est）//一定要记得函数的形参设置成指针，才能改变最终的值！
{
//数据维数 datasize
int datasize=9;
//未知参数维数 parasize
int parasize=2;
//迭代最大次数
int n_iters=50;
//LM算法的阻尼系数初值
float lamda=0.01;
//step1: 变量赋值
int updateJ=1;
*a_est=a0;
*b_est=b0;//给出a0b0的初值
//step2: 迭代
float* J=new float[datasize*parasize];
float* Jt=new float[parasize*datasize];
float* y_est=new float[datasize];
float* y_est_lm=new float[datasize];
float* d=new float[datasize];
float* d_lm=new float[datasize];
float* H=new float[parasize*parasize];
float* H_lm=new float[parasize*parasize];
float* H1=new float[parasize*parasize];
float* g=new float[parasize*1];
float* dp=new float[parasize*1];
int itij;
float ee_lma_lmb_lm;
for（it=1;it<=n_iters;it++）
{
if（updateJ==1）
{
//根据当前估计值，计算雅克比矩阵
for（i=0;i
*（J+i*parasize+0）=exp（-（*b_est）*（*（data+i）））;
*（J+i*parasize+1）=-（*a_est）*（*（data+i））*exp（-（*b_est）*（*（data+i）））;
}
// J=[exp（-b_est*data_1（i）） -a_est*data_1（i）*exp（-b_est*data_1（i））];
// 根据当前参数，得到函数值
for（i=0;i
*（y_est+i）=*a_est*exp（-（*b_est）*（*（data+i）））;
}
//计算误差
for（i=0;i
*（d+i）=*（obs+i）-（*（y_est+i））;
}
//计算（拟）海塞矩阵
transpose（JdatasizeparasizeJt）;
multiplication（JtparasizedatasizeparasizeJH）;
//若是第一次迭代，计算误差
if（it==1）
{
e=dot（ddatasized）;
}
}
//根据阻尼系数lamda混合得到H矩阵
for（i=0;i
for （j=0;j
if （i!=j）
{
*（H_lm+i*parasize+j）=*（H+i*parasize+j）;
}
else
*（H_lm+i*parasize+j）=*（H+i*parasize+j）+lamda;
}
}
inv2（H_lmH1）;
multiplication（Jtparasizedatasize1dg）;
multiplication（H1parasizeparasize1gdp）;
//依据步长更新ab
a_lm=*a_est+（*dp）;
b_lm=*b_est+*（dp+1）;
//计算新的可能估计值对应的y和计算残差e
for（i=0;i
*（y_est_lm+i）=a_lm*exp（-b_lm*（*（data+i）））;
}
for（i=0;i
*（d_lm+i）=*（obs+i）-（*（y_est_lm+i））;
}
e_lm=dot（d_lmdatasized_lm）;
//根据误差，决定如何更新参数和阻尼系数
if （e_lm
lamda=lamda/10;
*a_est=a_lm;
*b_est=b_lm;
e=e_lm;
updateJ=1;
}
else
{
updateJ=0;
lamda=lamda*10;
}
}
delete [] J;
J=NULL;
delete [] Jt;
Jt=NULL;
delete [] y_est;
y_est=NULL;
delete [] y_est_lm;
y_est_lm=NULL;
delete [] d;
d=NULL;
delete [] d_lm;
d_lm=NULL;
delete [] H;
H=NULL;
delete [] H_lm;
H_lm=NULL;
delete [] H1;
H1=NULL;
delete [] g;
g=NULL;
delete [] dp;
dp=NULL;
}属性 大小 日期 时间 名称
----------- --------- ---------- ----- ----
文件 74 2014-02-12 13:34 MyLM\1.txt
文件 33280 2014-02-13 14:37 MyLM\Debug\MyLM.exe
文件 408264 2014-02-13 14:37 MyLM\Debug\MyLM.ilk
文件 617472 2014-02-13 14:37 MyLM\Debug\MyLM.pdb
文件 1174 2014-02-13 14:37 MyLM\MyLM\Debug\cl.command.1.tlog
文件 8470 2014-02-13 14:37 MyLM\MyLM\Debug\CL.read.1.tlog
文件 460 2014-02-13 14:37 MyLM\MyLM\Debug\CL.write.1.tlog
文件 39309 2014-02-13 14:37 MyLM\MyLM\Debug\levernberg.obj
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
文件 2 2014-02-13 14:37 MyLM\MyLM\Debug\li
............此处省略32个文件信息
川公网安备 51152502000135号
评论
共有 条评论