用matlab拟和模型参数和计算参数误差

Matlab用以建立数学模型是一个很好的工具。对模型函数的评价，一个很重要的方法就是最小二乘(Least squares)由least mean squares这个方法得到。假如有点集P(X, Y)，每一个点 P(i) 由X(i), Y(i) , i = 1 ~ m组成；模型 Y_fit = F( A, X ), Y_fit(i) = F(A, X(i) ); 其中 A= A(1) A(2) … A(n)是模型的n个参数。least mean squares = (1/m) * sum ((Y(i) - Y_fit(i) ).^2)? ?(i = 1 ~ m)。一个好的模型，least mean squares就小；而另一方面，如何得到模型参数A，使得least mean squares有最小值，就是所谓的，最小二乘拟合(least squares curve fitting)了。简介：模型有线性和非线性之分。对于线性模型，求参数，其实就是求一步矩阵的逆(稍候我们可以看到)。而非线性模型，往往不能一步就得到结果，所以就需要多步逼近。就这样，在众多的多步逼近的方法中，最快收敛于最佳参数值的方法就比较垂青。这中间，最强的当然就是Newton 法：A: n+1 = A: n + (Hessen ( L ))^-1??*??grad(L)这里Hessen ( L )是被拟合的模型函数的least mean squares方法的Hessen矩阵。grad(L)是她的梯度矩阵。参数矩阵A的当前值是A:n和下一步值A: n+1。这个方法包含了一个求hessen矩阵的逆的运算。其实，这个方法难的不是这个逆，而是如何得到Hessen矩阵和梯度矩阵。梯度矩阵还好说，就是least mean squares方法的对各个参数的一介偏导数。而Hessen矩阵包含了一介偏导数的组合(主要是相乘)，和二介偏导数。当然，许多模型的二介偏导数相对于一介偏导数的组合是一个比较小的量，特别是线性模型，就没有二介偏导(所以，线性模型可以直接求出参数)。于是，新的方法就利用这个特点，将逼近限制在一介偏导数构成的伪Hessen矩阵上。这就诞生了两个比较著名的方法Gauss-Newton 法和Levenberg-Marquardt法。Gauss-Newton 法直接用Jacobian 行列式代替 Hessen矩阵，用least squares值代替梯度(注意，不是least mean squares，因为当用Jacobian 行列式代替Hessen矩阵时，中间有一个自由度的差别)这里的拟合就变成了A: n+1 = A: n + (Jacobian ( L ))^-1??*??L? ?? ?? ? (对L的定义会在下文中给出)因为越是接近最佳值(或者临界值)，Jacobian ( L )就越是畸形，所以在实际的计算机运算中，求逆这一步都用所谓的帽子运算符 (假如 J= Jacobian ( L ) )；( Jacobian ( L ) )^-1? ?--->? ?( ( Trans(J) * J )^-1) * Trans(J)这里Trans()是转置运算。Levenberg-Marquardt法比Gauss-Newton好，可以说，在很多时候，和Newton 法是一样的(当参数互相独立的时候，而这个条件是好的数学模型基本的条件，虽然有时候不能100%保证，但是，基本上参数互相之间相关性不大)。与Newton 法相比，Levenberg-Marquardt法因为没有纯的二介偏导运算，速度有时候还比Newton 法要快。而且，Levenberg-Marquardt法引入了一个步长，使得该方法在计算机运算上，更容易控制。实现：下面我们简单说一下这个Levenberg-Marquardt法，和拟合过程中的参数误差估计：假设 L = sum ((Y - Y_fit ).^2)? ?(Y= [Y (1), Y(2)… Y(m)]是测量值，Y_fit = F( A, X )=[ F( A(1), A(2)…A(n), X(1)), F( A(1), A(2)…A(n), X(2))…F( A(1), A(2)…A(n), X(m))] ), 是模型值。其中的X(1)~(m)就是对应于Y的测量值，A(1)~(n)是模型的参数)使L最小的参数A就是我们要求的参数。当L线性时，L的最小值，就是要求满足L对于各个参数的偏导数=0的点(相当于解一个多元一次方程组)。当L非线性时，上述的条件依然有效，就是不能解。于是就有了逼近法。考虑到非线性模型的二介偏导数不全为零，如果给定一个接近最佳参数值的估