收藏 2019-10-17

问题1（随便给一组二维数据）

完整代码如下

setwd("E:/AllClass/junior1/RegressionAnalysis/unit2")#设置文件路径
#保留原路径setwd("C:/User/10854/documents")
#以下利用理论方法使用一元回归模型
#导入数据
data<-read.csv("2-7.csv")#书本2.15，表数据2-7
x<-data[,1]
y<-data[,2]
n<-length(x)
split.screen(c(1,3))
screen(1)
plot(x,y,pch=16)
title(main="数据散点图")
#求均值与回归变量lxx，lyy，lxy
meanx<-mean(x)
meany<-mean(y)
lxx<-sum((x-meanx)^2)
lyy<-sum((y-meany)^2)
lxy<-sum((x-meanx)*(y-meany))
#回归系数估计
beta_1<-lxy/lxx#beta_1
beta_0<-meany-beta_1*meanx#beta_0
screen(2)
plot(x,y,pch=16)
points(x,beta_0+beta_1*x,type="l")
title(main="回归图")
#预测值与平方和
y_hat<-beta_0+beta_1*x
sse<-sum((y_hat-y)^2)#残差平方和
ssr<-sum((y_hat-meany)^2)#回归平方和
sst<-ssr+sse#总离差平方和
#回归误差ε的方差sigma估计
sigma_hat<-sqrt(1/(n-2)*sse)
#对bet_0、beta_1的95%区间估计
alpha<-0.05
#beta_0,beta_1的分布标准差
sd.beta_0<-sqrt((1/n+(meanx^2)/lxx))*sigma_hat
sd.beta_1<-sqrt(sigma_hat^2/lxx)
beta_1l<-beta_1-qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信下限
beta_1u<-beta_1+qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信上限
beta_0l<-beta_0-qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信下限
beta_0u<-beta_0+qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信上限
#remark
#qt是求出置信度1-α对应的统计量值t（1-α）
#dt是求出统计量对应的置信度值，即p值(这里用不上t分布)
#dt返回概率密度，pt返回分布函数，qt返回分位数函数，rt生成随机数。
#qf\df都是同理，对应的是F分布
#计算xy决定系数
R<-ssr/sst
#回归方程的显著性检验
#法一方差分析F检验
f<-(ssr/1)/(sse/(n-2))
p1<-pf(f,1,n-2)#F为统计量、1为第一个自由度，n-2为第二个自由度
#法二回归系数的beta_1的t检验
t1<-beta_1/sd.beta_1#t统计量
p2<-pt(t1,n-2)
#法三相关系数r的t检验
r<-lxy/(sqrt(lxx*lyy))
t2<-sqrt(n-2)*r/sqrt(1-r^2);
p3<-pt(t2,n-2) #p值
#残差图
screen(3)
e<-y_hat-y  #残差
n<-length(e) 
sigma_u<-seq(2*sigma_hat,2*sigma_hat,length.out=n)  #残差2σ原则
sigma_l<-seq(-2*sigma_hat,-2*sigma_hat,length.out=n)
plot(x,e,pch=16,ylim=c(5,-5))
points(x,sigma_u,type="l")  #画2σ上下界
points(x,sigma_l,type="l")
title(main="残差图")
#预测广告费用为1000万元时，销售达多少
x0<-1000
y0<-beta_0+beta_1*x0
#因变量新值得95%置信区间
h00<-1/n+((x0-meanx)^2)/lxx
y0_l<-y0-qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信下限
y0_u<-y0+qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信上限
#近似置信区间
y0_l_ <- y0-2*sigma_hat
y0_u_ <- y0+2*sigma_hat
#因变量新值得平均值的95%置信区间
y0_l_E<-y0-qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信下限
y0_u_E<-y0+qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信上限
#以下利用R函数回归

(1)画散点图

1.1问题求解

1.1.1输入

#导入数据
data<-read.csv("2-7.csv")#书本2.15，表数据2-7
x<-data[,1]
y<-data[,2]
n<-length(x)
# split.screen(c(1,3)) 
# screen(1)
plot(x,y,pch = 16)
title(main="数据散点图")

1.1.2输出

(2) ${x}$ 与 $y$ 之间是否大致呈线性关系

由第一问的散点图，大致正相关，并且呈线性关系

(3)用最小二乘估计求回归方程

3.1问题分析

先求样本均值
$\bar{x}=\Sigma{x_i},\bar{y}=\Sigma{y_i}$

$L_{xx}=\Sigma_{i=1}^{n}(x_i-\bar{x})^2$

$L_{yy}=\Sigma_{i=1}^{n}(y_i-\bar{y})^2$

$L_{xy}=\Sigma_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})$

则回归系数的估计如下：
$\hat\beta_1=\frac{L_{xx}}{L_{yy}}$

