多峰曲线拟合

#/13step5：若k=n，令x=x，转step2.0kstep6：d=—g+dstep6：kkk—1step7:若果dTg>0，令x=x，转step2；否则转step3。kk0k在第三步中，要求米用一维搜索获取最佳步长九，本文米用0.618法来作精k确线性搜索。2.1.20.618法0.618法是一种用于求单峰函数的极小值的分割法。其基本思想是通过取试探点和进行函数值的比较，使包含极小点的搜索区间不断缩短，当区间长度缩短到一定程度时，区间上各点的函数值均接近极小值，从而各点可以看作为极小点的近似。这种方法仅需计算函数值，不涉及导数，又称直接法。下面介绍0.618法的技术细节。设包含极小点的初始搜索区间为［a,b］，设^(a)=f(x+ad)在［a,b］上kk是凸函数。0.618法的基本思想是在搜索区间［a,b］上选取两个对称点a，a，12其中a是区间［a,b］上的黄金分割点，a是a的对称点，且a<a，比较这两21212点处的函数值申(a)和申(a)的大小，若申(a)>申(a)，则说明极小值在区间1212［a,b］内，故取新区间为［a,b］，若申(a)<^(a)，则说明极小值在区间［a,a］11122内，故取新区间为［a,a］。新区间的长度为原区间长度的0.618倍。新区间包含2原区间中两个对称点中的一点，我们只要再选一个对称点，并利用这两个新对称点处的函数值继续比较，重复这个过程，最后即可确定极小点a*。下面给出具体算法步骤。stepl:选取初始数据。确定初始搜索区间［a,b］和精度要求8>0。计算最11初两个试探点九，卩，11九=a+0.328(b—a)1111卩二a+0.618（b—a）1111计算甲（九）和申（卩），令k=1。11step2:比较目标函数值。若申（九）>申（卩），转step3，否则转step4。kkstep3：若b-X<5，则停止计算，输出卩；否则，令TOC\o"1-5"\h\zkkka=X,b=b,X=卩，k+1kk+1kk+1k申（X）=申（卩）,卩=a+0.618（b-a）k+1kk+1k+1k+1k+1计算甲（卩），转step2。k+1step4：若卩-a<5，则停止计算，输出X，否则，令：kkka=a,b=卩,卩=Xk+1kk+1kk+1k申（卩）=申（X）,X=a+0.382（b—a）k+1kk+1k+1k+1k+1计算申（X），转step2。k+12.1.3计算效果本文采用C#作为编程语言来实现上述算法。当进行多次调试之后，发现共轭梯度法并不能解决问题。现总结原因如下。1、迭代初始值难以确定。共轭梯度法在开始之前要求设定初始迭代值。该方法对初始值的设定非常敏感，如果初始值设置的不够好，这种方法的收敛将会变得非常慢，甚至不能收敛。由于每次进行电泳的DNA片段有可能不同，峰值信息业会不断变化，因此没有一个好的方法来为算法提供初始值。2、在一维线性搜索方法中，首先要求给出搜索区间。这个搜索区间也是由于问题的实际背景决定的。在本文中，不能确定峰值信息的取值区间，因此一维搜索只能采用纯数学意义上的搜索区间，这样搜索负担将会变得非常沉重，迭代计算变得非常耗时。

经过再三权衡，笔者决定放弃这种方法。可能笔者对共轭梯度法的认识还不够，也可能这种方法不适合本问题的求解。因此笔者决定采用《数值分析》中的多项式拟合算法来获取峰值信息。2.2分段多项式拟合2.2.1多项式拟合的不足多项式拟合算法简单，易于实现，是工程应用中常用的拟合方法。多项式拟合对于单峰、变化平缓的曲线拟合效果较好，但对于多峰曲线拟合来说，则差强人意。图4是对电泳数据作40阶多项式拟合的结果。从图4可以看到，多项式函数不能表现曲线的突变趋势，因而便不能直接用来拟合多峰曲线。标題I一廿生茁1I一廿生茁1050010001500200025003M0timififEra图4多项式拟合电泳曲线为了改变这种现状，本文采用分段多项式拟合技术，实践证明，将数据人为分段之后，对每一段单独进行拟合，效果很好［3］。分段的目的就是要把曲线平坦的部分和有峰的部分区分开，这样才会得到好的拟合效果。下面介绍针对于本问题的一种分段方式。2.2.2分段标准对于整段曲线来说，有峰的地方所占的比例不大，因此整段曲线的平均值将接近于背景暗电流值。在经过试验以后，发现曲线的平均值如图5所示。在图5

