[Rscript]商场销售预测-使用gbm和xgboost（包括调参）

其他 2018-07-17 20:55:26 阅读次数: 0 
#加载包
library(ggplot2)
library(grid)
library(gridExtra)
library(caret)
library(randomForest)
library(xgboost)
library(Hmisc)
library(dplyr)
library(rpart)
library(rpart.plot)
library(rattle)
getwd()
setwd('D:\\Users\\yang\\Desktop\\商场销售预测')
#加载数据
train.data <- read.csv("train.csv",stringsAsFactors = F)
test.data <- read.csv("test.csv",stringsAsFactors = F) 
test.data$Item_Outlet_Sales <- NA             #预测变量(测试数据的变为NA)，该产品在对应商品的销售额
#合并数据集，填补变量需要
all_data <- rbind(train.data,test.data) 
#描述统计量
describe(all_data)
summary(all_data)
#计算商品重量均值，默认移除缺失值
mmm1 <- aggregate(Item_Weight~Item_Identifier,data = all_data,FUN = mean)
#将均值赋值给对应缺失的重量值
for (i in which(is.na(all_data$Item_Weight))) {
  all_data$Item_Weight[i] <- mmm1$Item_Weight[mmm1$Item_Identifier == all_data$Item_Identifier[i]]
} 
#检查缺失值，处理后没有缺失值 
sum(is.na(all_data$Item_Weight))
#检查商店大小与商店类型的关系
table(all_data$Outlet_Size)
mmm2 <- aggregate(Item_Outlet_Sales~Outlet_Identifier+Outlet_Type+Outlet_Size,data = all_data,FUN = mean)
mmm2
#构建一个决策树来补插缺失值
#非缺失部分作为训练集，缺失部分作为测试集
fit <- rpart(factor(Outlet_Size)~Outlet_Type,data = all_data[all_data$Outlet_Size!="",],method = "class")
pred <- predict(fit,all_data[all_data$Outlet_Size=="",],type = "class") 
all_data$Outlet_Size[all_data$Outlet_Size==""] <- as.vector(pred)
#可以看到每个商店对应的大小都被补充完整了 
table(all_data$Outlet_Size,all_data$Outlet_Identifier)
#计算每个商品的销售量，向上取整 
all_data$Item_Sales_Vol <- ceiling(all_data$Item_Outlet_Sales/all_data$Item_MRP)
#低脂商品写法有误差
table(all_data$Item_Fat_Content) 
#将商品正确归类
all_data$Item_Fat_Content <- as.character(all_data$Item_Fat_Content)
all_data$Item_Fat_Content[all_data$Item_Fat_Content %in% c("LF","low fat")] <- "Low Fat"
all_data$Item_Fat_Content[all_data$Item_Fat_Content %in% c("reg")] <- "Regular" 
table(all_data$Item_Fat_Content)
#商品的类别可以进一步归类
str(all_data$Item_Type)
#通过查看商品ID，可以看到ID开头是商品的大类，DR饮品,FD食物,NC非消耗品，DR和FD为日常消耗品。ID后是A-Z
all_data$Item_Attribute <-  factor(substr(all_data$Item_Identifier,1,2))
table(all_data$Item_Attribute)
#由于食物才有低脂的情况，非消耗品则需要不同区分，生成一个新类别NF
all_data$Item_Fat_Content[all_data$Item_Attribute == "NC"] <- "NF"
table(all_data$Item_Fat_Content)
#新建变量Item_Attribute：DR饮品,FD食物,NC非消耗品三个类别
#由于存在非消耗品，对Item_Fat_Content新增一个因子，Non-Food 非食品

#去掉'OUT010‘'OUT019'后按照Item_Identifier各类均值填补百分比为0的异常值
sum(all_data$Item_Visibility==0)                        #879个缺失值
all_data1<-all_data[all_data$Outlet_Identifier!='OUT010' & all_data$Outlet_Identifier!='OUT019',]
sum(all_data1$Item_Visibility==0)                 #去掉'OUT010‘'OUT019'后776个缺失值
mmm3 <- aggregate(Item_Visibility~Item_Identifier,data = all_data1,FUN = mean)
for (i in which(all_data1$Item_Visibility==0)) {
  all_data1$Item_Visibility[i] <- mmm3$Item_Visibility[mmm3$Item_Identifier == all_data1$Item_Identifier[i]]
} 
sum(all_data1$Item_Visibility==0)
####填补'OUT010‘'OUT019'对应的缺失值
all_data2<-all_data[all_data$Outlet_Identifier=='OUT010' | all_data$Outlet_Identifier=='OUT019',]
sum(all_data2$Item_Visibility==0)                 #'OUT010‘'OUT019'103个缺失值
mmm4 <- aggregate(Item_Visibility~Outlet_Identifier,data = all_data2,FUN = mean)
for (i in which(all_data2$Item_Visibility==0)) {
  all_data2$Item_Visibility[i] <- mmm4$Item_Visibility[mmm4$Outlet_Identifier == all_data2$Outlet_Identifier[i]]
} 
sum(all_data2$Item_Visibility==0)
##分开填补后再合并数据集
all_data <- rbind(all_data1,all_data2)
str(all_data)
sum(all_data$Item_Visibility==0)                  #填补后不存在缺失值
aggregate(Item_Visibility~Outlet_Identifier,data = all_data,FUN = sum)        #查看填补后各个商店商品占比百分比之和

