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

recbio

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非线性时，上述的条件依然有效，就是不能解。于是就有了逼近法。考虑到非线性模型的二介偏导数不全为零，如果给定一个接近最佳参数值的估计参数A:current，最佳参数（使得L最小的参数）A:min可由Newton 法给出：
A:min = A:current + Hessen (L)^-1 * grad(L)。
grad(L)是一个一维向量，第j个向量g(j)是L对参数A(j)的一介偏导数，所以，有n个元素。
假设F( A, X )对A(j)的一介偏导数为D_Y_fit(j)，这样：
g(j)= - 2 * sum[(Y - Y_fit)* D_Y_fit(j)]
Hessen (L)是一个二维的矩阵，第i行第j列的值
H(i, j)=2*sum[D_Y_fit(i)* D_Y_fit(j)- (Y - Y_fit)* DD_Y_fit(j, i)]
这里DD_Y_fit(j, i)是模型函数F( A, X )先对A(j)参数求一介偏导再对A(i)参数求第二介偏导。

在Levenberg-Marquardt法中，首先考虑到(Y - Y_fit)* DD_Y_fit(j, i)这一项，在接近最佳值的时候，(Y - Y_fit)很小；其次，假如参数的互相之间相关性很小，这就使得DD_Y_fit(j, i)在i, j不相等的时候值很小甚至等于0（参数之间相互正交）。再次，考虑到2和-2是常数，可以互相抵消。于是，LM法用H’(i, j)（几何值相当于Hessen矩阵的一半）来代替Hessen矩阵，具体的计算方法如下：
H’(i, j)= sum[D_Y_fit(i)* D_Y_fit(j)] (i不等于j 时)
H’(i, i)= (1 + r) * sum[D_Y_fit(i)* D_Y_fit(i)]
假设D_A = ( A:min - A:current)是最佳参数和当前估计参数的差，是一个一维向量。
LM法就是
1) 假设一个r比如可以从1开始，解方程 D_A * H’ = (-1/2)grad(L)
细心的同学已经发现了，这是一个线性方程，意味着在matlab里面一步就可以解决：
D_A = H’ \ ((-1/2)grad(L))
没有错，但是！（见后面）
2) 比较L(A:current + D_A ) 和L(A:current )，这个时候假如
a) L(A:current + D_A ) < L(A:current )，就保留A:current + D_A，并减小r（比如r=r/10）
b) 否则就增加 r（比如r=r*10）
3) 重复解方程 D_A * H’ = (-1/2)grad(L)
这个方法可以有两个结束指标a) r很大时，说明方程到了平坦区。b) 每一次解方程以后，会比较一次L，假如L还在减小，但是减小的幅度很小时，可以结束。
实际应用中，没有用r很小来结束，因为L变化的减小比r快。当然，r很小也是一个结束条件，只是，r还没有很小时，L的变化已经很小了。说明拟合已经十分接近终点了。

有一点要注意，就是接近终点时计算Hessen矩阵（Newton 法）或者H’（ Levenberg-Marquardt法）的inv会变得十分不稳定，就算matlab说用QR decomposition，也可能出现错误。如果你用r=1e+20计算时，基本上matlab backslash 是得不到什么结果的。建议用函数inv()来计算。或者用经过“最大绝对值优先”的Gauss-Jordan法。他们都比较稳定。

关于拟合的误差：
最小二乘的拟合误差可以用1/2Hessen矩阵来估计，Levenberg-Marquardt法中可以用r=0的时候的H’的inv来计算。具体的：
A的无偏标准误差（覆盖68%的置信区间）= +/- sqrt( L/(自由度)* diag(H’^-1))
这里，L就是用最后得到的拟合参数计算的方差和，自由度的定义是
自由度 = 拟合样本数(x-y 的个数，m) - 拟合参数个数 (n)

