统计软件应用课程设计

1、Monte Carlo 在线性回归中的简单应用班级：1324001姓名：徐萍萍学号：201320400118指导老师：闫霏霏电话号码：2016年4月17日摘要：本文针对确定运动员的耗氧量与其他一些因素的关系数据，在SAS中进行回归分析得到对数据拟合较好的线性模型。再用Monte Carlo随机过程产生残差项并代入线性方程中，分别假定残差项不符合均值为零、正态分布、异方差三大假设，从而检验回归理论。关键字：蒙特卡罗、多元线性回归、残差项基本假设问题重述：在运动生理学的研究中，为了确定运动员的耗氧量与其他一些因素的关系，在一个实验中对31个人测量了年龄（age），体重（weight），跑完1.5英

2、里用的时间（runtime），静态时的心率（rstpulse）,跑动时的心率（runpulse）,跑步时的最大心率（maxpulse）,每公斤体重每分钟的耗氧量（oxy）。实测数据（oxy.txt）见下表，试以oxy为因变量，估计该变量对于问题中所有其他变量的直线回归方程。并用Monte Carlo验证回归理论中残差项三大基本假设。(数据在附录中)第1章 、基本理论一、蒙特卡罗方法蒙特卡罗（Monte Carlo）方法概述 蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解

3、的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。蒙特卡罗方法的基本原理 由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡罗法正是基于此思路进行分析的。设有统计独立的随机变量Xi(i=1，2，3，k)，其对应的概率密度函数分别为fx1，fx2，fxk，功能函数式为Z=g(x1，x2

4、，xk)。 首先根据各随机变量的相应分布，产生N组随机数x1，x2，xk值，计算功能函数值 Zi=g(x1，x2，xk)(i=1，2，N)，若其中有L组随机数对应的功能函数值Zi0，则当N时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。蒙特卡罗方法的收敛性设所求的量是随机变量的数学期望E(x)，那么Monte Carlo方法通常使用随机变量的简单子样的算术平均值，即 作为所求量X的近似值。由柯尔莫哥罗夫（Kolmogorov）大数定理可知，即当N充分大时，有成立的概率等于1，亦即可以用作为所求量X的估计值。根据中心极限定理，如果随机变量的标准差不为零，那么Monte Ca

5、rlo方法的误差为式中，为正态差，是与置信水平有关的常量。Monte Carlo方法的收敛速度的阶为，误差是由随机变量的标准差S和抽样次数N决定的。二、多元线性回归模型多元线性回归模型的一般形式 设随机变量y与一般变量的线性回归模型为： 记为式中，是P+1个未知参数，称为回归常数，称为回归系数。y称为被解释变量（因变量），是P个可以精确测量并控制的一般变量，称为解释变量（自变量）。P=1时，式即一元线性回归模型；P2时，我们就称式为多元线性回归模型。是随机误差，与一元线性回归一样，对随机误差我们常规定：称 为理论回归方程。 对一个实际问题，如果我们获得n组观测数据，则线性回归模型式可表示为：写

6、成矩阵形式为： 记位式，式中 X是一个阶矩阵，称为回归设计矩阵或资料矩阵。在实验设计中X的元素是预先设定并可以控制的，人的主观因素可作用其中，因而称X为设计矩阵。多元线性回归模型的基本假定 回归方程式有如下基本假定：解释变量x1,x2,xp是确定性变量，不是随机变量，且要求rank(X)=p+1<n.这里的rank(X)=p+1<n，表明设计矩阵X中的自变量列之间不相关，样本量的个数应大于解释变量的个数，X是一满秩矩阵。随机误差项具有零均值和等方差，即 这个假定常称为高斯-马尔柯夫条件。，即假设观测值没有系统误差，随机误差项的平均值为零。随机误差项的协方差为零，表明随机误差项在不同