#Outlet_Establishment_Year与假设的客户流量相关，成立越久知名度越高，这是2013年收集的数据
all_data$Outlet_Established_period <- 2013-all_data$Outlet_Establishment_Year    #生成变量成立年数
table(all_data$Outlet_Established_period)           #共有9个不同的值，10个商店，所以变量采用因子型
#将所有分类型变量转为因子型
cols <- c("Item_Fat_Content","Item_Type","Outlet_Location_Type","Outlet_Type","Outlet_Established_period","Item_Attribute","Outlet_Identifier") 
for (i in cols) {
  all_data[,i] <- factor(all_data[,i])
}

#分割数据集
train <- all_data[!is.na(all_data$Item_Outlet_Sales),]
test <- all_data[is.na(all_data$Item_Outlet_Sales),]
write.csv(train,'train_fill.csv')
write.csv(test,'test_fill.csv')
write.csv(all_data,'all_data_fill.csv')

ggplot(train, aes(x =  Item_Attribute, y = Item_Sales_Vol))+ 
  geom_boxplot(aes(fill=Outlet_Type)) +
  theme(axis.title.x =element_text(size=14), axis.title.y=element_text(size=14),legend.text = element_text(size=12),legend.title = element_text(size=14),axis.text = element_text(size=12))
ggplot(train, aes(x =  Item_MRP, y = Item_Sales_Vol),colour=train$Outlet_Type)+
  geom_point()+geom_smooth( method = "loess",se=F)+
  scale_colour_brewer(palette = 'Set1')+facet_wrap(~Outlet_Type)
ggplot(train, aes(x =  Item_Type, y = Item_Sales_Vol),colour=train$Outlet_Type)+
  geom_boxplot(aes(fill=Outlet_Type))+geom_smooth( method = "loess",se=F)+
  scale_colour_brewer(palette = 'Set1')+facet_wrap(~Outlet_Type)

###读取数据 
setwd('C:/Users/lenovo/Desktop/Big_mart_sales_III/填补后的数据')
train <- read.csv('train_fill.csv')
test <- read.csv('test_fill.csv')
#设置种子，创建训练数据的训练子集和测试子集
set.seed(1234)
ind <- createDataPartition(train$Item_Sales_Vol,p=.7,list = FALSE)
train_val <- train[ind,]
test_val <- train[-ind,]
#需要建模的变量，尝试采用不同的变量情况
#创建模型评估RMSE函数 
model.rmse <- function(pred,act){ sqrt(sum((act - pred)^2)/length(act)) }
###回归分析
myformula <- Item_Sales_Vol ~  Outlet_Type + Item_Visibility + Outlet_Location_Type + 
  Item_MRP + Item_Type + Outlet_Established_period + Outlet_Size 
