随机模型方法及应用1

一元线性回归

第一节 大数定律与数理统计的若干知识

§1﹒1 大数定律及中心极限定理

大数定律（low of large numbers）及中心极限定理（central limit theorem）不仅为概率论（theary of probability）提供统计方面的理论保证，而且也为数理统计（mathematical statistics）的理论和方法奠定了坚实的理论基础。

1﹒1﹒1 ЧебЫШв不等式

设随机变量的方差D存在且有限，则对0，有

PED2 （２．１）

1﹒1﹒2 Bernoulli大数定律

n重实验中事件A出现的频率实验中出现的概率p，即0，

vn，依概率收敛于事件A在每次nvnlimPp1 （２．２） nn1﹒1﹒3 ЧебЫШв大数定律

设n是相互的随机变量序列，且

nDc，n1，2，

（２．３）

11

其中c是常数，则对0，有

1n1n limPkEk1 （２．４）

nnk1nk11﹒1﹒4 Хинчин大数定律

设n是相互的随机变量序列，且

En，n1，2， 则对0，有

n1nnlimPk1 k11﹒1﹒5 Lèvy－Lindeberg中心极限定理

设n是、同分布的随机变量序列，且

E，D2nn0，n1，2， 则xR，有

nknt2limPk1x1nn2x2edt １﹒１﹒6 De Moivre-Laplace中心极限定理

设n是、同分布的随机变量序列，且

n~B1，p，(0p1，n1，2，) 则xR，有

（２．７）

12

（２．５）

（２．６）（２．８）（２．９）nnpk1klimPxnnp1p12xet22 dt （２．１０）

§1﹒2 基本统计量和常用统计分布

在数理统计中，统计量（statistic）及其分布被广泛用于参数估计（parameters estimation）和假设检验等统计推断（statistical inference）的过程中，

1﹒2﹒1

统计量的定义及常用统计量

定义２．１ 设1，2，，n是总体的一个样本（sample），

T1，2，，是样本1，2，，n的不含任何未知参数的函数，则称T1，2，，为一个统计量；如果x1，x2，，xn是样本1，2，，n的一个观测值，那么称Tx1，x2，，xn是统计量

T1，2，，n的一个观测值。

定义２．２ 设1，2，，n是来自总体的一个容量为n的样本，常用的统计量有

１、 样本均值（sample mean）：

1nk （２．１１）

nk1MATLAB: mean(x)

２、 样本方差（sample variance）：

1nSknk122 （２．１２）

13

３、 样本标准差（sample standard deviation）：

S1nnk1k2 （２．１３）

４、 修正的样本方差（repaired sample variance）：

S*21nkn1k12 （２．１４）

MATLAB: var(x)

５、修正的样本标准差（repaired sample standard deviation）：

S*1nkn1k12 （２．１５）

1．２．２ 常用统计分布

１、2分布：设随机变量（random variable）1，2，，n相互、同分布，且k~N0，1，k1，2，，n，则随机变量

2k1n2k （２．１６）

所服从的分布称为自由度是n的2分布，记做2~2n。 如果随机变量2~2n，那么有 １）2的概率密度为：

0，x0nx1122xe，x0 （２．１７） pxnn222MATLAB: chi2cdf(x,n) 并且有

14

E2n （２．１８） D22n （２．１９）

x2 Distribution Density0.160.140.120.10.080.060.040.020from left side to right side: n = 5, 10, 20, 5001020304050607080 图２．１

２）定理２．１ （2分布的可加性）设随机变量k2~2nk，

k1，2，，n，则

n~nk （２．２０）

k1k122k2n３）定理２．２ 如果随机变量2~2n，那么

2n2nN0，1

L （２．２１）

４）定理２．３（Fisher） 如果随机变量2~2n，那么

15

22n1N0，1 （２．２２）

2L２、t分布：设随机变量，相互，且~N0，1，~2n，则随机变量

t （２．２３） n所服从的分布称为自由度是n的t分布，记做t~tn，并称其为自由度为n的t变量。

１） 设随机变量t~tn，则其概率密度为：

n1n1222x，x （２．２４） px1nnn2t Distribution Density0.40.350.30.250.20.150.10.050-5------- n = 5------- n = 50-4-3-2-1012345 图２．２

２）定理２．４ 设随机变量t~tn，则

16

