统计软件应用课程设计

Monte Carlo 在线性回归中的简单应用

班级：1324001 姓名：徐萍萍 学号：201320400118 指导老师：闫霏霏 电话号码：

2016年4月17日

摘要：本文针对确定运动员的耗氧量与其他一些因素的关系数据，在

SAS中进行

回归分析得到对数据拟合较好的线性模型。再用Monte Carlo随机过程产生残差项并代入线性方程中，分别假定残差项不符合均值为零、正态分布、异方差三大假设，从而检验回归理论。

关键字：蒙特卡罗、多元线性回归、残差项基本假设

问题重述：在运动生理学的研究中，为了确定运动员的耗氧量与其他一些因素的

关系，在一个实验中对31个人测量了年龄（age），体重（weight），跑完1.5英里用的时间（runtime），静态时的心率（rstpulse）,跑动时的心率（runpulse）,跑步时的最大心率（maxpulse）,每公斤体重每分钟的耗氧量（oxy）。实测数据（oxy.txt）见下表，试以oxy为因变量，估计该变量对于问题中所有其他变量的直线回归方程。并用Monte Carlo验证回归理论中残差项三大基本假设。(数据在附录中)

第一章、基本理论

一、蒙特卡罗方法

㈠蒙特卡罗（Monte Carlo）方法概述

蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用蒙特卡罗命名。

㈡蒙特卡罗方法的基本原理

由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡罗法正是基于此思路进行分析的。 设有统计的随机变量Xi(i=1，2，3，…，k)，其对应的概率密度函数分别为fx1，fx2，…，fxk，功能函数式为Z=g(x1，x2，…，xk)。 首先根据各随机变量的相应分布，产生N组随机数x1，x2，…，xk值，计算功能函数值 Zi=g(x1，x2，…，xk)(i=1，2，…，N)，若其中有L组随机数对应的功能函数值Zi≤0，则当N→∞时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。

㈢蒙特卡罗方法的收敛性

设所求的量是随机变量的数学期望E(x)，那么Monte Carlo方法通常使用随机变量的简单子样

，12，，  ， n的算术平均值，即

1N-

i1Ni

作为所求量X的近似值。由柯尔莫哥罗夫（Kolmogorov）大数定理可知，

P（lim-Nx）1

_N-即当N充分大时，有NE()x成立的概率等于1，亦即可以用

N作为所求量X的估计值。

根据中心极限定理，如果随机变量的标准差不为零，那么Monte Carlo方法的误差为

N

式中，

O()为正态差，是与置信水平有关的常量。Monte Carlo方法的收敛速度的阶为

N1/2，误差是由随机变量的标准差S和抽样次数N决定的。

二、多元线性回归模型

㈠多元线性回归模型的一般形式

设随机变量y与一般变量

x,x,...,x12p的线性回归模型为：

y01xx12013x

,,记为①式中，

是P+1个未知参数，

10称为回归常数，

p，

…

1

p称为回归系数。y称为被解释变量（因变量），

x,x,...,x2是P个可以精确测量

并控制的一般变量，称为解释变量（自变量）。P=1时，①式即一元线性回归模型； P≧2时，我们就称①式为多元线性回归模型。是随机误差，与一元线性回归一样，

0E（）2var() 对随机误差我们常规定：称

E(y)01xx122...pxp

n

组观测数据

为理论回归方程。

对一个实际问题，如果我们获得

(xi1,xi2,...,xip:)i1,2,...,n）y（，则线性回归模型①式可表示为：

iyxxx10111212p1p1y201x212x22px2p2xnpnpyn01xn12xn2

写成矩阵形式为：

yX

记位②式，式中

y1yy2yn1  11121p1   21222pX   1   n1n2npxxxxxxxxx

012 1np

X是一个n(p1)阶矩阵，称为回归设计矩阵或资料矩阵。在实验设计中X的元素是预先设定并可以控制的，人的主观因素可作用其中，因而称X为设计矩阵。

㈡多元线性回归模型的基本假定

回归方程①式有如下基本假定：

⑴解释变量x1,x2,…,xp是确定性变量，不是随机变量，且要求rank(X)=p+1,n E(i)0 , i1,2,2 ， iji,j1,2,,ncov(,)ij0 , ij

这个假定常称为高斯-马尔柯夫条件。E(i)0，即假设观测值没有系统误差，随机误差项i的平均值为零。随机误差项的协方差为零，表明随机误差项在不同的样本点之间是不相关的（在正态假定下即为的），不存在序列相关，并且有相同的

精度。

i~N(0 , 2) , i1,2,,n,,  , n 相互⑶正态分布的假设条件为：12

对于多元线性回归的矩阵模型②式，这个条件可表示为：

~N（0 ， 2In）

第二章、模型建立

一、SAS中线性回归过程

这里在SAS软件中，对问题中数据（oxy.txt）做多元线性回归，运行程序如下： I：libname ep 'e:\\saslx';run; data ep.oxys;

infile 'e:\\sasdt\\oxy.txt';

input age weight oxy runtime rstpulse runpulse maxpulse @@; proc reg;

model oxy= age weight 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的影响都有统计意义。

⑵ 决定系数R-square=0.8487，表示所得回归防城对数据的线性拟合很好。它的平方根为因变量oxy与6个自变量的复相关系数。