$\hat\beta_0=\bar{y}-\beta_1\bar{x}$

得到回归方程如下：
$y_i=\hat\beta_0+\hat\beta_1x_i$

3.2问题求解

3.2.1输入

#求均值与回归变量lxx，lyy，lxy
meanx<-mean(x)
meany<-mean(y)
lxx<-sum((x-meanx)^2)
lyy<-sum((y-meany)^2)
lxy<-sum((x-meanx)*(y-meany))
#回归系数估计
beta_1<-lxy/lxx#beta_1
beta_0<-meany-beta_1*meanx#beta_0
# screen(2)
plot(x,y,pch=16)  #散点
points(x,beta_0+beta_1*x,type="l") #回归线
title(main="回归图")

3.2.2输出

得出各参数

统计量	统计值
$\bar{x}$	762
$\bar{y}$	2.85
$L_{XX}$	1297860
$L_{xy}$	$4653$
$L_{yy}$	18.525
$\hat\beta_0$	0.118
$\hat\beta_1$	0.00359

回归方程：
$y_i=0.118+0.00359x_i$

回归图：

(4)求回归标准误差 $\hat\sigma$

4.1问题分析

分别得出
$SSE=\Sigma(y_i-\hat{y})^2$
标准误差 ${\sigma}$ 的无偏估计为：
$\hat\sigma={\frac{1}{n-2}}\Sigma_{i=1}^{n}e_{i}^2={\frac{1}{n-2}}\Sigma\left(y_i-\hat{y}\right)^2$

4.2问题求解

4.2.1输入

sse<-sum((y_hat-y)^2)#残差平方和
#回归误差ε的方差sigma估计
sigma_hat<-sqrt(1/(n-2)*sse)

4.2.2输出

$SSE=1.843,\hat\sigma=0.48$

(5)给出 $\hat\beta_0$ 与 $\hat\beta_1$ 的置信度为95%的区间估计

5.1问题分析

$\hat\beta_1$ 与 $\hat\beta_0$ 的分布为
$\hat\beta_1\sim{N\left(\beta_1,\frac{\sigma^2}{L_{xx}}\right)} , \quad \hat\beta_0\sim{N\left(\beta_0,\left[\frac{1}{n}+\frac{\bar{x}^2}{L_{xx}}\right]\sigma^2\right)}$
由 $\hat\beta_1$ 分布构造了服从自由度为 $n-2$ 的 $t$ 分布统计量
$t=\frac{(\hat\beta_1-\beta_1)\sqrt{L_{xx}}}{\hat\sigma}$
因而
$P\left(\left|\frac{(\hat\beta_1-\beta_1)\sqrt{L_xx}}{\hat\sigma}\right|>t_{\alpha/2}\left(n-2\right)\right)$
得到 $\beta_1$ 的置信度为 $1-\alpha$ 的置信区间为( $\alpha=0.05$ )
$\left( \hat\beta_1-t_{\alpha/2}\frac{\hat\sigma}{\sqrt{L_{xx}}}, \hat\beta_1+t_{\alpha/2}\frac{\hat\sigma}{\sqrt{L_{xx}}} \right)$
对 $\hat\beta_0$ 同理得置信度为 $1-\alpha$ 的置信区间为( $\alpha=0.05$ )
$\left( \hat\beta_1-t_{\alpha/2}\sqrt{\left[\frac{1}{n}+\frac{\bar{x}^2}{L_{xx}}\right]}\hat\sigma, \hat\beta_1+t_{\alpha/2}\sqrt{\left[\frac{1}{n}+\frac{\bar{x}^2}{L_{xx}}\right]}\hat\sigma \right)$

5.2问题求解

5.2.1输入

alpha<-0.05 #置信度
#beta_0,beta_1的分布标准差
sd.beta_0<-sqrt((1/n+(meanx^2)/lxx))*sigma_hat
sd.beta_1<-sqrt(sigma_hat^2/lxx)
beta_1l<-beta_1-qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信下限
beta_1u<-beta_1+qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信上限
beta_0l<-beta_0-qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信下限
beta_0u<-beta_0+qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信上限

5.2.2输出

得到 $\beta_1$ 的置信度为0.05的置信区间为 $\left[0.0026,0.0046 \right]$

得到 $\beta_1$ 的置信度为0.05的置信区间为[-0.7,0.937]

(6) $x$ 与 $y$ 的决定系数

6.1问题分析

各平方和如下：
$SSE=\Sigma(y_i-\hat{y})^2$

$SSR=\Sigma(\hat{y}-\bar{y})^2$

$SST=SSR+SSE$

决定系数为
$r^2=\frac{SSR}{SST}$

6.2问题求解

6.2.1输入