1x2limpxe， x （２．２５） ntn2３、F分布：设随机变量，相互，且~2m，~2n，则随机变量

Fn （２．２６） m2称为自由度为m，n的F变量，所服从的分布称为的F分布，记做

F~Fm，n。

１）如果随机变量F~Fm，n，那么

1~Fn，m。 FF Distribution Density21.81.61.41.210.80.60.40.2000.511.522.5------- m = 10, n = 50------- m = 5, n = 1033.544.55------- m = 50, n = 200 图２．３

17

第二节 一元线性回归的若干问题

§2﹒1 简单线性回归分析

设随机变量y与随机变量x之间存在某种相关关系，对于x的取定的一组不完全相同的值x1,x2,xn，作实验得到n对观察结果：

xk,yk，k1,2,,n

（２７）

其中yk是随机变量y在xxk时的观测结果。

２﹒１﹒１ 简单线性回归模型及其基本理论假设

假设变量y与自变量x之间的相关关系可由下式表示：

yabx （２８）

其中~N0,回归模型。

由（２７）、（２８）可得

2，a和b是未知（回归）参数，称式（２８）为一元

ykabxkk （２９）

其中k~N0,2，且相互。

ˆ，ˆ和b当利用样本xk,yk，k1,2,,n，得到参数a和b的估计a 那么对于给定的x，取

ˆx （３０） ˆaˆby作为abx的估计，并称式（３０）为y关于x的线性回归方程，其图形称为回归直线。

２﹒１﹒２ 简单线性回归模型的基本特征

１、 由

ykabxkk，知yk是随机变量；

18

２、 EykEabxkkabxkEkabxk； ３、 DykDabxkkDk2； ４、 Covi,j0Covyi,yj0； ５、 ykEykykabxkk； ６、 yk~Nabxk,2。

２﹒１﹒３ 回归参数的最小二乘估计

１、最小二乘估计准则：

ˆminyabx2 （３１）ˆ,b Qa kka,bk1n５、 回归参数的最小二乘估计

ˆxaˆybnxkyknxy， （３２）ˆk1bn2xknx2k11n1n其中xxk，yyk。

nk1nk1３、回归参数的最小二乘估计的统计特性

ˆ都是y的线性组合；ˆ和b１） 线性性：akˆb ˆa和Eb２） 无偏性：Eaˆ的最小二乘估计都是a和b的所有线性无偏 ˆ和b３） 方差最小性：a 估计中方差最小的。

４） 2的估计、可决系数与相关系数定义：

 19

总偏差平方和STyky

2k1nnˆk，误差平方和 SEyky2k1nˆky回归平方和 SRy2k1总偏差平方和的分解：STSESR由于ESEn22，所以MSE称为平均误差平方和，MSR定义：

SE是2的一个无偏估计。MSEn2SR称为平均回归平方和。12SRSR （３３） STSRSE为可决系数，

2 （３４） 或

xk1nkxyky2xk1nkxyk1n （３５）

ky2为相关系数。

５）回归效果的显著性检验与方差分析表

由最小二乘法求得的线性回归方程是否具有实用价值，需要通过假设检验才能确定。如果线性假设符合实际，则b不应为零，因此，需要检验的假设为：

20

H0:b0，H1:b0 （１） F检验法 采用统计量：FMSR，当H0:b0为真时，F~F1,n2； MSE对于给定的显著性水平，如果FF11,n2，则应拒绝H0，认为线性回归效果显著；如果FF11,n2，则应接受H0，即认为线性回归效果不显著。

这一分析过程可由方差分析表给出：

方差分析表

误差来源 回归R 误差E 总和T 自由度 1 平方和S 均方和MS F SR MSR n2 n1 SE ST 表１ MSE MSRMSE （２） t检验法：

采用统计量

tˆbMSElxx， （３６）

21

其中lxxxkx，当H0:b0为真时，t~tn2；对于给定的

2k1n显著性水平，如果ttn2，则应拒绝H0，认为线性回归效果显

2著；如果ttn2，则应接受H0，即认为线性回归效果不显著。

2６） 回归参数的假设检验和参数估计 （１）回归参数b的假设检验和区间估计 记

ˆMSEs2blxxMSExk1n （３７）

2kx那么

ˆbb ~tn2 （３８）

ˆsbˆ0b假设检验H0:b0，H1:b0的统计量为t，因此，对于

