R软件做非参数回归（样条回归、多项式回归）（代码+数据）

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

    样条回归、局部多项式回归。
    代码全是根据算法自编函数，可以自行设置参数，尤其是窗宽的的选择更加自由。压缩包里有这部分讲义、自编代码、练习题。


贴出部分源代码

---

1. #------------------定义函数---------------------#<?xml:namespace prefix = "o" ns = "urn:schemas-microsoft-com:office:office" />
   
2. kern.epan <- function(u){
   
3. return(3/4*(1-u^2)*(abs(u)<=1));
   
4. }
   
5. 
   
6. LLS.fun <- function(x,X,Y,h)
   
7. {
   
8.   z   = 1/h*kern.epan((x-X)/h);     
   
9.   s1  = sum(z*(x-X));
   
10.   s2  = sum(z*(x-X)^2);
    
11.   w   = z*(s2-(x-X)*s1);
    
12. 
    
13.   return(sum(w*Y)/sum(w));  
    
14. }
    
15. 
    
16. CV1.fun=function(X,Y,h)
    
17. {
    
18.    y1=rep(0,length(Y));
    
19.    for(i in 1:length(Y)){
    
20.       y1[i]=LLS.fun(X[i],X[-i],Y[-i],h);
    
21.       }
    
22.    return(mean((Y-y1)^2));   
    
23. }
    
24. 
    
25. CI.boot.fun=function(resids,xpoints,X,yfit,B,h)
    
26. {
    
27. # bootstrap方法求置信区间
    
28.    result=matrix(0,length(xpoints),B);
    
29.    for(i in 1:length(xpoints)){
    
30.        for(b in 1:B){
    
31.            epsilon=sample(resids,replace=TRUE);
    
32.            newY=yfit+epsilon;
    
33.            #X2=X^2; X3=X^3; X4=X^4;
    
34.            #fit1 <- lm(Y~X+X2+X3+X4);
    
35.            #coefE=c(fit1$coeff);
    
36.            #sigmaE=sqrt(var(fit1$residuals));
    
37. 
    
38.            #CK=1.719 # Epanechnikov kernel used.
    
39.            #temp=cbind(2,3*2*X,4*3*X^2)%*%as.vector(coefE[-(1:2)]);
    
40.            #den=sum(temp^2);
    
41.            #h.ROT=CK*(sigmaE^2/den)^(1/(2*1+3));
    
42. 
    
43. 
    
44.            #newfit=LLS.fun(xpoints[i],X,newY,h=h.ROT);
    
45.            newfit=LLS.fun(xpoints[i],X,newY,h=h);
    
46.            
    
47.            result[i,b]=newfit;
    
48.            }
    
49.        }
    
50.       
    
51.    return(result);
    
52. 
    
53. }

复制代码

局部多项式法

---

1. #------------------局部多项式法---------------------#
   
2. #------------------Homework.one---------------------#
   
3. library(SemiPar);
   
4. data(age.income);
   
5. attach(age.income);
   
6. X=age;
   
7. Y=log.income;
   
8. plot(X,Y); #画原始数据的散点图
   
9. 
   
10. X2=X^2; X3=X^3; X4=X^4;
    
11. fit1<-lm(Y~X)
    
12. lines(X,fit1$fitted.values,lty=1,col=1,lwd=3)
    
13. fit2<-lm(Y~X+X2)
    
14. lines(X,fit2$fitted.values,lty=2,col=2,lwd=3)
    
15. fit3<-lm(Y~X+X2+X3)
    
16. lines(X,fit3$fitted.values,lty=3,col=3,lwd=3)
    
17. fit4<-lm(Y~X+X2+X3+X4)
    
18. lines(X,fit4$fitted.values,lty=4,col=4,lwd=3)
    
19. legend(20,15.1,c("线性","二次","三次","四次"),lty=c(1,2,3,4),col=c(1,2,3,4))
    
20. 
    
21. # Rule-of-Thumb Global Bandwidth Selector;
    
22. X2=X^2; X3=X^3; X4=X^4;
    
23. fit1 <- lm(Y~X+X2+X3+X4);
    
24. coefE=c(fit1$coeff);
    
25. resids=fit1$residuals;
    
26. sigmaE=sqrt(var(resids));
    
27. CK=1.719 # Epanechnikov kernel used.
    
28. temp=cbind(2,3*2*X,4*3*X^2)%*%as.vector(coefE[-(1:2)]);
    
29. den=sum(temp^2);
    
30. h.ROT=CK*(sigmaE^2/den)^(1/(2*1+3));
    
31. print(h.ROT)  #选择窗宽
    
32. 
    
33. #final fitting;
    
34. xfine=seq(0,70,length=50);
    
35. ypred1 <- rep(0,length(xfine));
    
36. for(i in 1:length(xfine)){
    
37.     ypred1[i] <- LLS.fun(xfine[i],X,Y,h=8);   
    
38. }
    
39. plot(X,Y);
    
40. lines(xfine,ypred1,lty=1,col=3);#拟合数据
    
41. 
    
42. 
    
43. #求g(x)的置信区间
    
44. yfit=rep(0,length(X));
    
45. for(i in 1:length(X)){
    
46. yfit[i] <- LLS.fun(X[i],X,Y,h=8);
    
47. }
    
48. 
    
49. temp=Y-yfit;
    
50. resids=temp-mean(temp);
    
51. result1=CI.boot.fun(resids=resids,xpoints=xfine,X=X,yfit=yfit,B=1000,h=8)
    
52. CI.upp=rep(0,length(xfine));
    
53. CI.low=rep(0,length(xfine));
    
54. for(i in 1:length(xfine)){
    
55.    temp=quantile(result1[i,],probs=c(0.025,0.975),na.rm = T);
    
56.    CI.low[i]=temp[1]; CI.upp[i]=temp[2];
    
57. }
    
58. plot(X,Y)
    
59. lines(xfine,ypred1,lty=1,col=2);                        
    
60. lines(xfine,CI.upp,lty=3,col=3);  
    
61. lines(xfine,CI.low,lty=4,col=4);
    
62. legend(20,15.1,c("Estmate","CI.upp","CI.low"),lty=c(1,3,4),col=c(2,3,4))

复制代码

资料下载：R软件做非参数回归拟合数据，有代码和数据

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏21 回帖

关键词：非参数回归 样条回归 非参数 多项式 r软件 非参数回归 样条回归 代码 非参数回归 样条回归 局部多项式回归 R

感谢楼主，感谢分享，发帖不易，用心良苦。

谢谢楼主分享，非常感谢！

学习了，谢谢！
但是你的代码有些看不太懂，比如B-spline回归那部分，有些函数没有太理解什么意思，请问有这方面的书籍推荐一下吗，特别是关于B-spline回归的？

xiexielouzhu

好帖，感谢分享

楼主是东财校友吗？也在做分位数回归，可以进一步交流下吗？

学习了，谢谢！

能不能分享一个啊，我论坛币不够了，1309731862@qq.com

一心一半888 发表于 2018-4-22 14:49
能不能分享一个啊，我论坛币不够了，1309731862@qq.com

我下载了发你邮箱