sse<-sum((y_hat-y)^2)#残差平方和
ssr<-sum((y_hat-meany)^2)#回归平方和
sst<-ssr+sse#总离差平方和
#计算xy决定系数
R<-ssr/sst

6.2.2输出

$SSE=1.84,SSR=16.68,SST=18.525,r^2=0.9$

(7)对回归方程做方差分析

7.1问题求解

7.1.1输入

f<-(ssr/1)/(sse/(n-2))
p1<-pf(f,1,n-2)#F为统计量、1为第一个自由度，n-2为第二个自由度
#1-p1为p值

7.2.2输出

一元线性回归方差分析表如下：

方差来源	自由度	平方和	均方	$F$ 值	$P$ 值
回归	1	$SSR=16.68$	$SSR/1=16.68$	$\frac{SSR/1}{SSE(n-2)}=72.40$	$p=2.8*10^{-5}$
残差	$n-2$	$SSE=1.84$	$SSE/(n-2)=0.23$
总和	$n-1$	$SST=18.525$

(8)做回归系数 $\beta_1$ 的显著性检验

8.1问题分析

$\hat\beta_1$ 的分布为
$\hat\beta_1\sim{N\left(\beta_1,\frac{\sigma^2}{L_{xx}}\right)}$
假设检验原假设为 $H_0:\beta_1=0$ ,构造检验统计量 $t$ 服从自由度为 $n-2$ 的 $t$ 分布
$t=\frac{\hat\beta_1\sqrt{L_{xx}}}{\hat\sigma}$
并计算对应的 $p$ 值

8.2问题求解

8.2.1输入

#回归系数的beta_1的t检验
t1<-beta_1/sd.beta_1  #t统计量
p2<-pt(t1,n-2) #（1-p2）为p值

8.2.2输出

$t=8.50,p=1.40*10^{-5}$

(9)做相关系数的显著性检验

9.1问题分析

相关系数为 $r$
$r=\frac{L_{xy}}{\sqrt{L_{xx}L_{yy}}}=\hat\beta_1\sqrt{\frac{L_{xx}}{L_{yy}}}$
构造检验统计量
$t=\frac{{\sqrt{n-2}}{r}}{\sqrt{1-r^2}}\sim{t(n-2)}$
当 $\vert{t}\vert>t_{\alpha/2}(n-2)$ 时，认为简单回归系数显著不为零

9.2问题求解

9.2.1输入

#相关系数r的t检验
r<-lxy/(sqrt(lxx*lyy))
t2<-sqrt(n-2)*r/sqrt(1-r^2);
p3<-pt(t2,n-2) #（1-p3）为p值

9.2.2输出

$t=8.50,p=1.40*10^{-5}$

(10)对回归方程作残差图并做相应的分析

10.1问题求解

10.1.1输入

#残差图
# screen(3)
e<-y_hat-y  #残差
n<-length(e) 
sigma_u<-seq(2*sigma_hat,2*sigma_hat,length.out=n)  #残差2σ原则
sigma_l<-seq(-2*sigma_hat,-2*sigma_hat,length.out=n)
plot(x,e,ylim=c(5,-5))
points(x,sigma_u,type="l")  #画2σ上下界
points(x,sigma_l,type="l")
title(main="残差图")

10.1.2输出

10.2问题分析

计算出残差后，及自变量x为横轴，以残差为总纵轴画散点图得到残差图，从残差图上看出，产茶时围绕着 $e=0$ 随机波动，并且波动范围在方差估计 $\hat\sigma$ 的两倍。所以可以判定模型的基本假定是满足的。

(11)该公司预计下一周签发新保单 $x_0=1000$ 张，需要的加班时间时多少？

11.1问题求解

11.1.1输入

x0<-1000
y0<-beta_0+beta_1*x0

11.1.2输出

$y=3.7$

(12)给出 $y_0$ 的置信度为95%的精确预测区间和近似预测区间

12.1问题分析

可以计算得到 $\hat{y}_0$ 的分布为
$\hat{y}_0\sim N\left( \beta_0+\beta_1x_0,(\frac{1}{n}+\frac{(x_0-\bar{x})^2}{L_{xx}})\sigma^2\right)$
所以枢轴量为
$y_0-\hat{y}_0\sim N\left(0,(1+h_{00})\sigma^2\right),h_{00}=(\frac{1}{n}+\frac{(x_0-\bar{x})^2}{L_{xx}})$
统计量为
$t=\frac{y_0-\hat{y_0}}{\sqrt{1+h_{00}}\hat\sigma}\sim t(n-2)$
得精确置信区间为
$\hat{y_0} \pm t_{\alpha/2}(n-2)\sqrt{1+h_{00}}\hat\sigma$
近似预测区间为
$\hat{y}_0\pm2\hat\sigma$

