#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\link-cvtres.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link-cvtres.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.14048-cvtres.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.14048-cvtres.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.14048.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.14048.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.2232-cvtres.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.2232-cvtres.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.2232.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.2232.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5560-cvtres.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5560-cvtres.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5560.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5560.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5616-cvtres.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5616-cvtres.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5616.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5616.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5672-cvtres.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5672-cvtres.write.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5672.read.1.tlog

?????文件??????????2??2014-02-13?14:37??MyLM\MyLM\Debug\link.5672.write.1.tlog

............此處省略32個文件信息