7、的样本点之间是不相关的（在正态假定下即为独立的），不存在序列相关，并且有相同的精度。正态分布的假设条件为：对于多元线性回归的矩阵模型式，这个条件可表示为：第二章、模型建立一、SAS中线性回归过程 这里在SAS软件中，对问题中数据（oxy.txt）做多元线性回归，运行程序如下：I：libname ep 'e:saslx'run; data ep.oxys; infile 'e:sasdtoxy.txt'input age weight oxy runtime rstpulse runpulse maxpulse ;proc reg;model oxy= age w

8、eight runtime rstpulse runpulse maxpulse;run;II：proc stepwise data=ep.oxys;*逐步回归法;model oxy=runtime age weight runpulse maxpulse rstpulse;run;quit; 首先用程序I做回归分析，得输出结果。从输出结果(图一）中可知： 方差分析中F=22.43，P<0.0001回归方程显著。但应注意所检验的统计假设是六个自变量的系数都为零。若拒绝该假设则表示至少有一个系数不为零，即至少有一个自变量对因变量oxy的影响是有意义的。所以，此时不能说所有自变量对oxy的影

9、响都有统计意义。 决定系数R-square=0.8487，表示所得回归防城对数据的线性拟合很好。它的平方根为因变量oxy与6个自变量的复相关系数。 Parameter Estimates部分是检验常数项或自变量系数等于零的假设。由检验结果可知，当显著水平为0.05，体重weight(p=0.1869)和休息时的心率rstpulse(p=0.7473)这两个变量应删除。一般来说，某个自变量是否对因变量有统计意义，不仅决定于这个变量本身，还跟回归方程中同时存在的其他自变量有关。若把对因变量无统计意义的自变量留在回归方程，会增大参数估计或误差，因此这些变量应剔除。 由以上分析可知，应采用逐步回归法S

10、TEPWISE过程，即采用II中程序。在逐步回归过程的结果表中，可以看到引进变量的顺序、每一步引进的变量对R-square的贡献和总的R-square的变化，C（P）值的变化及引进每个变量时的检验概率。其中C（P）值最小的最后那个模型为最适合的回归模型。最后得估计每公斤体重每分钟耗氧量（oxy）的线性方程为：y=0.27051x1-0.34811x2-2.76758x3-0.19773x4+98.14789其中x1、x2、x3、x4依次代表maxpulse、runpulse、runtime、age,体重weight和休息时的心率rstpulse被剔除。图一二、蒙特卡罗方法回归检验 由Monte

11、 Carlo思想，由正态随机数产生误差项e代入线性方程中，方程中自变量系数为在回归理论中估计得出的系数。这里通过分别假定残差项不符合均值为零、正态分布、异方差三大假设，从而检验回归理论。基于以上思想，在SAS中编写程序III如下：data ep.oxys;infile 'e:sasdtoxy.txt'input x4 weight y x3 rstpulse x2 x1 ;drop weight rstpulse y;e=0.001*rannor(10);y1=0.27051*x1-0.34811*x2-2.76758*x3-0.19773*x4+e;run;proc reg

12、data=ep.oxys;model y1=x1 x2 x3 x4 /noint ;run;quit; 读入数据，由前面分析体重weight和休息时的心率rstpulse两个变量应剔除，故读入数据后将其丢掉。由标准正态随机数rannor(random normal)产生误差项随机样本，前面乘以0.001是为避免数量级影响产生的误差项样本，从而影响估计的线性方程的精度。回归模型中去掉常数项，因为当显著水平为0.05时，包含常数项是截距项显然不能通过检验（见图二）。由结果中（图三）可以得到估计的线性方程为: y1=0.27053x1-0.34814x2-2.76749x3-0.19774x4+e选

13、择其中x1的系数为例，即与进行比较，计算估计系数与给定的真实系数的相对误差，则，可知误差较小，则利用Monte Carlo估计线性模型效果较好。图四中为产生的部分数据及残差项e的部分值。图二图三图四将III中程序划线两行程序修改如下：e=1+0.001*rannor(10);y1=0.27051*x1-0.34811*x2-2.76758*x3-0.19773*x4+e; 这里是对残差项均值为零的假设的检验,令e=1+0.001*rannor(10)，即使得残差均值为1。运行程序可得结果（图五）得估计的线性方程为: y1=0.27698*x1-0.35038*x2-2.76895*x3-0.1

14、9200*x4+e 上式中以自变量x1的系数为例，计算相对误差为2.391%，误差明显增大。图五将III中程序划线两行程序修改如下：e=1*uniform(10);y1=0.27051*x1-0.34811*x2-2.76758*x3-0.19773*x4+e; 这里是对残差项服从标准正态分布的检验，令e=1*uniform(10)产生符合均匀分布的残差随机样本，即不符合基本假设。运行程序可得结果（图六）中估计的线性方程为:y1=0.28343*x1-0.35837*x2-2.82705*x3-0.18253*x4+e上式中同样以x1的系数为例，计算相对误差为4.776% 误差明显增大。即用符

15、合均匀分布的残差随机样本，得到的线性方程估计系数偏差较大。图六将III中程序划线两行程序修改如下：if _n_<=15 then e=0.001*rannor(10);else e=1*uniform(10);y1=0.27051*x1-0.34811*x2-2.76758*x3-0.19773*x4+e; 这里是对同方差假设的检验，利用条件语句对前后样本赋予不同的残差来达到异方差的效果。运行程序可得结果（图七）中估计的线性方程为:y1=0.25787*x1-0.34389*x2-2.80289*x3-0.15246*x4+e计算出x1(maxpulse)系数估计值与真实值的相对误差为4

16、.671%，误差较与符合所有假设时的0.0073%有较大变化。图七第3章 、结论与体会 综合以上分析，先利用数据在SAS中得到拟合较为精确的线性回归方程，再利用Monte Carlo作为一个工具用于产生随机数。改变模型，分别产生不符合线性回归三个基本假设的残差项随机样本运用于新的模型中，得到的估计系数相较与符合基本假设的估计系数相对误差有较大变化。由此说明，要得到精确的线性回归模型必须符合基本假设，即达到了检验回归理论的目的。Monte Carlo（随机模拟）还可应用于很多有关于随机现象的近似计算方面，它还可以应用于概率计算、求积分的数值近似解、模拟置信区间理论、模拟假设检验。 在这次课程设计

17、过程中，首先很感谢老师的提点和同学们的帮助。自己也感觉受益良多，大到对Monte Carlo方法与回归理论的理解，小至论文的格式都进步不少，相信对将来毕业论文设计会有一定帮助。参考文献【1】何晓群、刘文卿 应用回归分析M 中国人民大学出版社，2015.【2】董大钧 统计分析应用M 电子工业出版社，2014.【3】林晓辰 蒙特卡罗方法 MBA智库百科，2008.附录age weight oxy runtime rstpulse runpulse maxpulse44 89.47 44.609 11.37 62 178 182 40 75.07 45.313 10.07 62 185 185 44

18、 85.84 54.297 8.65 45 156 168 42 68.15 59.571 8.17 40 166 172 38 89.02 49.874 9.22 55 178 180 47 77.45 44.811 11.63 58 176 176 40 75.98 45.681 11.95 70 176 180 43 81.19 49.091 10.85 64 162 170 44 81.42 39.442 13.08 63 174 176 38 81.87 60.055 8.63 48 170 186 44 73.03 50.541 10.13 45 168 168 45 87.66

19、37.388 14.03 56 186 192 45 66.45 44.754 11.12 51 176 176 47 79.15 47.273 10.60 47 162 164 54 83.12 51.855 10.33 50 166 170 49 81.42 49.156 8.95 44 180 185 51 69.63 40.836 10.95 57 168 172 51 77.91 46.672 10.00 48 162 168 48 91.63 46.774 10.25 48 162 164 49 73.37 50.388 10.08 67 168 168 57 73.37 39.407 12.63 58 174 176 54 79.38 46.080 11.17 62 156 165 52 76.32 45.441 9.63 48 164 166 50 70.87 54.625 8.92 48 146 155 51 67.25 45.118