R语言数据分析系列之八，r语言数据分析
R语言数据分析系列之八，r语言数据分析
R语言数据分析系列之八
—— by comaple.zhang
再谈多项式回归，本节再次提及多项式回归分析，理解过拟合现象，并深入cross-validation(交叉验证),regularization(正则化)框架，来避免产生过拟合现象，从更加深入的角度探讨理论基础以及基于R如何将理想照进现实。
本节知识点，以及数据集生成
1， ggplot2进行绘图;
2， 为了拟合更复杂的数据数据集采用sin函数加上服从正太分布的随机白噪声数据;
3， poly函数进行高次拟合,transform进行数据集变换;
4， 理解过拟合，以及如何利用正则化框架避免过拟合.
测试数据生成：
library(ggplot2) x <- seq(0,1,by=0.01) y <- sin(2*pi*x)+rnorm(length(x),0,0.1) df <- data.frame(x,y) ggplot(df,aes(x=x,y=y),main="sin(x)")+ geom_point()
高次函数拟合
先看一下直线拟合的效果：
fit1 <- lm(y ~ x,df) df <- transform(df,yy = predict(fit1)) ggplot(df,aes(x=x,y=y)) + geom_point()+geom_smooth(aes(x=x,y=yy),data=df)
#利用函数poly进行高次拟合，该函数可以将自变量自动进行高次变化，degree参数控制着最高次项的次数如下：
po <- poly(x,degree=3) fit3 <- lm(y ~ poly(x,degree=3),df) df <- transform(df,yy3 = predict(fit3)) ggplot(df,aes(x=x,y=y)) + geom_point()+geom_line(aes(x=x,y=yy3),data=df,col='blue')
我们继续提高最高次的次数：
fit26 <- lm(y ~ poly(x,degree=26),df) df <- transform(df,yy26 =predict(fit26)) ggplot(df,aes(x=x,y=y)) + geom_point()+geom_line(aes(x=x,y=yy26),data=df,col='blue')
从图中我们可以明显看到，并不是次数越高越好，那么如何确定合适得次数呢，一般认为次数越高，我们认为模型越复杂，复杂的模型能很好的拟合原始数据但泛化能力不一定最好，我倾向于找到一个简单，并且拥有良好的泛化能力的模型。我们可以首先从数字角度去分析，我们调用summary函数看结果可以得知这些参数中只有 1，3，5次方项通过了t验证，其他项目均未通过。那么我们如何保证模型的可用性呢，下面就介绍方法cross-validation和正则化框架。
模型验证，cross-validation
为了能够更好的选择适合的模型，我们引模型评价方法RMSE(均方根误差)
均方根误差可以用来比较两个模型的好坏。
我们用R来实现这个评价方法：
rmse <- function(y,ry)
{
return(sqrt(sum((y-ry)^ 2))/length(ry))
}
cross-validation的思想是将一份数据集，分成两份，一份用来训练模型叫训练集df$train，一份用来测试叫测试集df$test,下边的函数将数据集切分为两部分：
split <- function(df,rate)
{
n <- length(df[,1])
index <- sample(1:n,round(rate * n))
train <- df[index,]
test <- df[-index,]
df <- list(train=train,test=test,data=df)
return(df)
}
下面我们来看一下在不同的次数变换下，rmse的分布对比情况：
performance_Gen <- function(df,n){
performance <- data.frame()
for(index in 1:n){
fit <- lm(y ~ poly(x,degree=index),data = df$train)
performance <- rbind(performance,data.frame(degree =index,type='train',rmse=rmse(df$train['y'],predict(fit))))
performance <- rbind(performance,data.frame(degree = index,type='test',rmse=rmse(df$test['y'],predict(fit,newdata=df$test))))
}
return(performance)
}
df_split <- split(df,0.5)
performance<- performance_Gen(df_split,10)
ggplot(performance,aes(x=degree,y=rmse,linetype=type))+geom_point()+geom_line()
从图中可以看出当degree=3时模型对测试集已经表现的很好，当degree增加时，rmse变得越来越大，这说明模型的误差越来越大，而模型对训练集的拟合却变现的非常好。所以这种方案是选择degree=3时的模型。那么有没有更好的方法呢，正则化框架为我们提供了另外一种避免过拟合的方法下面一起来看一下。
正则化框架
我们在前面说过多项式回归的损失函数为，请参见上一节（http://blog.csdn.net/comaple/article/details/44959695）：
那么我现在引入模型复杂度的概念来改进我们得模型损失函数，这个复杂度也是越小越好，即如果复杂度越大说明我们的模型是不好的，如果复杂度越小同时模型对训练集拟合的比较好，那说明我们的模型是好的，以此来避免过拟合和增加模型的泛化能力。于是我们的正则化框架改造如下：
这里的lambda其实是正则化框架的惩罚因子，lambda越大证明我们对模型复杂度越关心，lambda越小证明我们对模型复杂度不关心。f(w) 是一个关于参数w的函数。
1， L1范式，如果我们对w参数取绝对值累加，即为L1范式
2， L2范式，如果我们对w参数取平方和，那即为L2范式
3， Lp范式，如果我们对w参数取p次方累加和，即为Lp范式
当然范式的选择是根据经验和你的模型的作用来最终确定的，函数的形式也可以更加复杂，以便于进行特殊作用的处理。
这里以上面已经有数据集进行绘图说明
l1 <- abs(sort(coef(fit3))) df_l1 <-data.frame(x=sort(coef(fit3)),y=l1) ggplot(df_l1,aes(x=x,y=y))+geom_line()
l2 <- (coef(fit3)^2)
df_l2 <-data.frame(x=sort(coef(fit3)),y=l2)
ggplot(df_l2,aes(x=x,y=y))+geom_smooth()+ggtitle('L2norm')
L1范式可以进行参数选择，使得某些不重要的参数为零，这样就可以简化我们得模型，使得模型有更好的泛化能力
但是这里的lambda的值如何选取呢，不同得lambda代表着我们对模型复杂度的关心程度。这里又要用来cross-validation来验证那个lambda适合我们。
在R语言里glmnet包已经实现了正则化框架，我们这里安装并加载他。并调用glmnet函数，该函数只要调用一次就可以拿到所有可能的lambda值和对应的模型，那么我们可以上面的方法，把rmse的结果评测和lambda的关系通绘制出来，看看效果：
install.packages(glmnet)
library(glmnet)
getperform <- function(df){
fit<- with(df$train,glmnet(poly(x,degree=10),y))
lambdas <- fit$lambda
performance <- data.frame()
for( lambda in lambdas){
performance <- rbind(performance,data.frame(
lambda=lambda,
rmse=rmse(df$test['y'],with(df$test,predict(fit,poly(x,degree=10),s=lambda)))))
}
return(performance)
}
performance <- getperform(df_split)
ggplot(performance,aes(x=lambda,y=rmse))+geom_point()+geom_line()+ ggtitle('Lambda vs RMSE')
从图中可以看出，在lambda 0.06附近时rmse取最小值，故该数据集可以取lambda=0.06时对应的模型。
评论暂时关闭