中，红线表示了整段曲线的平均值。可以看到，这条红线并不能“经过”所有的峰，这样就不能够获取分组依据。datal标题44—I111datal标题44—I111111111111111111111111111105001000150020002M03D0<]time图5曲线平均值示意图考虑到去掉一些暗电流值，再将余下的数据作平均，则可能使得新的平均值“经过”所有的峰。本文采用这样一种策略：以第一次求得整条曲线的平均值作为标准，将每个数据与之比较，去掉所有小于平均值的数据点，对剩下的数据点再做一次平均，将新获得的平均值作为分组依据。图6显示了新平均值的位置。可以看到，新的平均值经过所有峰，这样就可以根据新的平均值与原曲线的“交点”来划分曲线了。标题22-J"1"11J■1"1■J1■■1■11!■!■11■1!]20-图6新平均值示意图对每一段曲线作40阶多项式插值，最后得到的插值效果如图7所示。由于

多项式插值的算法已经非常成熟，在此就不作详细介绍[4]。需要指明的是，本文

在实现多项式插值算法的时候，求解法方程组所采用的算法是Gauss按比例列主

元消去法。图8、图9为图7采取局部放大之后的图像，来显示多项式插值的细

节。标题2222AHI-n—1图7分段多项式拟合效果曷丁由1-4001500曷丁由1-4001500图8拟合细节dotal10.610.4145015001550time/rns图9拟合细节2.2.3峰位、峰面积计算在分段多项式插值完成后，每一个峰都可以用一个多项式函数来表达，那么对峰位的求解就转化为求这个多项式函数在闭区间内的最小值，对峰面积的求解就转化为对此多项式函数在闭区间内求定积分，这在算法的实现上都非常简单，使得程序设计变得非常轻松。图10下方文本框中显示了峰位、峰面积计算结果。标盏I加富伽~~I050010001500200025003000Eime/ms图10峰位、峰面积计算结果3软件介绍本软件采用微软的WPF程序设计框架设计而成。这使得界面设计变得灵活简单，界面元素也变得更加漂亮。图像显示控件采用一款开源图表控件ZedGraph这款控件编程简单，性能优良。整个软件界面分成三个部分。最上面是图像显示部分，用来显示曲线。左下角为操作部分，有两个按钮，第一个用来载入数据，数据文件格式为.txt格式，第二个按钮用来执行拟合操作。右下角有一个文本框，用来显示计算结果。具体构造可以参看图10。4结语本文为了获取曲线的峰值信息，尝试了多种不同的方式进行求解。实践证明，复杂的算法并不一定能够产生好的效果。就像本文之前介绍的共轭梯度法，用这个方法可以解决大型非线性最优化问题，当然也应该能够解决高斯拟合参数求解问题，但是在实际操作中由于此种方法对于初始迭代值的要求过高，没有能够成功实现，并且该算法计算复杂，耗时很长，对于一个拟合问题来说，的确有大材小用之嫌。所以，本文随后采用了熟悉的方法：多项式拟合。只是在拟合之前对数据作了分段处理。这样一来，可以很完美的解决问题，而算法设计又非常简单，且计算非常轻松。所以，在考虑问题解决方式的时候，对于同样能解决问题的方法，应当选择那些简单而易于实现的方法，这样才能得到最快的完成效率。参