非参数模型采用GCV选择带宽的核估计核局部多项式估计

gcv选择局部多项式的带宽

K <- function(x) {
0.75 * (1 - x^2) * (abs(x) <= 1)
}
gcvhLL <- function(h,x,y) {
n <- length(x)
S <- matrix(NA,nrow = n,ncol = n)
GCVh <- 0
for (ii in 1:n) {
u=x-x[ii]
Kh=K(u/h)
s1=t(Kh)%%u
s2=t(Kh)%%(u^2)
if(s2 < 0) s2=s2+0.0002
s=sum(Kh*(s2-s1u))
if(s<0.0001) s=s+0.001
weight=as.vector(Kh(s2-s1u)/s)
GCVh=GCVh+(y[ii]-weight%%y)^2
S[ii, ]=weight
}
I=diag(n)
S1=I-S
trace=0
for (j in 1:n) trace=trace + S1[j, j]
GCVh=(GCVh/n)/(trace/n)^2
return(GCVh)
}

选择带宽的函数

gcvbandLL <- function(x,y,c=1.1) {
x <- as.vector(x)
y <- as.matrix(y)
n <- dim(y)[1]
tmp <- range(x)
hmin <- 2*(tmp[2] - tmp[1])/n
hmax <- 0.5*(tmp[2]-tmp[1])
h <- hmin
IUP <- 0
H <- vector()
H[0] <- 0
H[1] <- h
mt <- vector()
mt[1] <- gcvhLL(h,x,y)
i <- 1
while (h <= hmax && IUP < 3){
i <- i + 1
h <- h*c
H[i] <- h
mt[i] <- gcvhLL(h, x, y)
if (!is.na(mt[i]) && !is.na(mt[i - 1]) && as.numeric(mt[i]) >= as.numeric(mt[i - 1])) IUP <- 0 else IUP <- IUP + 1 # 如果当前的一个大于前面的
}
h <- H[mt == min(mt)][1]
return(list(h = h, i = i, mt = mt))
}

gcv选择核估计的带宽

gcvhNW <- function(h, x, y) {
n <- dim(x)[1]
S <- matrix(NA, nrow = n, ncol = n)
GCVh <- 0
for (i in 1:n) {
u <- x - x[i]
Kh <- K(u / h)
s <- sum(Kh)
if (s < 0.001) s <- s + 0.001
weight <- as.matrix(Kh / s)
GCVh <- GCVh + sum((y[i] - t(weight) %*% y)^2)
S[i, ] <- weight
}
I <- diag(n)
S1 <- I - S
trace <- 0
for (j in 1:n) trace <- trace + S1[j, j]
GCVh <- (GCVh / n) / (trace / n)^2
return(GCVh)
}

选择带宽的函数

gcvbandNW <- function(x, y, c = 1.1) {
x <- as.matrix(x)
y <- as.matrix(y)
n <- length(y)
tmp <- range(x)
hmin <- 6 * (tmp[2] - tmp[1]) / n
hmax <- (tmp[2] - tmp[1])
h <- hmin
IUP <- 0
H <- vector()
H[0] <- 0
H[1] <- h
mt <- vector()
mt[1] <- gcvhNW(h, x, y)
i <- 1

while (h <= hmax && IUP < 3) {
    i <- i + 1
    h <- h * c
    H[i] <- h # 存放带宽的取值
    mt[i] <- gcvhNW(h, x, y) # 存放求得的值
    if (!is.na(mt[i]) && !is.na(mt[i - 1]) && as.numeric(mt[i]) >= as.numeric(mt[i - 1])) IUP <- 0 else IUP <- IUP + 1 # 如果当前的一个大于前面的
}
h <- H[mt == min(mt)][1]
return(list(h = h, i = i, mt = mt))

}

估计过程

K=function(x)
0.75*(1-x^2)(abs(x)<=1) #ev核
m=function(x)
exp(-x^2)
k=20 #格子点个数
grid=rep(0,k)
a=-1;b=1
#grid[i]=a+(i-1)(b-a)/(k-1)
for(i in 1:k)
grid[i]=a+(i-1)*(b-a)/(k-1)
m1=rep(0,k)
m2=rep(0,k)
n=100
M=500
rmse1=rep(0,M)
rmse2=rep(0,M)
#产生x,e
for(j in 1:M){
x=runif(n,-1,1)
e=rep(0,n)
for(i in 1:n)
e[i]=rnorm(1,0,abs(x[i])) #放的标准差

#产生y
y=m(x)+e
h2=gcvbandLL(x,y)$hh1=gcvbandNW(x,y)$h
est1=rep(0,k)
est2=rep(0,k)
for(i in 1:k){
x1=x[abs(x-grid[i])<=h1]
y1=y[abs(x-grid[i])<=h1]
kvec=K((x1-grid[i])/h1)
length(y1)=length(kvec)
s=sum(kvec)
if(s==0) s = s + 1e-04
est1[i]=sum(kvec*y1)/s
}
m1=m1+est1
rmse1[j]=sqrt(mean((est1-m(grid))^2))
est2=lle(x,y,h2,grid)[,1]
m2=m2+est2
rmse2[j]=sqrt(mean((est2-m(grid))^2))
}
#######
est.NW=m1/M
est.LL=m2/M
op=par(mfrow=c(1,2))
plot(grid,m(grid),type=“l”,xlab=“x”,ylab=“m(x)”,xlim=c(-1,1),ylim=c(0.4,1))
lines(grid,est.NW,lty=2,col=2)
lines(grid,est.LL,lty=3,col=3,lwd=2)
legend(x=‘topright’,legend=c(“实际”,“NW”,“LL”),lty=c(1,2,3),lwd=2,col=c(1,2,3))
title(sub=“a”)
boxplot(data.frame(rmse1,rmse2),xlab=" ",ylab=“rmse”,range=3.5,names=c(expression(rmse1),expression(rmse2)))
title=(sub=“b”)