ˆsb给定的显著性水平，如果ttn2，则应拒绝H0；否则接受H0。

2回归参数b的置信度为1001%的置信区间为：

ˆˆbtn2sb2 （３９） （３） 回归参数a的置信区间：回归参数a的置信度为1001%的 置信区间为：

aˆtn2saˆ （４０） 2

22

注：在做线性回归分析时，一般将分析结果记为：

ˆxˆaˆby （４１）

ˆsaˆsa７） 预测

ˆ0可由回归方计算，对于任何给定的xx0，y0abx0的点估计y在小样本情况下，y0的置信度为1001%的置信区间为：

21x0xyˆtn2MSE101nlxx2 （４２） 大样本时，y0的置信度为1001%的置信区间为：

yˆ0zMSE （４３） 12§１.３常

双曲线（Hyperbola）：

1bxay （４４） yxaxb１） 对数曲线（Logarithm Curve）

yablnx （４５）

２） 多项式曲线（Polynomial Curve）

ya0a1xamxm （４６） ３） 指数曲线（Exponent Curve）

yaex （４７）

５） S型曲线

23

y1 （４８） xabe上述曲线图示如下：

§１.３问

某建材实验室在作陶粒混凝土强度实验中，考察每立方米混凝土的水泥用量xkg对28天后的混凝土抗压强度ykg/cm2的影响，并测得如下数据：

xi150160170180190200210220230240250260 yi56.958.361.6.668.171.374.177.480.282.686.4.7 表２ （１） 求y关于x的线性回归方程，并问：每立方米混凝土中增加1kg水泥时，可提高的抗压强度是多少？

（２） 检验线性回归效果的显著性0.05； （３） 求回归参数b的区间估计10.95； （４） 求x022.5kg时，y的预测值与预测区间。

ˆ0.30399x10.28，3 回归直线方程： yˆ20.2393 总体方差点估计： 总体方差区间估计： 0.303993.7283105 回归效果显著性检验：

１、利用F检验法：F5.52251034.96F0.951,n2 ２、利用t检验法：t74.31372.2281t0.975n2

 24

x022.5kg时，y的预测值: yˆx22.516.3626

预测区间: 16.36262.0323 Linear Regression Model90 85 80 75 70----------- Regression Beeline 65 60 55 140160180200220240260图2.4

Linear Regression Model100

90

80 70 60----------- Regression Beeline 50 40 30 20

10050100150200250300

25 图2.5

附程序：Linear_Regression_Model.m

%Linear Regression Model clear all x=150:10:260;

y=[56.9 58.3 61.6 .6 68.1 71.3 74.1 77.4 80.2 82.6 86.4 .7]; L=length(x);

[P,S]=polyfit(x,y,1); z=P(1).*x+P(2); fity=polyval(P,x); Prey=polyval(P,22.5) plot(x,y,'r*',x,fity) Meanx=mean(x); Meany=mean(y); Lxx=(L-1)*var(x); St=(L-1)*var(y); Sr=0; Se=0; for k=1:L

Se=Se+(y(k)-z(k))^2; Sr=Sr+(z(k)-Meany)^2; end

MSe=Se/(L-2); r=Sr/St

F=Sr/MSe % if P(1)=0 F~F(1,n-2). Sb=MSe/Lxx

T=abs(P(1)/(MSe/Lxx)^.5) % if P(1)=0 T~t(n-2). Falgfa=finv(0.95,1,L-2) Talgfa=tinv(0.975,L-2) Nalgfa=norminv(0.975,0,1)

title('Linear Regression Model')

gtext('----------- Regression Beeline') hold on

%Solving the prediction Interval t=20:.5:260;

26

prey=polyval(P,t); if L<45

prey1=prey-Talgfa*(MSe*(1+1/L+(t-Meanx).^2/Lxx)).^.5; prey2=prey+Talgfa*(MSe*(1+1/L+(t-Meanx).^2/Lxx)).^.5; else

prey1=prey-Nalgfa*(MSe)^.5; prey2=prey+Nalgfa*(MSe)^.5; end

plot(t,prey,t,prey1,'r',t,prey2,'r') P=poly2sym(P); P=vpa(P,5)

%H=polytool(x,y,1,0.05,22.5)

1.50 1.65 1.80 1.95 2.10 2.25 2.40 2.55 2.70 2.85 3.00 3.15 3.30 3.45 3.60 3.75 3.90 4.05 4.20 4.35 4.50 4.65 4.80 4.95 5.10 5.25 5.40 5.55 5.70 5.85

