曲线曲面拟合.课件.ppt

1、非线性回归分析非线性回归分析（曲线曲面拟合）顾世梁2011年11月1 非线性回归分析的任务 非线性关系是最普遍的变数间量化关系，合适的非线性回归分析对研明变数间的数量关系有重要作用。非线性回归分析的广泛应用，将促使试验研究从定性向定量发展，由粗放向精细发展。线性关系形式单一，而非线性关系多种多样，选择合适的非线性模型并非易事。多项式也是一种（简单的一种）非线性关系，先前已有论述，本章仅讨论多项式以外的纯非线性关系。对于纯非线性回归分析，非线性回归统计数的估计、假设测验等均有很大难度。bx yaeb yaxabx yx2exp() yaxbxcx1exp()K yabx(1exp()mK yab

2、x() yf x,bf2222225 1217 212() (1)(1) 23512134 1126(1)()bxxxbxxxY b bx eb bxxx ebe 231234235671bb xb xb xYb xb xb x非线性回归分析的主要任务有下列4项：1) 建立合适的非线性模型；2) 估计非线性方程的统计数曲线曲面拟合；3) 合理的显著性测验；4) 方程的进一步利用（插值与外推）。2 非线性回归方程的选择主要有3种方法：1）解微分和偏微分方程组）解微分和偏微分方程组1exp()kyabxdyGRydxyyyGR/ /dsolve(Dy+y+c,)y=dsolve(Dy-b*y+c*

3、y2,y(0)=k/(1+a)syms c b k; y=subs(y,c b/k); pretty(y)cybRGR2/, 0yybcyybycy()yf x,b2）根据机理或基本数量关系推导）根据机理或基本数量关系推导 每一种函数关系都有一些基本特点，可以根据这些基本要素确定不同的方程。这些基本要素如零点(初值点)、峰值点（极大、极小）、拐点、渐近点等，应符合数据事实。byaxbxyae 3) 试算、比较与选择试算、比较与选择 当变数间的可能关系所知甚少，可对不同方程进行试拟合，比较分析后选出最佳关系模型。除了前述的关键点数据应与曲线、曲面有好的吻合外，也应保证数据在前、中、后段都能较好地