lm<-lm(myformula,data = train_val)
step.lm<-step(lm)
pred <- predict(step.lm,test_val) 
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
#对训练集预测
pred.test <- predict(step.lm,test)
#销量*单价得到销售额预测结果
pred_result <- data.frame(Item_Identifier = test.data$Item_Identifier, Outlet_Identifier = test.data$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred_result
ggplot()+geom_point(aes(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales))+ggtitle('线性回归')

####回归分析test_val测试结果
pred_result1 <- data.frame(Item_Identifier = test_val$Item_Identifier,Item_Outlet_Sales=pred*test_val$Item_MRP) 
pred_result1
x<-NULL
x$regression<-pred_result1$Item_Outlet_Sales
quant<-NULL
quant$regression<-pred
  
####
myformula <- Item_Sales_Vol ~  Outlet_Type + Item_Visibility + Outlet_Location_Type + 
  Item_MRP + Item_Type + Outlet_Established_period + Outlet_Size       #1090.73
#决策树建模
library(zoo)
library(party)
fit_ctree<-ctree(myformula,data = train_val)
print(fit_ctree)
plot(fit_ctree)
plot(fit_ctree,type='simple')
pred <- predict(fit_ctree,test_val) 
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
#画好看点的决策树的图，为何用fancyRpartPlot?
fit.tr <- rpart(myformula,data = train_val,method = "anova")
summary(fit.tr)
rpart.plot(fit.tr)
library(rattle)
fancyRpartPlot(fit.tr)

##7个和11个预测的销售额结果误差相同
#对训练集预测
pred.test <- predict(fit.tr,test)
#销量*单价得到销售额预测结果
pred_result <- data.frame(Item_Identifier = test.data$Item_Identifier, Outlet_Identifier = test.data$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred_result
plot(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
ggplot()+geom_point(aes(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales))+ggtitle('决策树')

####决策树test_val测试结果
pred_result2 <- data.frame(Item_Identifier = test_val$Item_Identifier,Item_Outlet_Sales=pred*test_val$Item_MRP) 
pred_result2
x$tr<-pred_result2$Item_Outlet_Sales
quant$tr<-pred


##随机森林
set.seed(2345) 
fit.rf <- randomForest(myformula,data = train_val,ntree=500)
summary(fit.rf)
pred <- predict(fit.rf,test_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
# 变量重要性
library(pROC)
varImpPlot(fit.rf, main = "变量重要性")
# AUC ROC
pre<-predict(fit.rf)
true.value <- train_val$Item_Sales_Vol
modelroc <- multiclass.roc(true.value,pre)  #AUC=0.7626
#创建用于上传评分的测试结果rf.csv 
pred.test <- predict(fit.rf,test)
pred_result <- data.frame(Item_Identifier = test.data$Item_Identifier, Outlet_Identifier = test.data$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred_result

ggplot()+geom_point(aes(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales))+ggtitle('随机森林')

####随机森林test_val测试结果
pred_result3 <- data.frame(Item_Identifier = test_val$Item_Identifier,Item_Outlet_Sales=pred*test_val$Item_MRP) 
pred_result3
x$rf<-pred_result3$Item_Outlet_Sales
quant$rf<-pred

####
#控制器,5折交叉验证
Ctrl <- trainControl(method="repeatedcv",number=5,repeats=5)
set.seed(3456)
fit.gbm <- train(myformula,data = train_val,trControl=Ctrl,method="gbm",verbose=FALSE) 
barplot()
summary(fit.gbm) 
pred <- predict(fit.gbm,test_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)    #1087.254

pred.test <- predict(fit.gbm,test)
pred.result <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred.result

ggplot()+geom_point(aes(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales))+ggtitle('GBM')

####交叉验证test_val测试结果
pred_result4 <- data.frame(Item_Identifier = test_val$Item_Identifier,Item_Outlet_Sales=pred*test_val$Item_MRP) 
pred_result4
x$gbm<-pred_result4$Item_Outlet_Sales
quant$gbm<-pred

####xgboost
#构建稀疏矩阵
mymatrix <- function(train){
  matrix_num <- train[,c("Item_Visibility","Item_MRP")]
  matrix_num <- cbind(matrix_num,
                      model.matrix(~Outlet_Type-1,train),
                      model.matrix(~Outlet_Location_Type-1,train),
                      model.matrix(~Outlet_Size-1,train),
                      model.matrix(~Item_Type-1,train),
                      model.matrix(~Outlet_Established_period-1,train))
  return(data.matrix(matrix_num))
}
#获取每个数据集的稀疏矩阵
xgb.train_val <- mymatrix(train_val)
xgb.test_val <- mymatrix(test_val) 
xgb.test <- mymatrix(test)
#生成用于xgboots模型的DMatrix
#注意，预测变量集和响应变量是要分开的
dtrain_val <- xgb.DMatrix(data =xgb.train_val,label=train_val$Item_Sales_Vol)
dtest_val <- xgb.DMatrix(data = xgb.test_val,label=test_val$Item_Sales_Vol)
dtest_sub <- xgb.DMatrix(data = xgb.test)

library(xgboost)
#初步xgboost建模
model <- xgboost(data = dtrain_val,nround = 5)
summary(model) 
pred <- predict(model,dtest_val) 
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales) 
xgb.importance(colnames(xgb.train_val),model)
pred.test <- predict(model,dtest_sub)
pred.result <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred.result
##增加循环次数为10，但为了防止过度拟合，将默认的树模型最大深度6调整为5
model_tuned <- xgboost(data = dtrain_val,nround = 10,max.depth = 5)
summary(model_tuned)
pred <- predict(model_tuned,dtest_val) 
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
xgb.importance(colnames(xgb.train_val),model_tuned)
pred.test <- predict(model_tuned,dtest_sub)
pred.result <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred.result

#在模型调整循环11次，最大深度为4
for (i in 13:19) {
  for (j in 3:4) {
    model_tuned <- xgboost(data = dtrain_val,nround = i,max.depth = j)
    pred <- predict(model_tuned,dtest_val) 
    print(i)
    print(j)
    print(model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales))
  }
}
##i=8,j=3时1077.234

model_tuned <- xgboost(data = dtrain_val,nround = 8,max.depth = 3)
summary(model_tuned)
pred <- predict(model_tuned,dtest_val) 
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
xgb.importance(colnames(xgb.train_val),model_tuned)
pred.test <- predict(model_tuned,dtest_sub)
pred.result <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP) 
pred.result

ggplot()+geom_point(aes(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales))+ggtitle('xgboost')

####交叉验证test_val测试结果
pred_result5 <- data.frame(Item_Identifier = test_val$Item_Identifier,Item_Outlet_Sales=pred*test_val$Item_MRP) 
pred_result5
x$xgboost<-pred_result5$Item_Outlet_Sales
quant$xgboost<-pred

#查看模型中各个变量的权重
xgb.importance(colnames(xgb.train_val),model_tuned)

####################支持向量回归
library(DMwR)
library(nnet)
train_val
norm.data<-train_val[,c(5,6,7,10,11,12,16,14)] #数据标准化#c(5,6,7,10,11,12,16,14)
library(e1071)
model.svm<-svm(Item_Sales_Vol~.,norm.data)
preds<-predict(model.svm,norm.data)
plot(preds,norm.data$Item_Sales_Vol)
library(rminer)
set.seed(1234)
svr<-fit(Item_Sales_Vol~.,norm.data,model='svm')
#利用模型进行预测
norm.preds<-predict(svr,norm.data)
#绘制预测值与真实值之间的散点图
plot(norm.preds,norm.data$Item_Sales_Vol)


norm.preds<-predict(svr,test_val[,c(5,6,7,10,11,12,16,14)])
ggplot()+geom_point(aes(norm.preds*test_val$Item_MRP,test_val$Item_Outlet_Sales))+ggtitle('svr')

#计算相对误差
model.rmse(norm.preds*test_val$Item_MRP,test_val$Item_Outlet_Sales)
plot(norm.preds*test_val$Item_MRP,test_val$Item_Outlet_Sales)
####交叉验证test_val测试结果
pred_result6 <- data.frame(Item_Identifier = test_val$Item_Identifier,Item_Outlet_Sales=norm.preds*test_val$Item_MRP) 
pred_result6
x$svr<-pred_result6$Item_Outlet_Sales
quant$svr<-norm.preds


x$act<-test_val$Item_Outlet_Sales
quant$act<-test_val$Item_Sales_Vol
write.table(quant,'quant.txt')
z<-read.table('z.txt')
quant<-read.table('quant.txt')
qq<-matrix(0,nrow=17892,ncol=3)
qq[,2]<-c(rep(quant$act,7))
qq[,1]<-c(rep(1,2556),rep(2,2556),rep(3,2556),rep(4,2556),rep(5,2556),rep(6,2556),rep(7,2556))
qq[,3]<-c(quant$regression,quant$tr,quant$rf,quant$gbm,quant$xgboost,quant$svr,quant$act)
colnames(qq)<-c('group','act','pre')
qq<-as.data.frame(qq)
qq$group<-as.factor(qq$group)

ggplot(qq,aes(x=act,y=pre,colour=group))+
  geom_point()+geom_smooth( method = "loess",se=F)+
  scale_colour_brewer(palette = 'Dark2')

ggplot(xxx,aes(x=act,y=pre,colour=group))+
  geom_point()+geom_smooth( method = "loess",se=F)+
  scale_colour_brewer(palette = 'Set1')

################画出每组分别的图像
ggplot(xxx,aes(x=act,y=pre)) +
  geom_point() +
  geom_smooth(span = 0.8) +
  facet_wrap(~group)

猜你喜欢

转载自blog.csdn.net/tomocat/article/details/80956398
今日推荐
周排行
Leetcode简单题61~80
解决zookeeper磁盘IO高的问题
多线程相关方法详解
Maven-setting.xml文件详解
Maven 项目的 classpath 理解
渊亭科技大数据笔试题
配置JVM内存分配
计算机网络个人学习笔记 （三）网络层 ：第三部分 连载
js中两个等号(==)和三个等号(===)的区别
用C程序自动打开电脑上的程序
每日归档
更多
2024-09-18(0)
2024-09-17(0)
2024-09-16(0)
2024-09-15(0)
2024-09-14(0)
2024-09-13(0)
2024-09-12(0)
2024-09-11(0)
2024-09-10(0)
2024-09-09(0)

代码天地 © 2018 codetd.com

境外服务器   最新电视剧2022