1.77 2.07 2.26 2.41 2.61 2.69 2.87 2.75 2.98 3.08 3.04 3.32 3.28 3.62 3.48 3.70 3.74 3.71 3.65 3.93 3.91 4.14 4.19 4.37 4.31 4.25 4.41 4.33 4.48 4.53 4.5 4Nonlinear Regression Model53.527 4.5342.53.5322.5----------- Regression Curve

Nonlinear Regression Model

76543----------- Regression Curve

210-1024681012ˆ1.05661.9553lnx， 回归曲线方程： yˆ20.2393 总体方差点估计： 总体方差区间估计： 0.303993.7283105 回归效果显著性检验：

28

１、利用F检验法：F5.52251034.96F0.951,n2 ２、利用t检验法：t74.31372.2281t0.975n2

ˆx022.5kg时，y的预测值: yx22.516.3626

预测区间: 16.36262.0323 1.9553*x+1.0566

附程序：Nonlinear_Regression_Model01.m

%Linear Regression Model clear all

s=1.5:0.15:5.85;

y=[1.77 2.07 2.26 2.41 2.61 2.69 2.87 2.75 2.98 3.08 3.04 3.32 3.28 3.62 3.48...

3.70 3.74 3.71 3.65 3.93 3.91 4.14 4.19 4.37 4.31 4.25 4.41 4.33 4.48 4.53];

L=length(s); x=log(s);

[P,S]=polyfit(x,y,1); LP=length(P); z=zeros(1,L); for i=1:LP

z=z+P(i).*x.^(LP-i); end

fity=polyval(P,x); plot(s,y,'r*',s,fity)

title('Nonlinear Regression Model') gtext('----------- Regression Curve') Meanx=mean(x); Meany=mean(y); Lxx=(L-1)*var(x); St=(L-1)*var(y); Sr=0;Se=0; for k=1:L

29

Se=Se+(y(k)-z(k))^2; Sr=Sr+(z(k)-Meany)^2; end Se=Se

MSe=Se/(L-2)

F=Sr/MSe % if P(1)=0 F~F(1,n-2). T=abs(P(1)/(MSe/Lxx)^.5) % if P(1)=0 T~t(n-2). Falgfa=finv(0.95,1,L-2) Talgfa=tinv(0.975,L-2) hold on l=0.5:.5:12; t=log(l);

prey=polyval(P,t); if L<45

prey1=prey-Talgfa*(MSe*(1+1/L+(t-Meanx).^2/Lxx)).^.5; prey2=prey+Talgfa*(MSe*(1+1/L+(t-Meanx).^2/Lxx)).^.5; else

prey1=prey-Nalgfa*(MSe)^.5; prey2=prey+Nalgfa*(MSe)^.5; end

plot(l,prey,l,prey1,'r',l,prey2,'r') P=poly2sym(P); P=vpa(P,5)

附程序：Nonlinear_Regression_Model02.m

%Linear Regression Model clear all

x=1.5:0.15:5.85;

y=[1.77 2.07 2.26 2.41 2.61 2.69 2.87 2.75 2.98 3.08 3.04 3.32 3.28 3.62 3.48...

3.70 3.74 3.71 3.65 3.93 3.91 4.14 4.19 4.37 4.31 4.25 4.41 4.33 4.48 4.53] L=length(x);

30

n=input('The Degree n of The Polynomial Pn(x) n = '); [P,S]=polyfit(x,y,n); LP=length(P); z=zeros(1,L); for i=1:LP

z=z+P(i).*x.^(LP-i); end

fity=polyval(P,x); plot(x,y,'*',x,fity) Meany=mean(y); Lxx=(L-1)*var(x); St=(L-1)*var(y); Sr=0;Se=0; for k=1:L

Se=Se+(y(k)-z(k))^2; Sr=Sr+(z(k)-Meany)^2; end Se=Se

MSe=Se/(L-2)

F=Sr/MSe T=abs(P(1)/(MSe/Lxx)^.5) Falgfa=finv(0.95,1,L-2) Talgfa=tinv(0.975,L-2) P=poly2sym(P); P=vpa(P,5)

% if P(1)=0 F~F(1,n-2). % if P(1)=0 T~t(n-2). 31