总结：
matlab计算非线性最小二乘的方法大致就是这样一个过程：
1) 确定拟合函数F( A, X )，写出函数对各个参数的一介偏导数的函数D_Y_fit(1~n)（如果手算有困难，可以用Mathematica之类的支持符号运算的软件帮忙，据说也可以用matlab的符号运算工具箱，不过自己没有试过），如果偏导数比较简单，还可以再写出二介偏导数DD_Y_fit(1~n, 1~n)（共n*n个，用于Newton 法）。
2) 用一组估计值作为A:current，计算L = sum ((Y - F( A, X ) ).^2)  (least square)
3) 根据1)得到的函数，用上面的定义，用估计值A:current构建(-1/2)grad(L) 矩阵G。第j个元素g'(j)可以由下面的方法得到：
g'(j)= sum[(Y - Y_fit)* D_Y_fit(j)]
和H’矩阵（Levenberg-Marquardt法）：
H’(i, j)= sum[D_Y_fit(i)* D_Y_fit(j)] (i不等于j 时)
H’(i, i)= (1 + r) * sum[D_Y_fit(i)* D_Y_fit(i)]
或者1/2Hessen矩阵（Newton 法）
H(i, j)= sum[D_Y_fit(i)* D_Y_fit(j)- (Y - Y_fit)* DD_Y_fit(j, i)]
因为A有了，X Y有了，所以上面的矩阵都是数字矩阵。对于LM法，r可以从1开始。
4)解方程
D_A= H’ \ G (LM法) 或者 D_A= H \ G (Newton 法)
用A:next = A:current + D_A 重新计算2)
注意，计算2)之前，可以对A:next作一个范围检查，对于越过边界条件的值可以修正。方法：i ) 保留原来的值；ii )用边界条件；iii)对于第j个不合条件的值，取一个小的D_A(j)（比如A:next(j) = A:current(i) + D_A(i)/10 , 重复直到满足条件）。
5) 如果是LM法，根据上面文章，改变r值，重复3)
a) L(A:next) < L(A:current )，就保留A:next (A:current = A:next)，并减小r（比如r=r/10）
b) 否则就增加 r（比如r=r*10）放弃A:next

如果是Newton 法
a) L(A:next) < L(A:current )，就保留A:next (A:current = A:next)，
b) 否则，就结束！说明A:current已经是最佳值了。
通常，Newton 法结束的时候，D_A会等于0。


6) 在重复3)之前，可以用文章提到的结束条件判断一下。
7) 完全结束了可以根据需要，求一个误差。不要的话，就结束好了。

用matlab大致就是这样。

其实写这篇文章的时候，自己正好在做关于非线性拟合的模型。发现matlab的fit不能给出参数的误差，就用以前的c++的方法改了一下，放在matlab里面。顺便总结了上面的文字，当然，只是自己的理解，如果有什么不完整的，还希望大家讨论。

下次我会帖c++的非线性拟和的方法，主要是Levenberg-Marquardt方法。大家如果有比较棘手的函数需要拟和，不妨帖一个，看看是否可以作为例子，用c++解决。

有任何问题，也欢迎大家回帖或者给自己消息！

freecat

想问一下
确定拟合函数F( A, X )这一步
如果函数非线性而且非常复杂，而且不能用一个函数来描述，也就是说F未知
这种情况，用什么方法求极值

[ 本帖最后由 freecat 于 2007-6-22 23:42 编辑 ]

recbio

这个问题非常好。其实自己也不是很清楚这个中间的一般方法。

通常的，如果没有规律，且非线性。如果是拟合某一个连续的区间，多项式拟合总是可以的(类似原函数的泰勒展开)。就是参数的意义不十分明确了。

如果是不连续的，那么就只有自己建立拟合规则了。

实在不行，还可以用"枚举"，建立评价函数，利用计算机的速度来解决。如果实在太多了，就用Monte Carlo法，然后看运气了。
一般来说，不能100%的解决，还是可以在一个比较高的概率范围内解题的。

当然，具体的方案还需视实际情况而定。很高兴进一步讨论。

狐狸老大

建立评价函数，后用启发式的迭代优化方法求极值，是一个行之有效的方法，虽然这种方法的收敛性质未必一定稳定，但确实可以在一个比较高的概率范围内解题的。

freecat

原帖由 recbio 于 2007-6-23 22:13 发表
这个问题非常好。其实自己也不是很清楚这个中间的一般方法。

通常的，如果没有规律，且非线性。如果是拟合某一个连续的区间，多项式拟合总是可以的(类似原函数的泰勒展开)。就是参数的意义不十分明确了。

...

我自己的论文也是卡在了这里，对于一个所谓评价函数的求解，倒是可以用所谓启发式迭代求近似最优解，但是在函数和参数的一一对应关系不明确之前，用传统的梯度迭代法却遇到了困难。

qqsweet

LZ是个热心人$送花$