西北农林科技大学实验报告学院名称：理学院 专业年级：2013级信计1班 姓 名： 学 号 课 程：数学模型与数学建模 报告日期：2013年12月1日拟合模型与回归分析实验目的配合《数学建模与数学模型》的第3章“常见的模型及其组建”，介绍如何运用数学软件进行模型组建，并结合数学理论分析求解模型。

拟合模型的组建是通过对有关变量的观测数据（散点图）的观察、分析。

结合问题背景，运用数学分析，选择当前恰当的数学表达方式得到的。

拟合的目的是寻找一条光滑曲线y=ψ(x),能够很好地表现受随机因素干扰的观测数据(){}ni i i y x 1,=所反映的规律。

原则上尽量选择简单的数学公式表达规律，在简单的数学表达式中选择拟合效果好的。

一、赛跑成绩与赛跑距离1 实验题目赛跑成绩与赛跑距离2 实验问题陈述下面的表2.1.1给出了1997年以前6个不同距离的中短距离赛跑成绩的世界纪录：3 实验内容解 共分4个步骤，分别叙述如下。

步骤1 在坐标系上画出观测数据的散点图。

>> X=[100 200 400 800 1000 1500];>> Y=[9.95 19.72 43.86 102.4 133.9 212.1]; >> plot(X,Y,'*')步骤2 根据散点图，取线性拟合模型y=a+bx.步骤3 利用数据（x i ,y i ）估计模型参数a,b 。

就是在寻找超定方程（方程个数多于未知数的个数）Ad =y ′的近似解d =(a,b)′,其中⎪⎪⎪⎪⎭⎫⎝⎛=n x x A ...1...11，⎪⎪⎪⎪⎭⎫ ⎝⎛=n y y ...y ′1 称X=（x 1,x 2,....,x n ）′为设计矩阵。

采用最小二乘法确定参数的估计值∧a ,∧b ,也就是求拟合残差平方和∑=--=ni i i bx a y Q 12)(的最小值（a,b)。

下面利用MATLAB 指令完成参数估计。

>> A=[ones(size(X))',X']; >> d=A\Y';>> z=d(1)+d(2).*X; ;得到线性模型：y=-9.99+0.145x. 步骤4 分析拟合效果，做拟合图。

>> plot(X,Y,'*',X,z,'LineWidth',2) >> Q=sum((Y-z).^2)简单的根据拟合残差图和拟合的残差平方和Q=81.76看，拟合的效果不是特别糟糕，但是，结果不符合实际。

根据拟合得到的模型，当x<68.89m时，跑步时间y<0，显然不正确。

实际上当跑步距离为0时，所需要的时间也为0.在前面选择模型时没有考虑到实际问题的这一基本要求，因此导致矛盾的结果。

修正模型，要求拟合函数满足条件y（0）=0，并根据散点图特点，取幂函数模型：y=ax b.为了利用线性拟合指令，令z=lny,u=lnx,a*=lna,则幂函数拟合问题转变为线性拟合z=a*+bu.>> A=[ones(size(X))',log(X)'];D=A\log(Y)';d0=[exp(D(1)),D(2)];fun=inline('d(1).*X.^(d(2))','d','X');>> Q1=sum((Y-fun(d0,X)).^2);于是得到幂函数模型y=0.048x1.145,结果比较符合实际。

