2023年最小二乘拟合实验报告

南昌工隹学院《必算方法》实验赧告课程系专班学生课程系专班学生业信息与计算科学级12级一班姓名魏志辉学号《最小二乘求解》1引言在科学实验和生产实践中，经常要从一组实验数据(".,%)(，=1,2,…，加出发,寻求函数y=f(x)的一个近似表达式y=4)(x),称为经验公式，从几何上来看，这就是一个曲线拟合的问题。多项式的插值虽然在一定限度上解决了由函数表求函数近似表达式的问题，但用它来解决这里的问题，是有明显的缺陷的。一方面，由实验提供的数据往往有测试误差。假如规定近似曲线y=6（x）严格地通过所给的每个数据点（乙,匕），就会使曲线保存本来的测试误差,因此当个别数据的误差较大的时候，插值的效果是不抱负的。另一方面，当实验数据较多时,用插值法得到的近似表达式，明显缺少实用价值。在实验中，我们经常用最小二乘法来解决这类问题。定义&=g）-匕为拟合函数在七处的残差。为了是近似曲线能尽量反映所给数据点的变化趋势，我们规定I心I尽也许小。在最小二乘法中，我们选取以外,使得偏差平方和最小,州於即£耳=Z〔以巴）-=m行，这就是最小二乘法的原理。/=i/=i2实验目的和规定运用matlab编写.m文献，规定用最小二乘法拟定参数。以下一组数据中x与y之间存在着y=四加的关系，运用最小二乘法拟定式中的参数a和b,并计算相应的军方误差与最大偏差。数据如下:Xi234567891()y0.8982.383.071.842.021.942.222.774.024.76X111213141516171819y5.466.5310.916.522.535.750.661.681.83算法原理与流程图（1）原理最小二乘是规定对于给定数据列（乙,匕）（，=L2,…，勿），规定存在某个函数类①={.（X）,.（X）,…纥（X）}（〃＜/）中寻求一个函数："（＞）=W.（x）+瓶⑺+…+或（p\x）,使得夕*（＞）满足

Z［?*（巧）-之］2=叫心Z［以七）一之F。/=1①/=1根据以上条件可知，点（吊,吊,…，耳）是多元函数S&,4,…，%）=20%%（巧）-兀］2/=!Ar=0的极小点，从而办父，…，a：满足方程组，一=0（A=0,1,•••,n）da,AYa即qZ然（X1）%（巧）+qZ（Pk（X）（p\x）+…+”（Z）/（乙）=Z外（储）匕,AY/=12=12=1/=1记（力,#=£力（匕）双乙）,则上述方程组可表达成7=1综（然，°。）+a、®，，）+…+练（/，/）=（%，，），（k=0/，•・・，n）（%，%）防%）…（%，%）防%）…M，（P）

（件必）（件/）…（件外）..._（%,%）3,外）…（外，/）（%，%）防%）…M，（P）

（件必）（件/）…（件外）（%，%）防%）…M，（P）

（件必）（件/）…（件外）..._（%,%）3,外）…（外，/）M，f）M'f）,这个方程组成为法方程组，可以证明，当*

*.（%，f）_%（x）,%（X）,…0〃（x）线性无关时，它有唯一解。特别地，曲线拟合的一种常用情况为代数多项式，即取/（X）=1,%（才）=X,…?⑺=X"，则（化.，/）=产2=1/=1®,「）=Zx5i（k=O,l,--\n）7-1故相应的法方程组变为

Ex；y,Ex；y,这就是最小二乘法的原理。在解决本题时，为了简便起见，我们将指数转变成代数多项式去计算。在夕=ae及两边取对数，得到Iny=Ina+^,取y⑴=In匕x⑴=x,可见y⑴，x⑴是呈线性关系的。这样我们可以方便地运用最小二乘法求取参数。（2）流程图.\.陆入1'2,…‘⑼油'一i=l»2,…，m生乎”］层组的系数矩解注,租组俎为G'生乎”］层组的系数矩解注,租组俎为G'=0,1,••,n)整体流程图生成矩阵C流程图4程序代码及注释整体流程图括最小二乘拟合%&为线性拟合中的常数，b为一次项系数%t为均方误差，maxi为最大偏差function[a,b,maxi]=po1yfit(xO,y0,n)m=length(xO);p=length(yO);%xO和yO长度不等时,报错ifm~=pfprintf(*Error!Pleaseinputagain!\n');end%生成中间矩阵Cfori=l：mC(i,1)=1;forj=2:n+lC(i,j)=x0(i)*C(i,j-1);endend%生成系数矩阵AA=C**C；%将题目中的y的每项求A然对数后得到ylfori=l：my1(i)=log(yO(i));end生生成法方程组的右端向量BB=C'*yl1;告求解拟合系数X=A\B；%题中，a=exp(X(1)),b=X(2)a=exp(X(1));b=X(2);%先求偏差平方和，再求均方误差sum=0;fork=1：my2(k)=a*exp(b*x0(k));1(k)=y0(k)-y2(k)；sum=sum+1(k).A2；endt=sqrt(sum);方最大偏差为偏差矩阵中绝对值最大的一项maxi=max(max(abs(1)));end5算例分析1、测试示例»x=1:18;»y=[0.8982.383.071.842.021.942.222.774.024.765.466.5310.916.522.535.750.661.681.8];>>[abtmaxJ=polyfit(x,y,1)Error!Pleaseinputagain!a=0.7185b=0.2227t=31.6255max=22.06312、计算过程(I)一方面输入已知点>>x=1:19;»y=L0.8982.383.071.842.021.942.222.774.024.765.466.5310.916.522.535.750.661.681.81;(2)输出结果>>[abtmaxl=polyfit(x,y,1)a=0.6814b=0.2306t=38.3255max=27.3047其中,a,b,t,max分别为常数项，一次项系数，均方误差，最大偏差。6讨论与结论1、时间复杂度：>>tic;[abtmax]=polyfit(x,y,l);tocElapsedtimeis0.859861seconds.说明该算法具有一定的复杂性。2、直观展示：输入以下命令：»x=l:l9;»y=[0.8982.383.071.842.021.942.222.774.024.765.466.5310.916.522.535.750.661.681.8J；»fori=l:19y0(i)=1og(y(i));end»x0=0:0.01:20;»yl=0.6814*exp(0.2306*x0);»plot(x,y,'+')>>holdon»plot(x0