⑶ Parameter Estimates部分是检验常数项或自变量系数等于零的假设。由检验结果可知，当显著水平为0.05，体重weight(p=0.1869)和休息时的心率rstpulse(p=0.7473)这两个变量应删除。一般来说，某个自变量是否对因变量有统计意义，不仅决定于这个变量本身，还跟回归方程中同时存在的其他自变量有关。若把对因变量无统计意义的自变量留在回归方程，会增大参数估计或误差，因此这些变量应剔除。

由以上分析可知，应采用逐步回归法——STEPWISE过程，即采用II中程序。在逐步回归过程的结果表中，可以看到引进变量的顺序、每一步引进的变量对R-square的贡献和总的R-square的变化，C（P）值的变化及引进每个变量时的检验概率。其中C（P）值最小的最后那个模型为最适合的回归模型。最后得估计每公斤体重每分钟耗氧量（oxy）的线性方程为：

y=0.27051x1-0.34811x2-2.76758x3-0.19773x4+98.147

其中x1、x2、x3、x4依次代表maxpulse、runpulse、runtime、age,体重weight和休息时的心率rstpulse被剔除。 图一

二、蒙特卡罗方法回归检验

由Monte Carlo思想，由正态随机数产生误差项e代入线性方程中，方程中自变量系数为在回归理论中估计得出的系数。这里通过分别假定残差项不符合均值为零E(i)0、正态分布、异方差三大假设，从而检验回归理论。

⑴基于以上思想，在SAS中编写程序III如下： data ep.oxys;

infile 'e:\\sasdt\\oxy.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 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

选择其中x1的系数为例，即1与1进行比较，计算估计系数与给定的真实系数的相

ˆ1-10.270510.270530.27051ˆ对误差，则

10.0073%，可知误差较小，则利

用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.765*x3-0.19200*x4+e

ˆ1-10.270510.276980.270512.391%1

上式中以自变量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

ˆ1-10.270510.283430.270514.776%1

上式中同样以x1的系数为例，计算相对误差为4.776% 误差明显增大。即用符合均

匀分布的残差随机样本，得到的线性方程估计系数偏差较大。 图六

⑷将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.343*x2-2.802*x3-0.15246*x4+e

ˆ1-10.270510.257870.270514.672%1

计算出x1(maxpulse)系数估计值与真实值的相对误差为4.671%，误差较与符合所有

假设时的0.0073%有较大变化。 图七

第三章、结论与体会

综合以上分析，先利用数据在SAS中得到拟合较为精确的线性回归方程，再利用Monte Carlo作为一个工具用于产生随机数。改变模型，分别产生不符合线性回归三个基本假设的残差项随机样本运用于新的模型中，得到的估计系数相较与符合基本假设的估计系数相对误差有较大变化。由此说明，要得到精确的线性回归模型必须符合基本假设，即达到了检验回归理论的目的。Monte Carlo（随机模拟）还可应用于很多有关于随机现象的近似计算方面，它还可以应用于概率计算、求积分的数值近似解、模拟置信区间理论、模拟假设检验。 在这次课程设计过程中，首先很感谢老师的提点和同学们的帮助。自己也感觉受益良多，大到对Monte Carlo方法与回归理论的理解，小至论文的格式都进步不少，相信对将来毕业论文设计会有一定帮助。

参考文献

【1】何晓群、刘文卿 应用回归分析[M] 中国人民大学出版社，2015. 【2】董大钧 统计分析应用[M] 电子工业出版社，2014. 【3】林晓辰 蒙特卡罗方法 MBA智库百科，2008.

附录

age weight oxy runtime rstpulse runpulse maxpulse 44 .47 44.609 11.3762 178 182 40 75.07 45.313 10.07 62 185 185 44 85.84 54.297 8.65 45 156 168 42 68.15 59.571 8.17 40 166 172 38 .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 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 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 1 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 1 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 1 166 50 70.87 54.625 8.92 48 146 155 51 67.25 45.118 11.08 48 172 172 54 91.63 39.203 12.88 44 168 172 51 73.71 45.790 10.47 59 186 188 57 59.08 50.545 9.93 49 148 155 49 76.32 48.673 9.40 56 186 188 48 61.24 47.920 11.50 52 170 176 52 82.78 47.467 10.50 53 170 172

东华理工大学 学年课程设计报告评分表

学生姓名：徐萍萍学号：201320400118 班级：1324001 课程设计题目：Monte Carlo在线性回归中的简单应用 项目内容 满分 5 10 10 10 10 15 15 15 10 5 100 实 评 能结合所学课程知识、有一定的能力训练。符合选题要求 选 （3人一题） 题 工作量适中，难易度合理 能熟练应用所学知识，有一定查阅文献及运用文献资料能力 能 理论依据充分，数据准确，公式推导正确 力 水 能应用计算机软件进行编程、资料搜集录入、加工、排版、制平 图等 能体现创造性思维，或有独特见解 模型正确、合理，各项技术指标符合要求。 摘要叙述简练完整，假设合理、问题分析正确、数学用语准确、结论严谨合理；问题处理科学、条理分明、语言流畅、结构严谨、版面清晰 课程设计报告主要部分齐全、合理，符号统一、编号齐全。 格式、绘图、表格、插图等规范准确，符合课程设计报告要求 正文字数不少于2000字，不超过15000字 指导教师评语： 总 分 成 果 质 量 指导教师签名： 年 月 日