问题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)与之间是否大致呈线性关系
由第一问的散点图,大致正相关,并且呈线性关系
(3)用最小二乘估计求回归方程
3.1问题分析
先求样本均值
则回归系数的估计如下:
得到回归方程如下:
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输出
得出各参数
统计量 | 统计值 |
---|---|
762 | |
2.85 | |
1297860 | |
18.525 | |
0.118 | |
0.00359 |
回归方程:
回归图:
(4)求回归标准误差
4.1问题分析
分别得出
标准误差的无偏估计为:
4.2问题求解
4.2.1输入
sse<-sum((y_hat-y)^2)#残差平方和
#回归误差ε的方差sigma估计
sigma_hat<-sqrt(1/(n-2)*sse)
4.2.2输出
(5)给出与的置信度为95%的区间估计
5.1问题分析
与的分布为
由分布构造了服从自由度为的分布统计量
因而
得到的置信度为的置信区间为()
对同理得置信度为的置信区间为()
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输出
得到的置信度为0.05的置信区间为
得到的置信度为0.05的置信区间为[-0.7,0.937]
(6) 与的决定系数
6.1问题分析
各平方和如下:
决定系数为
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输出
(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输出
一元线性回归方差分析表如下:
方差来源 | 自由度 | 平方和 | 均方 | 值 | 值 |
---|---|---|---|---|---|
回归 | 1 | ||||
残差 | |||||
总和 |
(8)做回归系数的显著性检验
8.1问题分析
的分布为
假设检验原假设为 ,构造检验统计量服从自由度为 的 分布
并计算对应的 值
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输出
(9)做相关系数的显著性检验
9.1问题分析
相关系数为
构造检验统计量
当 时,认为简单回归系数显著不为零
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输出
(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为横轴,以残差为总纵轴画散点图得到残差图,从残差图上看出,产茶时围绕着 随机波动,并且波动范围在方差估计的两倍。所以可以判定模型的基本假定是满足的。
(11)该公司预计下一周签发新保单张,需要的加班时间时多少?
11.1问题求解
11.1.1输入
x0<-1000
y0<-beta_0+beta_1*x0
11.1.2输出
(12)给出的置信度为95%的精确预测区间和近似预测区间
12.1问题分析
可以计算得到 的分布为
所以枢轴量为
统计量为
得精确置信区间为
近似预测区间为
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)给出的置信度为95%的区间估计
13.1问题分析
根据 构造枢轴量(含)
统计量为
得置信水平95%得精确置信区间为
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输出
的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输出
回归方程如下
参数表如下:
参数 | 参数值 |
---|---|
10.28 | |
4.91 | |
0.92 | |
2.31 | |
[7.2,13.34] | |
[4.07,5.76] |
一元线性回归方差分析表如下:
方差来源 | 自由度 | 平方和 | 均方 | 值 | 值 |
---|---|---|---|---|---|
回归 | 1 | ||||
残差 | |||||
总和 |
回归方程的 的检验中,值为
数据的散点图,回归图,残差图如下
(2)结果分析
因为决定系数为0.92,具有比较强的线性相关性,且残差均在 内波动,而线性回归系数的检验通过,所以可以认为火灾发生地点与最近的消防站距离和火灾损失呈线性关系,符合模型