12.2问题求解

12.2.1输入

#因变量新值95%置信区间
h00<-1/n+((x0-meanx)^2)/lxx
y0_l<-y0-qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信下限
y0_u<-y0+qt(1-alpha/2,n-2)*sqrt(1+h00)*sigma_hat#预测新值得置信上限

12.2.2输出

精确置信区间为[2.519,4.887]

近似预测区间为[2.743,4.663]

所以两区间比较接近，可以用近似区间估计

(13)给出 $E(y_0)$ 的置信度为95%的区间估计

13.1问题分析

根据 $\hat{y}_0$ 构造枢轴量（含 $E(y_0)$ ）
$\hat{y_0}-E（{y}_0）\sim N\left(0,h_{00}\sigma^2\right),h_{00}=(\frac{1}{n}+\frac{(x_0-\bar{x})^2}{L_{xx}})$
统计量为
$t=\frac{y_0-\hat{y_0}}{\sqrt{h_{00}}\hat\sigma}\sim t(n-2)$
得置信水平95%得精确置信区间为
$\hat{y_0} \pm t_{\alpha/2}(n-2)\sqrt{h_{00}}\hat\sigma$

13.2问题求解

13.2.1输入

#因变量新值得平均值的95%置信区间
y0_l_E<-y0-qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信下限
y0_u_E<-y0+qt(1-alpha/2,n-2)*sqrt(h00)*sigma_hat#预测新值均值得置信上限

13.2.2输出

$E(y_0)$ 的95%置信区间为[3.283,4.122]

例题2.1(火灾)

（1）问题求解

根据上一题的程序，我们对例题2.1直接求解

1.1输入

setwd("E:/AllClass/junior1/RegressionAnalysis/unit2")#设置文件路径
#保留原路径setwd("C:/User/10854/documents")
#以下利用理论方法使用一元回归模型
#导入数据
data<-read.csv("2-1.csv")#书本火灾题目，表数据2-1
x<-data[,1];y<-data[,2]
n<-length(x)
split.screen(c(1,3))
screen(1)
plot(x,y)
title(main="数据散点图")
#求均值与回归变量lxx，lyy，lxy
meanx<-mean(x);meany<-mean(y)
lxx<-sum((x-meanx)^2)
lyy<-sum((y-meany)^2)
lxy<-sum((x-meanx)*(y-meany))
#回归系数估计
beta_1<-lxy/lxx#beta_1
beta_0<-meany-beta_1*meanx#beta_0
screen(2)
plot(x,y)
points(x,beta_0+beta_1*x,type="l")
title(main="回归图")
#预测值与平方和
y_hat<-beta_0+beta_1*x
sse<-sum((y_hat-y)^2)#残差平方和
ssr<-sum((y_hat-meany)^2)#回归平方和
sst<-ssr+sse#总离差平方和
#回归误差ε的方差sigma估计
sigma_hat<-sqrt(1/(n-2)*sse)
#对bet_0、beta_1的95%区间估计
alpha<-0.05
#beta_0,beta_1的分布标准差
sd.beta_0<-sqrt((1/n+(meanx^2)/lxx))*sigma_hat
sd.beta_1<-sqrt(sigma_hat^2/lxx)
beta_1l<-beta_1-qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信下限
beta_1u<-beta_1+qt(1-alpha/2,n-2)*sd.beta_1#beta_1置信上限
beta_0l<-beta_0-qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信下限
beta_0u<-beta_0+qt(1-alpha/2,n-2)*sd.beta_0#beta_0置信上限
#计算xy决定系数
R<-ssr/sst
#回归方程的显著性检验
#方差分析F检验
f<-(ssr/1)/(sse/(n-2))
p1<-pf(f,1,n-2)#F为统计量、1为第一个自由度，n-2为第二个自由度
#回归系数的beta_1的t检验
t1<-beta_1/sd.beta_1#t统计量
p2<-pt(t1,n-2)
#残差图
screen(3)
e<-y_hat-y
n<-length(e)
sigma_u<-seq(2*sigma_hat,2*sigma_hat,length.out=n)#残差2σ原则
sigma_l<-seq(-2*sigma_hat,-2*sigma_hat,length.out=n)
plot(x,e,ylim=c(12,-12))
points(x,sigma_u,type="l")
points(x,sigma_l,type="l")
title(main="残差图")

1.2输出

回归方程如下
$y_i=10.28+4.91x_i$
参数表如下：

参数	参数值
$\hat\beta_0$	10.28
$\hat\beta_1$	4.91
$r^2$	0.92
$\hat\sigma$	2.31
$\hat\beta_0区间估计$	[7.2,13.34]
$\hat\beta_1区间估计$	[4.07,5.76]