藥店銷量預測

本案例大綱：

1、xgboost理論介紹

2、R語言中xgboost相關函式的引數

3、案例背景

4、資料預處理

5、R語言的xgb模型實現程式碼

1、xgboost理論介紹

      這部分我直接把一些牛人寫的關於xgb的理論介紹引用過來了，大家可以直接看以下部落格連結資料，既有原理介紹又有程式碼的函式引數介紹：

       http://blog.csdn.net/a819825294/article/details/51206410

      http://blog.csdn.net/sb19931201/article/details/52557382

      http://blog.csdn.net/sb19931201/article/details/52577592

2、R語言中xgboost相關函式的引數

        R語言的XGBOOST包的引數包括三個方面的引數：常規引數、模型引數和任務引數。通用引數用於選擇哪一類分類器，是樹模型還是線性模型；模型引數取決於常規函式中選擇的模型型別；任務引數取決於學習的場景。

常規數：

模型引數：

（1）樹模型的引數

eta[default=0.3] 

學習率，一般設定小一些。
range: [0,1] 
gamma [default=0] 
後剪枝時，用於控制是否剪枝，值越大，演算法越保守。
range: [0,∞] 
max_depth [default=6] 
樹的最大深度
範圍: [1,∞] 
min_child_weight [default=1] 
這個引數預設是 1，是每個葉子裡面 h 的和至少是多少，對正負樣本不均衡時的 0-1 分類而言，假設 h 在 0.01 附近，min_child_weight 為 1 意味著葉子節點中最少需要包含 100 個樣本。這個引數非常影響結果，控制葉子節點中二階導的和的最小值，該引數值越小，越容易 overfitting。
range: [0,∞] 
max_delta_step [default=0] 
這個引數在更新步驟中起作用，如果取0表示沒有約束，如果取正值則使得更新步驟更加保守。可以防止做太大的更新步子，使更新更加平緩。
range: [0,∞] 
subsample [default=1] 
樣本隨機取樣，較低的值使得演算法更加保守，防止過擬合，但是太小的值也會造成欠擬合。
range: (0,1] 
colsample_bytree [default=1] 
列取樣，對每棵樹的生成用的特徵進行列取樣.一般設定為： 0.5-1
range: (0,1] 
lambda [default=1] 
權重L2正則化
alpha [default=0] 
權重L1正則化

（2）線性模型引數

lambda[default=0] 
權重L2正則化
alpha [default=0] 
權重L1正則化
lambda_bias 
L2 regularization term on bias, default 0(no L1 reg on bias because it is notimportant)

偏導L2正則化引數，預設為0（沒有偏導L1正則化引數）

任務引數：

objective [default=reg:linear ]   定義最小化損失函式型別，常用引數如下：
“reg:linear” –linear regression 
“reg:logistic” –logistic regression 
“binary:logistic” –logistic regression for binary classification, outputprobability 
“binary:logitraw” –logistic regression for binary classification, output scorebefore logistic transformation 
“count:poisson” –poisson regression for count data, output mean of poissondistribution 
max_delta_step is set to 0.7 by default in poisson regression (used tosafeguard optimization) 
“multi:softmax” –set XGBoost to do multiclass classification using the softmaxobjective, you also need to set num_class(number of classes) 
“multi:softprob” –same as softmax, but output a vector of ndata * nclass, whichcan be further reshaped to ndata, nclass matrix. The result contains predictedprobability of each data point belonging to each class. 
“rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwiseloss 
base_score [ default=0.5 ] 
the initial prediction score of all instances, global bias 所有例項的初始預測評分, 全域性偏差
eval_metric [ default according to objective ]   評估指標選擇
evaluation metrics for validation data, a default metric will be assignedaccording to objective( rmse for regression, and error for classification, meanaverage precision for ranking ) 
User can add multiple evaluation metrics, for Python user, remember to pass the metrics in as list of parameters pairsinstead of map, so that latter ‘eval_metric’ won’t override previous one 
The choices are listed below: 評估指標可選列表如下：
“rmse”: root mean square error 
“logloss”: negative log-likelihood 
“error”: Binary classification error rate. It is calculated as #(wrongcases)/#(all cases). For the predictions, the evaluation will regard theinstances with prediction value larger than 0.5 as positive instances, and theothers as negative instances. 
“merror”: Multiclass classification error rate. It is calculated as #(wrongcases)/#(all cases). 
“mlogloss”: Multiclass logloss 
“auc”: Area under the curve for ranking evaluation. 
“ndcg”:Normalized Discounted Cumulative Gain 
“map”:Mean average precision 
“[email protected]”,”[email protected]”: n can be assigned as an integer to cut off the top positionsin the lists for evaluation. 
“ndcg-“,”map-“,”[email protected]“,”[email protected]“: In XGBoost, NDCG and MAP will evaluate thescore of a list without any positive samples as 1. By adding “-” in theevaluation metric XGBoost will evaluate these score as 0 to be consistent undersome conditions. training repeatively 
seed [ default=0 ] 隨機種子
random number seed.

