求高手用R或者其他软件解决一元个非参数回归方面的问题

以下是图，这次把图也传上来

对39楼的朋友致以崇高的谢意！
深感论坛上高手如云和人品的伟大，向你们致敬！

我曾经做过一元非参数核密度估计，是用R做的，下面附有整个源代码，你只要输入你的数据就行，一起分享

个人觉得你真要把实证做好，就要懂那些计量理论，
如果你用非参数方法，至少要掌握非参数各种估计
方法的统计推断，在建立模型时必须要符合经济学
假设以及计量理论假设，否则，所得到的结果与
真实会相差很多，所以既然要用到计量模型，我们
就必须考虑这些问题。
37# zhangtao

刚发现我的程序中核函数编错了， 编成了三次权重函数，楼主可以把它改成正太核函数。

45楼的朋友，要不你把程序修改后重新上传比较好！说句说话，这个正态核函数我不会改。
非常感谢！
非常赞同44楼朋友的意见！

45# 可爱的河马
#标准正太核函数

kernel.normal <- function(u){ y <- ifelse(abs(u)>=1,0, (1 - abs(u)^3)^3)}
根据你的说法，把上面一行删除，程序运行出来的结果与以前一样，
结果如下：                      这是为什么？希望指教，谢谢！

> y<-c(1510.2,1700.6,2026.6,2577.4,3496.2,4283,4838.9,5160.3,5425.1,5854,6280,6859.6,7702.8,8472.2,9421.6,10493,11759.5,13785.8,15780.8,17174.7)
> x<-c(1278.89,1453.81,1671.73,2110.81,2851.34,3537.57,3919.47,4185.64,4331.61,4615.91,4998,5309.01,6029.88,6510.94,7182.1,7942.88,8696.55,9997.47,11242.85,12264.55)
>
>
> #最优窗宽h,
> library("KernSmooth")
> h<-dpill(x, y)
> h
[1] 4474.243
>
> #和回归函数
> kernel.smooth1 <- function(X, Y, kernel, h, plot.it=T){
+     x <- X
+     n<-length(X)
+     fx <- numeric(n)
+     for(j in 1:n){
+
+     fx[j] <- sum(kernel((x[j]-X)/h)*Y) /sum(kernel((x[j]-X)/h))
+     }
+     if(plot.it){
+     plot(X, Y, type="p")
+     lines(x, fx)
+     }
+     fx=fx
+   }
>
> f<-kernel.smooth1(X=x, Y=y, kernel=kernel.normal, h=h)
>
> #非参数回归的拟合值
> y.hat<-f
> y.hat
[1]  3075.116  3211.858  3373.403  3671.418  4144.037  4627.864  4941.489
[8]  5185.778  5327.592  5620.397  6039.715  6379.198  7084.285  7586.481
[15]  8434.089  9680.494 11034.642 13168.157 14707.591 15632.828
>
> #均方误差
> mse<-mean((y.hat-y)^2)
> mse
[1] 794223.4

我把核函数改成正太核， 另外均方误差是用回归值减y，平方后再平均，不知对不对？

#在R中运行前检查一下R是否安装了"KernSmooth"包，
#我输入了两个数据，后面的自己输入。
y<-c(1278.89, 1453.81, ......)
x<-c(1510.2, ......)

#标准正太核函数
kernel.normal <- function(u) (2*pi)^(-1/2)*exp(-u^2/2)
#最优窗宽h,
library("KernSmooth")
h<-dpill(x, y)
h

#和回归函数
kernel.smooth1 <- function(X, Y, kernel, h, plot.it=T){
    x <- X
    n<-length(X)
    fx <- numeric(n)
    for(j in 1:n){

    fx[j] <- sum(kernel((x[j]-X)/h)*Y) /sum(kernel((x[j]-X)/h))
    }
    if(plot.it){
    plot(X, Y, type="p")
    lines(x, fx)
    }
    fx=fx
  }

f<-kernel.smooth1(X=x, Y=y, kernel=kernel.normal, h=h)

#非参数回归的拟合值
y.hat<-f
y.hat

#均方误差
mse<-mean((y.hat-y)^2)
mse

y<-c(1510.2,1700.6,2026.6,2577.4,3496.2,4283,4838.9,5160.3,5425.1,5854,6280,6859.6,7702.8,8472.2,9421.6,10493,11759.5,13785.8,15780.8,17174.7)
> x<-c(1278.89,1453.81,1671.73,2110.81,2851.34,3537.57,3919.47,4185.64,4331.61,4615.91,4998,5309.01,6029.88,6510.94,7182.1,7942.88,8696.55,9997.47,11242.85,12264.55)
> #标准正太核函数
> kernel.normal <- function(u) (2*pi)^(-1/2)*exp(-u^2/2)
> #最优窗宽h,
> library("KernSmooth")
> h<-dpill(x, y)
> h
[1] 4474.243
>
> #和回归函数
> kernel.smooth1 <- function(X, Y, kernel, h, plot.it=T){
+     x <- X
+     n<-length(X)
+     fx <- numeric(n)
+     for(j in 1:n){
+
+     fx[j] <- sum(kernel((x[j]-X)/h)*Y) /sum(kernel((x[j]-X)/h))
+     }
+     if(plot.it){
+     plot(X, Y, type="p")
+     lines(x, fx)
+     }
+     fx=fx
+   }
>
> f<-kernel.smooth1(X=x, Y=y, kernel=kernel.normal, h=h)
>
> #非参数回归的拟合值
> y.hat<-f
> y.hat
[1]  5033.772  5090.737  5163.163  5314.155  5584.923  5855.052  6013.834
[8]  6128.199  6192.233  6319.659  6496.637  6645.604  7008.019  7263.194
[15]  7636.672  8083.398  8547.992  9390.117 10222.199 10903.259
>
> #均方误差
> mse<-mean((y.hat-y)^2)
> mse
[1] 8104366
MSE太大，从图上也可以看出，估计程序有问题，希望可爱的河马朋友再看看，谢谢！

程序应该没有问题，估计是窗宽的问题，我也不会估计窗宽，直接用的dpill函数生成的窗宽，你问问其他人，谁会估计窗宽。 或者自己给定一个窗宽值，例如：h<-2000    , 多试几个窗宽值。