利用最小二乘法进行数据拟合
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
利用最小二乘法进行数据拟合:
例1.
在某个低温过程中,函数y 依赖于温度()C θ的试验数据如下表:
已知经验公式的形式为2y a b θθ=+,根据最小二乘法原理编制MATLAB 程序求出,a b ,并做相应的理论分析。
解:在两个观测量中,往往总有一个量精度比另一个高得多,为简单起见把精
度较高的观测量看作没有误差,并把这个观测量选作x ,而把所有的误差只认为是y 的误差。设x 和y 的函数关系由理论公式
y =f (x ;c 1,c 2,……c m ) (0-0-1)
给出,其中c 1,c 2,……c m 是m 个要通过实验确定的参数。对于每组观测数据(x i ,y i )i =1,2,……,N 。都对应于xy 平面上一个点。若不存在测量误差,则这些数据点都准确落在理论曲线上。只要选取m 组测量值代入式(0-0-1),便得到方程组
y i =f (x ;c 1,c 2,……c m ) (0-0-2) 式中i =1,2,……,m.求m 个方程的联立解即得m 个参数的数值。显然N 在N>m 的情况下,式(0-0-2)成为矛盾方程组,不能直接用解方程的方法求得m 个参数值,只能用曲线拟合的方法来处理。设测量中不存在着系统误差,或者说已经修正,则y 的观测值y i 围绕着期望值 ()()[]⎪⎭⎪ ⎬⎫⎪⎩⎪⎨⎧--=2 2 212,......,,;exp 21 i m i i i i c c c x f y y p σσπ, 式中i σ 是分布的标准误差。为简便起见,下面用C 代表(c 1,c 2,……c m )。考虑各次 测量是相互独立的,故观测值(y 1,y 2,……c N )的似然函数 ( ) ()[]⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧--= ∑=N i i i N N C x f y L 12 2 21;21ex p (21) σσσσπ. 取似然函数L 最大来估计参数C ,应使 ()[]min ;1 1 2 2 =-∑=N i i i i C x f y σ (0-0-3) 取最小值:对于y 的分布不限于正态分布来说,式(0-0-3)称为最小二乘法准则。若 为正态分布的情况,则最大似然法与最小二乘法是一致的。因权重因子2 /1i i σω=,故式 (0-0-3)表明,用最小二乘法来估计参数,要求各测量值y i 的偏差的加权平方和为最小。 根据式(0-0-3)的要求,应有 ()[] () m k C x f y c c c N i i i i k ,...,2,10 ;1 ˆ1 22==-∂∂==∑σ 从而得到方程组 ()[]()()m k C C x f C x f y c c N i k i i i ,...,2,10;;1 ˆ1 2==∂∂-==∑σ (0-0-4) 解方程组(0-0-4),即得m 个参数的估计值 m c c c ˆ,...,ˆ,ˆ21,从而得到拟合的曲线方程 ()m c c c x f ˆ,...,ˆ,ˆ;21。 然而,对拟合的结果还应给予合理的评价。若y i 服从正态分布,可引入拟合的x 2 量, ()[]∑ =-=N i i i i C x f y x 1 2 2 2 ;1 σ (0-0-5) 把参数估计 ()m c c c c ˆ,...,ˆ,ˆˆ21=代入上式并比较式(0-0-3),便得到最小的x 2 值 ()[]∑ =-=N i i i i c x f y x 1 22 2 min ˆ;1 σ (0-0-6) 可以证明,2 m in x 服从自由度v =N-m 的x 2分布,由此可对拟合结果作x 2检验。 由x 2 分布得知,随机变量2m in x 的期望值为N-m 。如果由式(0-0-6)计算出2 m in x 接近N-m (例 如m N x -≤2 min ),则认为拟合结果是可接受的;如果 22min >--m N x ,则认为拟合结果与观测值有显著的矛盾。 把2y a b θθ=+化为线性二乘拟合,即: y a b θθ =+ 于是就可以写编程序了。先 写一个主程序zxecf.m ,然后写一个被调用程序phi_k.m Matlab 程序如下: function S = zxecf(x,y,n,w) global i;global j; if nargin < 4 w = 1; end if nargin < 3 n = 1; end phi2 = zeros(n+1); for i = 0:n for j = 0:n phi2(i+1,j+1) = sum((w.*phi_k(x,i)).*phi_k(x,j)); end end phif = zeros(n+1,1); for i = 0:n