数值分析论文--曲线拟合的最小二乘法

数值分析论文--曲线拟合的最小二乘法

曲线拟合的最小二乘法 姓 名： 徐志超 学 号： 2019730059 专 业： 材料工程 学 院：

材料科学与工程学院 科 目：

数值分析 曲线拟合的最小二乘法 一、 目的和意义 在物理实验中经常要观测两个有函数关系的物理量。

根据两个量的许多组观测数据来确定它们的函数曲线， 这就是实验数据处理中的曲线拟合问题。 这类问题通常有两种情况：

一种是两个观测量 x 与 y 之间的函数形式已知， 但一些参数未知， 需要确定未知参数的最佳估计值； 另一种是 x 与 y 之间的函数形式还不知道， 需要找出它们之间的经验公式。

后一种情况常假设 x 与 y 之间的关系是一个待定的多项式， 多项式系数就是待定的未知参数， 从而可采用类似于前一种情况的处理方法。

在两个观测量中， 往往总有一个量精度比另一个高得多， 为简单起见把精度较高的观测量看作没有误差， 并把这个观测量选作 x， 而把所有的误差只认为是y 的误差。

设 x 和 y 的函数关系由理式 y＝f（x； c1， c2， cm）

1 / 13

（0-0-1） 给出， 其中 c1， c2， cm 是 m 个要通过实验确定的参数。

对于每组观测数据（xi， yi） i＝1， 2， ， N。 都对应于 xy 平面上一个点。

若不存在测量误差， 则这些数据点都准确落在理论曲线上。

只要选取 m 组测量值代入式 （0-0-1），便得到方程组 yi＝f （x；c1，c2，cm） （0-0-2） 式中 i＝1， 2， ， m.求 m 个方程的联立解即得 m 个参数的数值。 显然Nm 时， 参数不能确定。

在 Nm 的情况下， 式（0-0-2） 成为矛盾方程组， 不能直接用解方程的方法求得 m 个参数值， 只能用曲线拟合的方法来处理。 设测量中不存在着系统误差，或者说已经修正， 则 y 的观测值 yi 围绕着期望值 f （x； c1， c2， cm） 摆动， 其分布为正态分

布

，

则

yi

的

概

率

密

度

为

式中是分布的标准误差。

为简便起见， 下面用 C 代表（c1， c2， cm）。

考虑各次测量是相互的， 故观测值（y1， y2， cN） 的似然

函

数

取似然函数 L 最大来估计参数 C， 应使

（0-0-3） 取

---------------------------------------------------------------最新资料推荐------------------------------------------------------

最小值：

对于 y 的分布不限于正态分布来说， 式（0-0-3） 称为最小二乘法准则。

若为正态分布的情况， 则最大似然法与最小二乘法是一致的。 因权重因子

， 故式（0-0-3） 表明， 用最小二

乘法来估计参数， 要求各测量值 yi的偏差的加权平方和为最小。 根据式（0-0-3） 的要求， 应有

从而得到方程组

（0-0-4） 解方程组（0-0-4）， 即

得 m 个参数的估计值mccc , 1,...,2， 从而得到拟合的曲线方程

。

然而， 对拟合的结果还应给予合理的评价。

若 yi 服从正态分布， 可引入拟合的 x2 量，

（0-0-5） 把参

数估计最

小

的

代入上式并比较式（0-0-3）， 便得到

x2

值

（0-0-6） 可

以证明，2minx服从自由度 v＝N-m 的 x2 分布， 由此可对拟合结果作 x2 检验。

3 / 13

由 x2 分布得知， 随机变量2minx的期望值为 N-m。 如果由式 （0-0-6）计算出2minx接 近 N-m （ 例 如mN

） ，则 认 为 拟 合 结 果 是 可 接 受 的 ；如 果， 则认为拟合结果与观测值有显著的矛盾。

在我们研究两个变量(x, y)之间的相互关系时， 通常可以得到一系列成对的数据(x1, y1、 x2, y2... xm , ym)； 将这些数据描绘在 x -y 直角坐标系中(如图 1), 若发现这些点在一条直线附近， 可以令这条直线方程如(式 1-1)。

Y 计= a0 + a1 X (式 1-1) 其中：

