共计 3714 个字符,预计需要花费 10 分钟才能阅读完成。
咱们晓得参数的置信区间的计算,这些都遵从肯定的散布 (t 散布、正态分布),因而在规范误前乘以相应的 t 分值或 Z 分值。但如果咱们找不到适合的散布时,就无奈计算置信区间了吗?侥幸的是,有一种办法简直能够用于计算各种参数的置信区间,这就是 Bootstrap 法。
本文应用 BOOTSTRAP 来取得预测的置信区间。咱们将在线性回归根底上探讨。
> reg=lm(dist~speed,data=cars)
> points(x,predict(reg,newdata= data.frame(speed=x)))
这是一个单点预测。当咱们想给预测一个置信区间时,预测的置信区间取决于参数估计误差。
预测置信区间
让咱们从预测的置信区间开始
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")
蓝色值是通过在咱们的观测数据库中从新取样取得的可能预测值。值得注意的是,在残差正态性假如下(回归线的斜率和常数估计值),置信区间(90%)如下所示:
predict(reg,interval ="confidence",
在这里,咱们能够比拟 500 个生成数据集上的值散布,并将教训分位数与正态假如下的分位数进行比拟,
> hist(Yx,proba=TRUE
> boxplot(Yx,horizontal=TRUE
> polygon(c( x ,rev(x I]))))
能够看出,教训分位数与正态假如下的分位数是能够比拟的。
> quantile(Yx,c(.05,.95))
5% 95%
58.63689 70.31281
+ level=.9,newdata=data.frame(speed=x))
fit lwr upr
1 65.00149 59.65934 70.34364
感兴趣变量的可能值
当初让咱们看看另一种类型的置信区间,对于感兴趣变量的可能值。这一次,除了提取新样本和计算预测外,咱们还将在每次绘制时增加噪声,以取得可能的值。
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ base=cars[indice,]
+ erreur=residuals(reg)
+ predict(reg,newdata=data.frame(speed=x))+E
在这里,咱们能够(首先以图形形式)比拟通过从新取样取得的值和在正态假如下取得的值,
> hist(Yx,proba=TRUE)
> boxplot(Yx) abline(v=U[2:3)
> polygon(c(D$x[I,rev(D$x[I])
数值上给出了以下比拟
> quantile(Yx,c(.05,.95))
5% 95%
44.43468 96.01357
U=predict(reg,interval ="prediction"
fit lwr upr
1 67.63136 45.16967 90.09305
这一次,右侧有轻微的不对称。显然,咱们不能假如高斯残差,因为有更大的正值,而不是负值。思考到数据的性质,这是有意义的(制动间隔不能是正数)。
而后开始探讨在供给中应用回归模型。为了取得具备独立性,有人认为必须应用增量付款的数据,而不是累计付款。
能够创立一个数据库,解释变量是行和列。
> base=data.frame(
+ y
> head(base,12)
y ai bj
1 3209 2000 0
2 3367 2001 0
3 3871 2002 0
4 4239 2003 0
5 4929 2004 0
6 5217 2005 0
7 1163 2000 1
8 1292 2001 1
9 1474 2002 1
10 1678 2003 1
11 1865 2004 1
12 NA 2005 1
而后,咱们能够从基于对数增量付款数据的回归模型开始,该模型基于对数正态模型
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.9471 0.1101 72.188 6.35e-15 ***
as.factor(ai)2001 0.1604 0.1109 1.447 0.17849
as.factor(ai)2002 0.2718 0.1208 2.250 0.04819 *
as.factor(ai)2003 0.5904 0.1342 4.399 0.00134 **
as.factor(ai)2004 0.5535 0.1562 3.543 0.00533 **
as.factor(ai)2005 0.6126 0.2070 2.959 0.01431 *
as.factor(bj)1 -0.9674 0.1109 -8.726 5.46e-06 ***
as.factor(bj)2 -4.2329 0.1208 -35.038 8.50e-12 ***
as.factor(bj)3 -5.0571 0.1342 -37.684 4.13e-12 ***
as.factor(bj)4 -5.9031 0.1562 -37.783 4.02e-12 ***
as.factor(bj)5 -4.9026 0.2070 -23.685 4.08e-10 ***
---
Signif. codes: 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residual standard error: 0.1753 on 10 degrees of freedom
(15 observations deleted due to missingness)
Multiple R-squared: 0.9975, Adjusted R-squared: 0.9949
F-statistic: 391.7 on 10 and 10 DF, p-value: 1.338e-11
>
exp(predict(reg1,
+ newdata=base)+summary(reg1)$sigma^2/2)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 2871.2 1091.3 41.7 18.3 7.8 21.3
[2,] 3370.8 1281.2 48.9 21.5 9.2 25.0
[3,] 3768.0 1432.1 54.7 24.0 10.3 28.0
[4,] 5181.5 1969.4 75.2 33.0 14.2 38.5
[5,] 4994.1 1898.1 72.5 31.8 13.6 37.1
[6,] 5297.8 2013.6 76.9 33.7 14.5 39.3
> sum(py[is.na(y)])
[1] 2481.857
这与链式梯度法的后果略有不同,但依然具备可比性。咱们也能够尝试泊松回归(用对数链接)
glm(y~
+ as.factor(ai)+
+ as.factor(bj),data=base,
+ family=poisson)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.05697 0.01551 519.426 < 2e-16 ***
as.factor(ai)2001 0.06440 0.02090 3.081 0.00206 **
as.factor(ai)2002 0.20242 0.02025 9.995 < 2e-16 ***
as.factor(ai)2003 0.31175 0.01980 15.744 < 2e-16 ***
as.factor(ai)2004 0.44407 0.01933 22.971 < 2e-16 ***
as.factor(ai)2005 0.50271 0.02079 24.179 < 2e-16 ***
as.factor(bj)1 -0.96513 0.01359 -70.994 < 2e-16 ***
as.factor(bj)2 -4.14853 0.06613 -62.729 < 2e-16 ***
as.factor(bj)3 -5.10499 0.12632 -40.413 < 2e-16 ***
as.factor(bj)4 -5.94962 0.24279 -24.505 < 2e-16 ***
as.factor(bj)5 -5.01244 0.21877 -22.912 < 2e-16 ***
---
Signif. codes: 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 46695.269 on 20 degrees of freedom
Residual deviance: 30.214 on 10 degrees of freedom
(15 observations deleted due to missingness)
AIC: 209.52
Number of Fisher Scoring iterations: 4
> predict(reg2,
newdata=base,type="response")
> sum(py2[is.na(y)])
[1] 2426.985
预测后果与链式梯度法失去的估计值吻合。克劳斯·施密特(Klaus Schmidt)和安吉拉·温什(Angela Wünsche)于 1998 年在链式梯度法、边际和最大似然预计中建设了与最小偏差办法的分割。