Module VarLMImplicit none!//Nparas為參數個數，Ndata為數據個數,Niters為最大迭代次數Integer ,Parameter :: Nparas = 2,Ndata = 9,Niters = 50,Fileid = 11!//臨時變量，order=1，繼續迭代，否則停止迭代Integer :: order!//循環變量，it為迭代次數循環變量Integer :: it,i,ii!//x_1為觀測數據自變量，y_1為觀測數據因變量，兩者均為已知數據Real(kind=8),Parameter :: x_1(Ndata) = [0.25,0.5,1.0,1.5,2.0,3.0,4.0,6.0,8.0]Real(kind=8),Parameter :: y_1(Ndata) = [19.306,18.182,16.126,14.302,12.685,9.978,7.849,4.857,3.005] !//a0與b0分別為參數猜測初始值,RealA、RealB為真值Real(kind=8),Parameter :: a0 = 0.,b0 = 0.,RealA = 20.5,RealB = 0.24!//誤差限度Real(kind=4) :: eps = 1e-8!//y_est為估計值，d為估計值和實際值y_1之間的誤差，rd為數組d的變形Real(kind=8) :: y_est(Ndata) = 0.,d(Ndata) = 0.,rd(Ndata,1) = 0.!//J為雅可比矩陣，JT為J的轉置矩陣，H為海森矩陣Real(kind=8) :: J(Ndata,Nparas) = 0.,JT(Nparas,Ndata) = 0.,H(Nparas,Nparas) = 0.!//Inv_H_1m為H_1m的逆矩陣Real(kind=8) :: H_Lm(Nparas,Nparas)=0.,Inv_H_Lm(Nparas,Nparas)=0.!//eye為單位矩陣，delta為步長矩陣Real(kind=8) :: eye(Nparas,Nparas) = 0.0,delta(Nparas,1) = 0.0!//a_est，b_est分別為反演參數Real(kind=8) :: a_est = 0.,b_est = 0.,e = 0.!//計算新對應的值Real(kind=8) :: y_est_Lm(Ndata) = 0.,d_Lm(Ndata) = 0.,e_Lm = 0.!//臨時變量，lamda為LM算法的阻尼系數初始值Real(kind=8) :: a_Lm,b_Lm,lamda = 0.01,v = 10.d0,tempEnd Module VarLMProgram LMInclude 'link_fnl_shared.h'Use VarLMUse LINRG_INTImplicit None!//生成單位矩陣Forall( i = 1:Nparas, ii = 1:Nparas, i==ii ) eye(i,ii) = 1.d0!//第一步：變量賦值order = 1a_est = a0b_est = b0!//第二步：迭代Write(*,"('------------------------------------------------')") loop1: Do it = 1,NitersIf ( order==1 ) Then !//根據當前估計值，計算雅可比矩陣y_est = a_est*Exp( -b_est*x_1 ) !//根據當前a_est，b_est及x_1，得到函數值y_estd = y_1 - y_est !//計算已知值y_1與y_est的誤差loop2: Do i = 1,Ndata !//計算雅可比矩陣。dy/da = Exp(-b*x)，dy/db = -a*x*Exp(-b*x)J(i,1) = Exp( -b_est*x_1(i) ) J(i,2) = -a_est*x_1(i)*Exp( -b_est*x_1(i) )End Do loop2 JT=Transpose(J) H = Matmul( JT,J ) !//計算海森矩陣End IfIf ( it==1 ) e = Dot_Product(d,d) !//若是第一次迭代，計算誤差epsilonH_Lm = H + lamda*eye !//根據阻尼系數lamda混合得到H矩陣!// 縧amda*I縇evenberg靠縧amda*diag(A'A)縈arquardt靠!//計算步長delta，并根據步長計算新的參數估計值Call LINRG( H_Lm,Inv_H_Lm ) !//使用imsl函數庫，計算H_Lm的逆矩陣Inv_H_Lmrd = Reshape( d,[Ndata,1] ) !//為了滿足內部函數Matmul的計算法則，對d的數組形狀進行改變delta = Matmul( Inv_H_Lm,matmul( JT,rd ) ) !//delta為增量a_Lm = a_est + delta(1,1)b_Lm = b_est + delta(2,1)!//如果||delta||<1e-8，終止迭代If ( Dot_Product(delta(:,1),delta(:,1))<eps ) Exit!//計算新的可能估計值對應的y和計算殘差ey_est_Lm = a_Lm*Exp( -b_Lm*x_1 )d_Lm = y_1 - y_est_Lme_Lm = Dot_Product( d_Lm,d_Lm ) !//e_LM等于||y_1 - y_est_LM||!//根據誤差，決定如何更新參數和阻尼系數!//迭代成功時將lamda減小，否則增大lamdaIf ( e_Lm<e ) Thenlamda = lamda/va_est = a_Lmb_est = b_Lme = e_Lmorder = 1Elseorder = 0lamda = lamda*vEnd IfWrite(*,"('a_est=',g0,2X,'b_est=',g0)") a_est,b_estEnd Do loop1Write(*,"('------------------------------------------------')") Write(*,"('停止迭代，總共迭代',g0,'次')") it - 1!//輸出正反演數據以及百分比誤差輸出到文件,總共100個數據點,計算區域為[0-50]Open(Fileid,file='原始數據與反演數據.dat',status='unknown')Do i = 0,100temp = i/2.d0Write(Fileid,'(f7.3,3f15.8)') temp,RealA*Exp( -RealB*temp ),a_est*Exp( -b_est*temp ),&Abs( a_est*Exp( -b_est*temp )-RealA*Exp( -RealB*temp ) )/RealA/Exp( -RealB*ii )*100End DoClose(Fileid)!//輸出反演參數a_est，b_est以及原始數據和擬合數據Write(*,"(/,'反演參數為：')")Write(*,"('a_est=',g0)") a_estWrite(*,"('b_est=',g0)") b_estWrite(*,"(/,'原始數據為：')") Write(*,'(3f9.4/)') y_1Write(*,"('擬合數據為：')") Write(*,'(3f9.4/)') a_est*Exp( -b_est*x_1 )Write(*,"('------------------------------------------------')") End Program LM

?

總結

以上是生活随笔為你收集整理的LM算法求解最小二乘问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。