a0、 a1 是任意实数 为建立这直线方程就要确定 a0 和 a1， 应用《最小二乘法原理》， 将实测值 Yi 与利用(式 1-1)计算值(Y 计=a0+a1X)的离差(Yi-Y 计)的平方和〔(Yi - Y计)2〕 最小为优化判据。

令: = (Yi - Y 计)2 (式 1-2) 把(式 1-1)代入(式 1-2)中得: = (Yi - a0 - a1 Xi)2 (式 1-3) 当(Yi-Y 计)平方最小时， 可用函数 对 a0、 a1 求偏导数， 令这两个偏导数等于零。 (式 1-4) (式 1-5) 亦即：

m a0 + (Xi ) a1 = Yi (式 1-6) (Xi ) a0 + (Xi2 ) a1 = (Xi, Yi) (式 1-7) 得到的两个关于 a0、 a1 为未知数的两个方程组， 解这两个方程组得出：

a0 = (Yi) / m - a1(Xi) / m (式 1-8) a1 = [Xi Yi - (Xi Yi)/ m] / [Xi2 - (Xi)2 / m)] (式 1-9) 这时把 a0、 a1 代

---------------------------------------------------------------最新资料推荐------------------------------------------------------

入(式 1-1)中, 此时的(式 1-1)就是我们回归的元线性方程即： 数学模型。

在回归过程中， 回归的关联式是不可能全部通过每个回归数据点(x1, y1、 x2, y2...xm,ym),为了判断关联式的好坏,可借助相关系数R， 统计量F， 剩余标准偏差S进行判断； R越趋近于 1 越好； F的绝对值越大越好； S越趋近于 0 越好。

R = [XiYi - m (Xi / m)(Yi / m)]/ SQR{[Xi2 - m (Xi / m)2][Yi2 - m (Yi / m)2]} (式 1-10) 在(式 1-1)中， m 为样本容量， 即实验次数； Xi、 Yi 分别任意一组实验 X、Y 的数值。 微积分应用课题一 最小二乘法 从前面的学习中, 我们知道最小二乘法可以用来处理一组数据, 可以从一组测定的数据中寻求变量之间的依赖关系, 这种函数关系称为经验公式. 本课题将介绍最小二乘法的精确定义及如何寻求 与 之间近似成线性关系时的经验公式. 假定实验测得变量之间的 个数据 , , , , 则在 平面上, 可以得到 个点 , 这种图形称为散点图, 从图中可以粗略看出这些点大致散落在某直线近旁, 我们认为 与 之间近似为一线性函数, 下面介绍求解步骤. 考虑函数 , 其中 和 是待定常数. 如果 在一直线上, 可以认为变量之间的关系为 . 但一般说来, 这些点不可能在同一直线上. 记 , 它反映了用直线 来描述 , 时, 计算值 与实际值 产生的偏差. 当然要求偏差越小越好, 但由于 可正可负, 因此不能认为总偏差 时, 函数 就很好地反映了 变量之间的关系,

5 / 13

因为此时每个偏差的绝对值可能很大. 为了 改进这一缺陷, 就考虑用 来代替 . 但是由于绝对值不易作解析运算, 因此, 进一步用 来度量总偏差. 因偏差的平方和最小可以保证每个偏差都不会很大. 于是问题归结为确定 中的常数 和 , 使 为最小. 用这种方法确定系数 , 的方法称为最小二乘法. 在 科 学 实 验 的 统 计 方 法 研 究 中 ，往 往 要 从 一 组 实 验 数 据

中， 寻找自变量 x 与因变量 y 之间的函数关系 。

由于观测数据往往不准确， 因此不要求

经过所有点

， 而只要求在 给 定 点ix 上 误 差 而 只 要 求 所

在 所 有 给 定 点ix 上 的 误 差

按某种标准最小。

若记

最小。

如果用最大范数， 计算上困难较大， 通常采用欧式范数为误差度量的标准。

的函数类型往往与实验的物理背景以及数据的实际分布

作

， 就是要求向量

的范数

有关， 它一般含有某些待定参数。 如果

是所有待定参数的线性函数， 那么相应的问题称

为线性最小二乘问题， 否则称为非线性最小二乘问题。 最小二乘法还是实验数据参数估计的重要工具。

这是因为这种方法比其他方法更容易理解，即使在其他方法失效

---------------------------------------------------------------最新资料推荐------------------------------------------------------