但是这样拟合得到的不是使得残差平方和∑=--=nii ibxayQ12)(达到最小的参数（a,b）,为了改进拟合效果，可以进一步利用MATLAB的非线性拟合指令。

由于非线性拟合求最小值点通常使用迭代逼近算法，需要先输入参数估计值作为初始值（a0,b）.因此选择前面通过线性化方法得到的参数拟合值作为下一步非线性拟合的参数初始估计值。

>> d=nlinfit(X,Y,fun,d0);>> Q2=sum((Y-fun(d,X)).^2);这样得到幂函数模型：y=0.0416x1.1678,残差平方和为Q2=6.1319,可见非线性拟合极大的改进了拟合效果。

注意，拟合模型通常也称为经验模型。

二、投资预测1 实验题目投资预测2 实验问题陈述研究某地区实际投资额与国民生产总值（GNP）及物价指数（ICP）的关系，以便根据对未来国民生产总值及物价指数的估计，预测未来的实际投资额。

附：以往20年数据表如表2.1.4所示。

3 实验内容解现在我们按回归分析方法讨论问题。

首先，表述问题，选择变量。

为确定实际投资额对国民生产总值和物价指数的依赖关系，取实际投资额为因变量y，国民生产总值和物价指数分别为自变量x1和x2.然后，进行数据描述分析。

有散点图可见y线性依赖x1和x2，而且变化趋势很相似，怀疑x1与x2之间存在共线性性质，画x1-x2散点图马上证实了这一点。

指令如下：>> x1=[596.7 637.7 691.1 756 799 873.4 944 992.7 1077.6 118 5.9 1326.4 1434.2 1549.3 1718 1918.3 2163.9 2417.8 2631.7 29 54.7 3073];>> x2=[0.7167 0.7277 0.7436 0.7676 0.7906 0.8254 0.8679 0.914 5 0.96011 1.0575 1.0575 1.1508 1.2579 1.3234 1.4005 1.5042 1 .6342 1.7842 1.9514 2.0688];>> y=[90.9 97.4 113.5 125.7 122.8 133.3 149.3 144.2 166.4 1 95 229.8 228.7 206.1 257.9 324.1 386.6 423 401.9 474.9 424.4];>> subplot(1,3,1),plot(x1,y,'*'),title('x1-y')>> subplot(1,3,2),plot(x2,y,'*'),title('x2-y')>> subplot(1,3,3),plot(x1,x2,'*'),title('x1-2')因此，实际投资额y可以表示成其中一个自变量的函数，选择国民生产总值，取线性模型y=a+bx做回归分析，指令如下：x=x1>> A=[ones(size(x1))',x1'];[d,bint,r,rint,stats]=regress(y',A);plot(r,'*'),axis([0,20,-60,60]),title('residual')结果如表2.1.5所示：表2.1.5 投资额与国民生产总值的回归结果参数估计值置信区间虽然，拟合优度R 2接近1，F 统计量概率值P ＜0.0001很小，但是参数估计的95%置信区间太大，而且含有零点，这意味着参数参数有可能取零值。

特别是残差序列图2.1.5出现异方差现象，残差散布的范围随着序列变化增大。

这与回归分析成立的前提“残差具有零均值和均方差”相矛盾。

考虑到投资额和国民生产总值这些数据（x t ,y t )都是来自同一个体的不同时间t 的观测值，不同时间的数据之间可能存在相关性，这种相关性简称为自相关性。

自相关性分析也称为自回归分析，是研究时间序列的常用方法。

但不是对所有时间序列数据都可以直接进行自回归分析，希望利用过去的数据预测未来的关系，就必须假设两个变量之间未来的依赖关系与过去的有着某种相似性，统计上定义这种相似性为时间序列的平稳性。

严格的说，称一个时间序列{r t }是平稳的，如果该序列满足：对任意的整数k ，任意的的时间点t 0，随机变量r t0,r t0+1,...,r t0+k 是独立同分布的。

也就是说该序列的均值和方差不会随时间的改变而变化。

从上面的残差图可见，对于k=0，残差序列{r t =y t -a-bx t }的方差随时间逐渐增大，它不是一个平稳过程，自相关性也非常不好，因此不能采用自回归模型。

重新考虑到作为时间序列，实际投资额对国民生产总值的依赖可能存在滞后，国民生产总值对实际投资额的部分影响可能隔几年后才显现出来。

经过多次试验，得到统计分析结果最佳的模型：y t =a+b 1x t-2+b2x t对这个模型进行模型回归分析，指令如下：>> A=[ones(size(x1(3:end)))',x1(1:end-2)',x1(3:end)']; >> [d,bt,r,rt,sts]=regress(y(3:end)',A); >> plot(r,'*'),grid结果表明，当年的国民生产总值与实际投资额是正相关的，前年的国民生产总值对实际投资额的影响是抑制的。

根据这个模型，只要知道国民生产总值就不难估计相应时间的实际投资额。

步骤5 预测预报建立拟合模型的一个主要目的是为了进行预测预报，在x 0处的拟合模型预测值的点估计为0∧∧∧0b a x y +=，多数情况下这个预报值是不可靠的。

科学的方法是给出区间估计，即给出95%置信度的预测区间|y -∧0y |≤L.为了确定区间半径L,需要知道随机变量y -∧0y 的概率分布。

以一元线性拟合模型为例说明设定置信区间的原理和方法。

引理2.1.1 设()()2σ,0~N r bx a y i i i =+-服从正态分布，则由最小二乘法得到的拟合模型x b a ∧∧∧y +=满足（1）⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛+-xx l x n a N a //1σ~22∧； （2）()xx l b N b /σ,~2∧； （3）()2-n χ~σ/2Q ；于是，在x 0处的拟合模型预测值是⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛+++=-xx l x n bx x N x b a y //1σ,~2200∧∧∧0 因为，观测值()20000σ,~bx a N r bx a y +++=所以⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛++--xx l x n N y y //11σ,0~22∧00 取统计量1)]-(n σQ/[√/x 1/n 1√σ/2-2∧00⎪⎪⎭⎫ ⎝⎛++⎪⎭⎫ ⎝⎛-=xx l y y T可以证明，)2(~-n t T ，即服从自由度为n-2的t 分布。