4、拟合；另外也应保证较高的拟合度（决定系数）、较小的离回归平方和以及较好的插值和外推。 通常，较少参数的曲线刚性有余、柔性不足，而参数较多的方程有较大的柔性。但参数太多往往会过参数化（over-parameterization），拟合的难度大大增加。3 参数估计2211()(| )nniiiiiiQYYYf X b目标函数： 当给定当给定Xi 与与 Yi (i=1,2,n)时，时，Q 也是也是b的的函数：函数： Q=F(b)。 拟合即为寻找拟合即为寻找opt=min(F(b)的过程。发的过程。发展稳定高效实现全局最优拟合的算法是非展稳定高效实现全局最优拟合的算法是非线性回归的关键，难度较大。线性

5、回归的关键，难度较大。1）线性化法）线性化法 对一些简单的方程，我们可以采用数据转换的对一些简单的方程，我们可以采用数据转换的方式将其化成线性方程，然后用一元或多元线性方式将其化成线性方程，然后用一元或多元线性回归的方式进行分析。如：回归的方式进行分析。如：221122121222112112Y = a+b X +b X +b X +b X +bX X lnlnbXYaeYabX，bXaY 1122334455Y =a+b X +b X +b X +b X +b X 其缺陷是该类方法仅适用于简单的方程，而绝大其缺陷是该类方法仅适用于简单的方程，而绝大多数纯非线性方程较复杂，不能用线性化方法进行

6、多数纯非线性方程较复杂，不能用线性化方法进行参数估计。参数估计。2）一些通用方法）一些通用方法(1) 梯度法(快速登山法, Gradient)；211()2()nniiQfYYYfbbb 2211()()nniiiiiiQYYYf11112212()02()0.2()0nininippQfYfbbQfYfQbbbQfYfbb 给定某一起始参数点：给定某一起始参数点： (0)(0)(0)(0)12(, , ,)jpbbbbb(0)若若=0, bj 在该点前后的变化不会使在该点前后的变化不会使Q变化变化0, bj 在该点的增加将使在该点的增加将使Q变大变大令令1()njijfYfbjQb00j0j

7、朝着使朝着使Q减小的方向减小的方向 因而因而(1)(0)(0)+jjjbbkjQb12()niQfYfbb ( 1)( )( )+lllbbk一个实例：b0=3，20，0.511exp()fKa-bX2exp()(1exp()f-K-bXaa-bX1exp()KfYa- bX2exp()(1exp()fKaX-bXba-bXXYf0Y-f0df/dKdf/dadf/db20.300.35896-0.058960.11965-0.01580.6320140.860.809340.050660.26978-0.029552.3639961.731.503200.226800.50107-0.037

8、54.4999882.202.195690.004310.73190-0.029434.70937102.472.64373-0.173730.88124-0.01573.13958122.672.85830-0.188300.95277-0.006751.62009142.802.94627-0.146270.98209-0.002640.738790()ifY - fb-.35275-4.813e-3.16480(2) 高斯法(Gauss)；(3) 高斯-牛顿法(Gauss-Newton)；( 1)( )( )+lllbbk(1)(0)(0)+bbk(1)3.0-.352752.96472

9、0 +.1 -4.813e-319.99950.5.1648050.51648Kab 以新的以新的b值再运行前述过程，反复迭代，值再运行前述过程，反复迭代，直至直至delta趋于趋于0，或，或Q已不再变小。已不再变小。(, )Yf X bf(0)b f 按多元按多元Taylor级数展开（略去二次及二次以上各项）：级数展开（略去二次及二次以上各项）： 0000i01212(, )iiiiiiijpjpffffYf X bffbbbb则目标函数可以转化为则目标函数可以转化为:220000121112(; )()nniiiiiiipiipfffQYf X bYfbbb(0)jjjbb00000001

10、201111122()nnnniiiiiiipiiiiiijjjpjjfffffffQYfbbbbbbbb001niijkijkffabb001()nijYiiijfaYfb令11 112 21121 122 2221 12 2p pYp pYpppp ppYaaaaaaaaaaaaA=K-1=A K得新的优化点得新的优化点: :k(0)b=b + 当当b与与b(0)有差异时，有差异时，应令应令b替代替代b(0)重新计算重新计算 (0)jjjbb由由 =0，或，或Q的前后的前后差异小于某一定值。差异小于某一定值。一个实例：b0=3，20，0.511exp()fKa-bX2exp()(1exp(

11、)fK-bXaa-bX1exp()KfYa- bX2exp()(1exp()fKaX-bXba-bXXYf0Y-f0df/dKdf/dadf/db20.300.35896-0.058960.11965-0.01580.6320140.860.809340.050660.26978-0.029552.3639961.731.503200.226800.50107-0.03754.4999882.202.195690.004310.73190-0.029434.70937102.472.64373-0.173730.88124-0.01573.13958122.672.85830-0.188300

12、.95277-0.006751.62009142.802.94627-0.146270.98209-0.002640.73879df/dKdf/dadf/dbY-f00.1197-0.01580.63201-0.058960.2698-0.02962.363990.050660.5011-0.03754.499980.22680.7319-0.02944.709370.004310.8812-0.01573.13958-0.173730.9528-0.00681.62009-0.18830.9821-0.00260.73879-0.14627XY, A = XX, K = XY = A K得新

13、的优化点得新的优化点: :( )1( 1)lllk()b =b+ 反复迭代反复迭代(4) 改良高斯牛顿法(Levenberg-Marquardt) 这是梯度法和高斯-牛顿法相结合的一种方法 。( ,1,2, )jkaj kpppA 很可能是奇异的，很可能是奇异的，需对此阵进行调整需对此阵进行调整: :*diag*A =A+(A)作用: 一可解决A阵奇异，无法求解之困；二是A阵对角线元素包含了较大的与求解相关的信息量，加快趋于全局优的进程。 (5) 极大似然法(maximum likelihood)。 大多数著名的统计软件如 SAS, Matlab, Sigmaplot等包含了基于这些算法的非线

14、性方程拟合模块。3）上述通用算法存在的问题：）上述通用算法存在的问题：(1) 需提供方程的导数或偏导数；(2) 需提供合适的初值；(3) 一般难于实现全局最优拟合。最后一点往往是最主要、最致命的缺陷。4）曲线、曲面拟合新算法）曲线、曲面拟合新算法 (Contraction-Expansion Algorithm) CE算法包含三个基本步骤：(1) 收缩步，缩小步长的搜索过程；收缩步，缩小步长的搜索过程；(2) 扩张步，扩大步长的搜索过程；扩张步，扩大步长的搜索过程；(3)调整步，中心点、临界值的重新调整。调整步，中心点、临界值的重新调整。(1) 收缩步收缩步(2) 扩张步扩张步(3) 中心点和

15、步长的确定中心点和步长的确定 全局最优拟合的能力和效率很大程度上取决于初全局最优拟合的能力和效率很大程度上取决于初始点和步长，初始步长一般总不是很合适的，必须由始点和步长，初始步长一般总不是很合适的，必须由寻优过程的信息反馈调整。记录在寻优搜索过程中的寻优过程的信息反馈调整。记录在寻优搜索过程中的度点（即满足一定要求的参数点）的数量和位置，算度点（即满足一定要求的参数点）的数量和位置，算出它们平均数和标准差（出它们平均数和标准差（Sj为为bj的二阶原点矩）：的二阶原点矩）：12/ (/)/, CnjjjCCjjjCCjjTbbnnsSTnnlfs(4) 调整临界值调整临界值C 若 C 很小，产

16、生的度点数量太少，若 C 很大，产生的度点数量太多，这些情形都将使算法的能力和效率受损。临界值 C 必须有反馈调节机制。若N是每一轮次的试算节点总数，nE是扩张步一个循环（由37个轮次组成）的度点数量，在一次循环后重新计算临界值（包括步长）。前后两次循环(v, v+1)使用不同的公式是为了减少循环过程波浪形C值的发生。在mod(v,2)=0时，需将nE清零。oldEnewoldEexp( 1.3 (.65)/), mod( ,2)=1exp( 1.3 (1.3)/), mod( ,2)=0CnNNvCCnNNv （5）CE 算法的主要优缺点：算法的主要优缺点： 不必提供导数与偏导数，利于通用程

17、序的编制； 无需提供适合的初值； 实现最优拟合的能力较强； 搜索效率不高，对多参数非线性问题难于实施。 (6) 缩张算法的一些改进：缩张算法的一些改进： 1 每一轮次的试算节点数(z)随p的增加而指数（爆炸）式增长。5步点时，z=5p=exp(1.60944p)；在3步点时，z=3p=exp(1.09861p)。因此在p7（5步点）或p13（3步点）时，算法负荷量已超出普通pc机的上限（每轮次试算节点数以1m计），该法不适宜用于参数数p15的非线性方程的拟合。在参数较多(p8)时，只在 p 维参数空间中均匀随机布点，试算节点数在基础条件下随p的增加而呈多项式（二次式）增长（z=300+25*p

18、2），这比指数式增长大为减少，使多参数复杂非线性问题的拟合成为可能。 2 与解析法中的改良高斯牛顿法相结合，在给定的参数初值（或中间值）点处，利用参数微小差量的差分方程获得方程对某一参数的近似偏导函数值，再将各（观察值）点的偏导函数值的乘积累加，得到近似的Jaccobi矩阵(A，或A*)和常数阵K，再由A=K，解出=A-1K（=A*-1K），当接近 0 或RSS(Q)小于收敛标准时结束。 f 依第依第j个参数个参数bj的近似偏导数为：的近似偏导数为：0(0)(0)(0)(0)(0)(0)jjj()()()()() ()2iijijijijjjjfffffbbbbbbb (0,0,0)jj (0

19、)()ijfb 是是Xi及参数点及参数点 bj(0)处仅第处仅第bj 参参数具微小差值时的回归值。数具微小差值时的回归值。是是bj 微小差值参数增量；微小差值参数增量；bi(0)基于数值微分基础的改良高斯基于数值微分基础的改良高斯-牛顿法：牛顿法： 0ijfb001niijkijkffabbA = K*-1=A K(0)b=b + k当当b与与b(0)有差异时，应令有差异时，应令b替代替代b(0)重新计算重新计算，当当接近接近0或小于收敛标准时结束。或小于收敛标准时结束。 001()nijiiijfKYfb构建构建矩阵矩阵4 非线性回归统计数的假设测验 Jaccobi阵A的逆阵C (C=A-1

20、)对角线元素为相应回归统计数标准化的方差，所谓标准化的方差是指离回归误差方差为1时的方差。因此，第j个回归统计数bj（与0的差异显著性）测验可用 t 测验：jjbbts1jbjjQsCnp5 曲线、曲面拟合的matlab命令b,R,J=Nlinfit(x, y,fun,b0)b, resnorm=Lsqcurvefit(fun, b0, x, y)Polyfit(x, y, n)，Tool中basic fittingNlintool(x, y, fun, b0)gnlin(X, y, b0)6 一些实例231234235671bb xb xb xYb xb xb x1exp()K yabx(1exp()mK yabx2222225 1217 212() (1)(1) 23512134 1126(1)()bxxxb xxxY b bx eb bxxx ebe Thank you !