关于算法:R语言中的多项式回归局部回归核平滑和平滑样条回归模型

46次阅读

共计 3306 个字符,预计需要花费 9 分钟才能阅读完成。

在规范线性模型中,咱们假如 ​。当线性假如无奈满足时,能够思考应用其余办法。

  • 多项式回归

扩大可能是假如某些多项式函数,

同样,在规范线性模型办法(应用 GLM 的条件正态分布)中,参数 ​ 能够应用最小二乘法取得,其中 ​ 在 ​。

即便此多项式模型不是真正的多项式模型,也可能依然是一个很好的近似值 ​。实际上,依据 Stone-Weierstrass 定理,如果 ​ 在某个区间上是间断的,则有一个对立的近似值 ​,通过多项式函数。

仅作阐明,请思考以下数据集

 db = data.frame(x=xr,y=yr)
plot(db)

与规范回归线

reg = lm(y ~ x,data=db)
abline(reg,col="red")

思考一些多项式回归。如果多项式函数的次数足够大,则能够取得任何一种模型,

reg=lm(y~poly(x,5),data=db)

然而,如果次数太大,那么会取得太多的“稳定”,

reg=lm(y~poly(x,25),data=db)

并且估计值可能不牢靠:如果咱们更改一个点,则可能会产生(部分)更改

 yrm=yr;yrm[31]=yr[31]-2 
lines(xr,predict(regm),col="red") 

  • 部分回归

实际上,如果咱们的趣味是部分有一个很好的近似值  ​,为什么不应用部分回归?

应用加权回归能够很容易地做到这一点,在最小二乘公式中,咱们思考

  • 在这里,我思考了线性模型,然而能够思考任何多项式模型。在这种状况下,优化问题是

​能够解决,因为

例如,如果咱们想在某个时候进行预测,思考 ​。应用此模型,咱们能够删除太远的观测值,

更个别的想法是思考一些核函数 ​ 给出权重函数,以及给出邻域长度的一些带宽(通常示意为 h),

这实际上就是所谓的 Nadaraya-Watson 函数预计器 ​。
在后面的案例中,咱们思考了对立核 
​,

然而应用这种权重函数具备很强的不连续性不是最好的抉择,尝试高斯核,

这能够应用

 w=dnorm((xr-x0))
reg=lm(y~1,data=db,weights=w) 

在咱们的数据集上,咱们能够绘制

 w=dnorm((xr-x0))
plot(db,cex=abs(w)*4)
lines(ul,vl0,col="red")
axis(3)
axis(2)
reg=lm(y~1,data=db,weights=w)
u=seq(0,10,by=.02)
v=predict(reg,newdata=data.frame(x=u))
lines(u,v,col="red",lwd=2) 

在这里,咱们须要在点 2 进行部分回归。上面的水平线是回归(点的大小与宽度成比例)。红色曲线是部分回归的演变

让咱们应用动画来可视化曲线。

然而因为某些起因,我无奈在 Linux 上轻松装置该软件包。咱们能够应用循环来生成一些图形

 name=paste("local-reg-",100+i,".png",sep="")
png(name,600,400)
 

for(i in 1:length(vx0)) graph (i)

而后,我应用

当然,能够思考部分线性模型,

 return(predict(reg,newdata=data.frame(x=x0)))}

甚至是二次(部分)回归,

 lm(y~poly(x,degree=2), weights=w) 

当然,咱们能够更改带宽

请留神,实际上,咱们必须抉择权重函数(所谓的核)​​。然而,有(简略)办法来抉择“最佳”带宽 h。穿插验证的想法是思考

 ​ 是应用部分回归取得的预测。

咱们能够尝试一些实在的数据。

library(XML)
 
data = readHTMLTable(html) 

整顿数据集,

 plot(data$no,data$mu,ylim=c(6,10))
segments(data$no,data$mu-1.96*data$se, 

咱们计算标准误差,反映不确定性。

 for(s in 1:8){reg=lm(mu~no,data=db, 
lines((s predict(reg)[1:12] 

所有节令都应该被认为是齐全独立的,这不是一个很好的假如。

 smooth(db$no,db$mu,kernel = "normal",band=5) 

咱们能够尝试查看带宽较大的曲线。

db$mu[95]=7
  
plot(data$no,data$mu

lines(NW,col="red")

样条平滑

接下来,探讨回归中的平滑办法。假如 ​,​ 是一些未知函数,但假设足够平滑。例如,假如 ​ 是间断的,​ 存在,并且是间断的,​ 存在并且也是间断的等等。如果 ​ 足够平滑,能够应用泰勒展开式。因而,对于 

也能够写成

第一局部只是一个多项式。

应用 黎曼积分,察看到

因而,

咱们有线性回归模型。一个天然的想法是思考回归 ​,对于 ​ 

给一些节点 ​。

plot(db)

如果咱们思考一个节点,并扩大阶数 1,

 B=bs(xr,knots=c(3),Boundary.knots=c(0,10),degre=1)
 
lines(xr[xr<=3],predict(reg)[xr<=3],col="red")
lines(xr[xr>=3],predict(reg)[xr>=3],col="blue")

能够将用该样条取得的预测与子集(虚线)上的回归进行比拟。

 lines(xr[xr<=3],predict(reg)[xr<=3


lm(yr~xr,subset=xr>=3) 

这是不同的,因为这里咱们有三个参数(对于两个子集的回归)。当要求间断模型时,失去了一个自由度。察看到能够等效地写

 lm(yr~bs(xr,knots=c(3),Boundary.knots=c(0,10) 

回归中呈现的函数如下

当初,如果咱们对这两个重量进行回归,咱们失去

 matplot(xr,B

abline(v=c(0,2,5,10),lty=2)

如果加一个节点,咱们失去

预测是

 lines(xr,predict(reg),col="red")

咱们能够抉择更多的节点

 lines(xr,predict(reg),col="red")

咱们能够失去一个置信区间

 polygon(c(xr,rev(xr)),c(P[,2],rev(P[,3]))

points(db) 

如果咱们放弃先前抉择的两个节点,但思考泰勒的 2 阶的开展,咱们失去

 matplot(xr,B,type="l")
abline(v=c(0,2,5,10),lty=2)

如果咱们思考常数和基于样条的第一局部,咱们失去

 B=cbind(1,B)
lines(xr,B[,1:k]%*%coefficients(reg)[1:k],col=k-1,lty=k-1)

如果咱们将常数项,第一项和第二项相加,则咱们失去的局部在第一个节点之前位于左侧,

k=3
lines(xr,B[,1:k]%*%coefficients(reg)[1:k] 

通过基于样条的矩阵中的三个项,咱们能够失去两个节点之间的局部,

 lines(xr,B[,1:k]%*%coefficients(reg)[1:k] 

最初,当咱们对它们求和时,这次是最初一个节点之后的右侧局部,

k=5 

这是咱们应用带有两个(固定)节点的二次样条回归失去的后果。能够像以前一样取得置信区间

 polygon(c(xr,rev(xr)),c(P[,2],rev(P[,3]))

points(db)
lines(xr,P[,1],col="red") 

应用函数 ​,能够确保点的连续性 ​。

再一次,应用线性样条函数,能够减少连续性束缚,

 lm(mu~bs(no,knots=c(12*(1:7)+.5),Boundary.knots=c(0,97),


lines(c(1:94,96),predict(reg),col="red")

然而咱们也能够思考二次样条,

 abline(v=12*(0:8)+.5,lty=2)

lm(mu~bs(no,knots=c(12*(1:7)+.5),Boundary.knots=c(0,97), 

正文完
 0