部分线性非参数的模拟基于NW估计和局部多项式估计

非参数部分线性模型的模拟 NW估计、局部多项式估计

setwd("e:/统计模型") #下面三个函数代码的位置，在另外一个文章中有此函数的代码
source("gcvbandLL.R")
source("gcvbandNW.R")
source("lle.R")

K=function(x)
	0.75*(1-x^2)*(abs(x)<=1)
g=function(u)
	50*((u-0.5)^2)
k=20
grid=rep(0,k)
a=0;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)
grid
g(grid)
m1=rep(0,k)
m2=rep(0,k)
n=100
M=5
rmse1=rep(0,M)
rmse2=rep(0,M)

beta_NW=matrix(NA,M,2)#两种方法beta的估计
beta_LL=matrix(NA,M,2)
guNW=matrix(NA,M,k)
guLL=matrix(NA,M,k)

# NW估计
#产生x,u,e
for(j in 1:M){
	x1=rnorm(n,0,1)
	x2=rnorm(n,0,1)
	e=rnorm(n,0,0.6)
	x=cbind(x1,x2)
	#x=as.matrix(x)
	real_beta=c(0.8,0.6)
	u=runif(n,0,1)
	#u=as.matrix(u)


#产生y
	y=g(u)+x%*%real_beta+e
	#y=as.matrix(y)
	h_NW1=gcvbandNW(u,y)$h
	h_NW2=gcvbandNW(u,x)$h
	g1NW=c()
	for(i1 in 1:n){
		u1=u[abs(u-u[i1])<=h_NW1] #认识
		y1=y[abs(u-u[i1])<=h_NW1]
		kvec=K((u1-u[i1])/h_NW1)
		length(y1)=length(kvec)
		s=sum(kvec)
		if(s==0) s = s + 0.0001
		g1NW[i1]=sum(kvec*y1)/s
	}
	g2NW=matrix(NA,n,2)
	for(i2 in 1:n){
		u1=u[abs(u-u[i2])<=h_NW2]#u1不是na 认识
		x1=x[abs(u-u[i2])<=h_NW2,]
		kvec=K((u1-u[i2])/h_NW2)
		length(y1)=length(kvec)
		s=sum(kvec)
		if(s==0) s = s + 0.0001
		g2NW[i2,]=sum(kvec*x1)/s #放在一行
	}

######先把beta估计出来
	beta_NW[j,]=solve(t(x-g2NW)%*%(x-g2NW))%*%t(x-g2NW)%*%(y-g1NW)
	g1NW=c()
	for(g1 in 1:k){
	
		u1=u[abs(u-grid[g1])<=h_NW1]#u1是na 不认识u了
		y1=y[abs(u-grid[g1])<=h_NW1]
		kvec=K((u1-grid[g1])/h_NW1)
		length(y1)=length(kvec)
		s=sum(kvec)
		if(s==0) s = s + 0.001
		g1NW[g1]=sum(kvec*y1)/s
	}
	g2NW=matrix(NA,k,2)
	for(g2 in 1:k){
		u1=u[abs(u-grid[g2])<=h_NW2]#u1s是na 
		x1=x[abs(u-grid[g2])<=h_NW2]
		kvec=K((u1-grid[g2])/h_NW2)
		length(y1)=length(kvec)
		s=sum(kvec)
		if(s==0) s=s+0.001
		g2NW[g2,]=sum(kvec*x1)/s #放在一行
	}

	guNW[j,]=g1NW-g2NW%*%beta_NW[j,]
	rmse1[j]=sqrt(mean((guNW[j,]-g(grid))^2))
}

# 局部线性
M=5 #模拟次数
kk=20 #格子点
#产生x,u,e
for(j1 in 1:M){
	x1=rnorm(n,0,1)
	x2=rnorm(n,0,1)
	e=rnorm(n,0,0.6)
	x=cbind(x1,x2)

	real_beta=c(0.8,0.6)
	u=runif(n,0,1)



#产生y
	y=g(u)+x%*%real_beta+e
	h_LL1=gcvbandLL(u,y)$h
	h_LL2=gcvbandLL(u,x)$h
	g1LL=c()
	g1LL=lle(u,y,h_LL1,u)[,1]
	g2LL=matrix(NA,n,2)
	g2LL=lle(u,x,h_LL2,u)
######先把beta估计出来
	beta_LL[j1,]=solve(t(x-g2LL)%*%(x-g2LL))%*%t(x-g2LL)%*%(y-g1LL)
##########格子点估计gu
	g1LL=c()
	g1LL=lle(u1,y1,h_LL1,grid)[,1]
	g2LL=matrix(NA,kk,2)
	g2LL=lle(u,x,h_LL2,grid)
	guLL[j1,]=g1LL-g2LL%*%beta_LL[j1,]
	rmse2[j1]=sqrt(mean((guLL[j1,]-g(grid))^2))
}

## 求偏差，画图
colMeans(beta_NW)-real_beta # bias
c(sd(beta_NW[,1]),sd(beta_NW[,2])) # sd
colMeans((beta_NW-real_beta)^2) # MSE

colMeans(beta_LL)-real_beta # bias
c(sd(beta_LL[,1]),sd(beta_LL[,2])) # sd
colMeans((beta_LL-real_beta)^2) # MSE


#####################
est.NW=colMeans(guNW)
est.LL=colMeans(guLL)

op=par(mfrow=c(1,2))
plot(grid,g(grid),type="l",xlab="u",ylab="g(u)",xlim=c(0,1),ylim=c(0,15))
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")