3、案例背景

        Rossmann在7個歐洲國家擁有3,000家藥店。目前，羅斯曼店經理的任務是提前六週預測其日銷量。商店銷售受到諸多因素的影響，包括促銷，競爭，學校和國家假日，季節性和地點。成千上萬的個人經理根據其獨特的情況預測銷售量，結果的準確性可能會有很大的變化。

      Kaggle所提供的資料的欄位如下：

        由於本案例主要講解xgboost模型，所以對於資料預處理和特徵工程都做得比較少。只做了兩方面的處理，第一，Kaggle官網把商店的一些屬性資料與訓練集、測試集分開放，在不同檔案，所以要把store資料集與train、test資料集按列合併；第二，將資料按照xgb要求的格式進行轉換，R語言的xgboost包的xgb.Matrix是轉換格式的包。

5、程式碼實現

資料下載地址：https://www.kaggle.com/c/rossmann-store-sales/data

library(readr)
library(xgboost)
library(lubridate)
train<-read.csv('D:/R語言kaggle案例實戰/Kaggle第五節課/train.csv')
test<-read.csv('D:/R語言kaggle案例實戰/Kaggle第五節課/test.csv')
store<-read.csv('D:/R語言kaggle案例實戰/Kaggle第五節課/store.csv')#store是對店鋪屬性的補充
train<-merge(train,store)#將兩個資料集按列合併
test<-merge(test,store)#將兩個資料集按列合併
train$Date<-as.POSIXct(train$Date)#將日期字元變成時間格式
test$Date<-as.POSIXct(test$Date)#將日期字元變成時間格式
train[is.na(train)]<-0#將空值置為零
test[is.na(test)]<-0
train<-train[which(train$Open=='1'),]#選擇開門的且銷售額不為0的樣本
train<-train[which(train$Sales!='0'),]
train$month<-month(train$Date)#提取月份
train$year<-year(train$Date)#提取年份
train$day<-day(train$Date)#提取日
train<-train[,-c(3,8)]#刪除日期列和缺失值較多的列
test<-test[,-c(4,7)]#刪除日期列和缺失值較多的列
feature.names<-names(train)[c(1,2,5:19)]#這一步主要使測試集和訓練集的結構一致。
for(f in feature.names){
  if(class(train[[f]])=="character"){
    levels<-unique(c(train[[f]],test[[f]]))
    train[[f]]<-as.integer(factor(train[[f]],levels = levels))
    test[[f]]<-as.integer(factor(test[[f]],levels = levels))
  }
}
tra<-train[,feature.names]
RMPSE<-function(preds,dtrain){ #定義一個評價函式，Kaggle官方給的評價函式作為xgboost中的評價函式。
  labels<-getinfo(dtrain,"label")
  elab<-exp(as.numeric(labels))-1
  epreds<-exp(as.numeric(preds))-1
  err<-sqrt(mean((epreds/elab-1)^2))
  return(list(metric="RMPSE",value=err))
}
h<-sample(nrow(train),10000)#進行10000次抽樣
dval<-xgb.DMatrix(data=data.matrix(tra[h,]),label=log(train$Sales+1)[h])#用於以下構建watchlist 
dtrain<-xgb.DMatrix(data=data.matrix(tra[-h,]),label=log(train$Sales+1)[-h])#構建xgb特定的矩陣形式
watchlist<-list(val=dval,train=dtrain)#構建模型引數的watchlist,watchlist是用於監聽每次模型執行時的模型效能情況。
param<-list(objective="reg:linear",
            booster="gbtree",
            eta=0.02,
            max_depth=12,
            subsample=0.9,
            colsample_bytree=0.7,
            num_parallel_tree=2,
            alpha=0.0001,
            lambda=1)
clf<-xgb.train(  params=param,
                 data=dtrain,
                 nrounds = 3000,
                 verbose = 0,
                 early.stop.round=100,
                 watchlist = watchlist,
                 maximize = FALSE,
                 feval = RMPSE
  
)
ptest<- predict(clf,test,outputmargin=TRUE)