的情况下， 用最小二乘法还能提供解答， 而且从统计学的观点分析， 用该方法求得各项估计具有最优统计特征， 因此这一方法也是系统识别的重要基础。

线性最小二乘问题可以借助多元微分学知识通过求解法方程组得到解答。

用最小二乘法求拟合曲线时， 首先要确定

的形式。

这不单纯是数学问题， 还与所研究问题的运动规律以及所得观测数据图， 确定

有关； 通常要从问题的运动规律以及给定数据描的形式， 并通过实际计算选出较好的结果。

都

为了使问题的提法更有一般性， 通常把最小二乘法中的考虑为加权平方和

这里

上的加权函数， 它表示不同点

是

处的数据比重不同。

二、 计算实例 从随机的数据中找出其规律性， 给出其近似表达式的问题， 在生产实践和科学实验中大量存在， 通常利用数据的最小二乘法求得拟合曲线。

在某冶炼过程中， 根据统计数据的含碳量与时间关系， 试求含碳量 y 与时间 t的( t 分） 0 5 10 15 20 25

2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4. 拟合曲线。

7 / 13

1、 掌握曲线拟合的最小二乘法； 2、 最小二乘法亦可用于解超定线代数方程组； 3、 探索拟合函数的选择与拟合精度间的关系 本 题 要 求 我 们 用线 进 行 拟 合 ，这 里

对 曲

故

由于

可以利用此式算出

拟

合

曲

线

的

ia

，

即

所 以 求 得

，

，

，，

误差为

， 而均方误差

---------------------------------------------------------------最新资料推荐------------------------------------------------------

为下图可见实际测出值与拟合值的

差别， 下表可见拟合出的每一点的误差以及均方误差。

00. 511. 522. 533. 5(0. 0001)44. 550102030405060时间(s)含碳量真实值拟合值 t y 拟合值 误差 误差平方 0 0 0 0 0 5 1.27 1.20215 -0.067850.004603623 10 2.16 2.1662 0.0062 3.844E-05 15 2.86 2.91855 0.05855 0.003428103 20 3.44 3.4856 0.0456 0.00207936 25 3.87 3.375 0.02375 0.0005063 30 4.15 4.1694 0.0194 0.00037636 35 4.37 4.335 -0.031050.0009102 40 4.51 4.4288 -0.0812 0.00659344 45 4.58 4.46535 -0.114650.013144622 50 4.02 4.475 0.455 0.207025 55 4. 4.48415 -0.155850.0242222 均方误差 0.512932 三、 结构程序设计 在本题使用 Visual Studio c# .NET 编译程序。 //在窗体的 Load 事件里调用 InitialDeal 和 Deal 函数来处理数据。 private

void

ResultReport_Load(object

sender,

System.EventArgs e) { this.InitialDeal(); this.Deal(); } //初始化各个变量。

private void InitialDeal() { int i = 0; while(tString.Length

0)

{

t[i++]

=

Convert.ToDouble(tString.Substring(0,tString.IndexOf(,))) * tUnit; tString = tString.Remove(0,tString.IndexOf(,) +

9 / 13

1); } num = i; i = 0; while(yString.Length 0) { y[i++] = Convert.ToDouble(yString.Substring(0,yString.IndexOf(,))) * yUnit; yString = yString.Remove(0,yString.IndexOf(,) + 1); } i = 0; while(i num) { a11 += y[i]*y[i]; i++;

}

this.labelA11.Text

=

a11.ToString(#.00E0;(#.00E0);0.00); i = 0; while(i num) { a12 += y[i]*y[i]*y[i]; i++; } a21 = a12; this.labelA12.Text = a12.ToString(#.00E0;(#.00E0);0.00); this.labelA21.Text = a21.ToString(#.00E0;(#.00E0);0.00); i = 0; while(i num) { a13 += y[i]*y[i]*y[i]*y[i]; i++; } a31 = a13; a22 = a13; this.labelA31.Text = a31.ToString(#.00E0;(#.00E0);0.00); this.labelA13.Text = a13.ToString(#.00E0;(#.00E0);0.00); this.labelA22.Text = a22.ToString(#.00E0;(#.00E0);0.00); i = 0; while(i num) { a23 += y[i]*y[i]*y[i]*y[i]*y[i]; i++; } a32 = a23; this.labelA23.Text = a23.ToString(#.00E0;(#.00E0);0.00); this.labelA32.Text = a32.ToString(#.00E0;(#.00E0);0.00); i =

0;

while(i

num)

{

a33

+= }

y[i]*y[i]*y[i]*y[i]*y[i]*y[i]; i++;

this.labelA33.Text = a33.ToString(#.00E0;(#.00E0);0.00); i = 0; while(i num) { b1 += y[i]*t[i]; i++; } this.labelB1.Text = b1.ToString(#.00E0;(#.00E0);0.00); i = 0;

---------------------------------------------------------------最新资料推荐------------------------------------------------------

while(i num) { b2 += y[i]*y[i]*t[i]; i++; } this.labelB2.Text = b2.ToString(#.00E0;(#.00E0);0.00); i = 0; while(i num) { b3 += y[i]*y[i]*y[i]*t[i]; i++; } this.labelB3.Text = b3.ToString(#.00E0;(#.00E0);0.00); } //进行最小二乘法处理。 private

void

Deal()

{

double

ab11,ab12,ab13,ab21,ab22,ab23,ab31,ab32,ab33 = 0; double aDiterminal = 0; aDiterminal = a11*a22*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31 - a23*a32*a11 - a12*a21*a33; ab11 = (a22*a33 - a32*a23)/aDiterminal; ab12 = (a13*a32 - a12*a33)/aDiterminal; ab13 = (a12*a23 - a13*a22)/aDiterminal; ab21 = (a23*a31 - a21*a33)/aDiterminal; ab22 = (a11*a33 - a13*a31)/aDiterminal; ab23 = (a13*a21 - a11*a23)/aDiterminal; ab31 = (a21*a32 - a31*a22)/aDiterminal; ab32 = (a12*a31 - a11*a32)/aDiterminal; ab33 = (a11*a22 - a12*a21)/aDiterminal; // 声明一个 3 行 3 列的两维数组， 用于存放矩阵 A // 初始化矩阵 A double[,] matrixA = new double[3,3]{ {ab11,ab12,ab13}, {ab21,ab22,ab23}, {ab31,ab32,ab33} }; // 声明一个 3 行 1 列的两维数组， 用于存放矩阵 B // 初始化矩阵 B double[,] matrixB = new double[3,1]{{b1},{b2},{b3}}; // 声明一个 3 行 1 列的两维数组， 用于存放矩阵 A 和矩阵 B 的乘积 double[,]

11 / 13

matrixC = new double[3,1]; // 计算矩阵 A 的逆转 // 计算矩阵 A 和矩阵 B 的乘积 for(int i = 0; i 3; i++) { for(int j = 0; j 1; j++) { // 初始化矩阵 C matrixC[i,j] = 0; // 计算矩阵 A 和矩阵 B 的乘积， 并把值存放在矩阵 C 中 for(int k = 0; k 3; k++) { matrixC[i,j] += matrixA[i,k] * matrixB[k,j]; } } } this.labela1Result.Text

matrixC[0,0].ToString(#.00E+0;#.00E-0;0.00); this.labela2Result.Text

matrixC[1,0].ToString(#.00E+0;#.00E-0;0.00); this.labela3Result.Text

= = =

matrixC[2,0].ToString(#.00E+0;#.00E-0;0.00); } 四、 结果讨论与分析 最小二乘法如果想将曲线拟合的比较完美， 必须应用适当的模拟曲线， 如果模拟曲线选择不够适当， 那么用最小二乘法计算完后， 会发现拟合曲线误差比较大， 均方误差也比较大， 而如果拟合曲线选择适当， 那么效果较好。

因此， 需要对已知点根据分布规律选取多个可能的近似拟合曲线， 算出后比较误差与均方误差， 得到最佳拟合曲线。

但是如果已知点分布非常不规律， 无法观察或是无法正确观察出其近似曲线， 那么根本无法使用最小二乘法进行曲线拟合， 我们只能使用其它方法进行逼近， 如最佳一致逼近多项式。 如果本题我们错误的使用

对本题进行模拟，

---------------------------------------------------------------最新资料推荐------------------------------------------------------

那样误差与均方误差都 非 常 大 ，肯 定 无 法 得 到 好 的 效 果 ，相 比 较 而 言 ，本 题 所 选 取 的

可以达到预期效果。

13 / 13