<jupyter_start><jupyter_text>DataApplication Lab baseline 比赛: M5 时序预测本 NoteBook 是 Kaggle M5 时序预测比赛的 Baseline，主要使用机器学习的建模方式进行时序的建模和预测。- 特征工程来参考：https://www.kaggle.com/kneroma/-m5-first-public-notebook-under-0-50- groupKfold参考：https://www.kaggle.com/ragnar123/simple-lgbm-groupkfold-cv- **BaseLine步骤**：                            1. 数据分析 EDA                        2. 特征工程                        3. 模型调参                        4. stacking>**提示：**Code 和 Markdown 区域可通过 **Shift + Enter** 快捷键运行。此外，Markdown可以通过双击进入编辑模式。- 缺陷： 1. groupKfold参考kernel的线上分数是0.53，但是现在我们single model部分提交的线上结果分数达到4，要回头认真比对groupKfold 2. 在调参过程中也要更改folds=groupKfold。--- --- 1、数据分析为了节约时间，我们直接对我们后面建模有用的结果进行分析，关于数据的详细分析可以参考 [EDA](https://www.kaggle.com/headsortails/back-to-predict-the-future-interactive-m5-eda)- 查看 sales 数据前几行- 查看 sales 数据聚合结果趋势- 查看 sales 数据标签分布<jupyter_code>import numpy as np
from scipy import sparse
import pandas as pd
import lightgbm as lgb
import re
import string
import time
import seaborn as sns
import itertools
import sys

from sklearn import preprocessing, pipeline, metrics, model_selection
from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer
import gc
from datetime import datetime, timedelta
from sklearn.model_selection import GroupKFold
from sklearn import metrics

import matplotlib.pyplot as plt

pd.set_option('display.max_colwidth',100)

%matplotlib inline 
from google.colab import drive
drive.mount('/content/drive')
!cp -r /content/drive/My\ Drive/colab/kaggle_M5/m5-forecasting-accuracy/ /content/
sale_data = pd.read_csv('/content/m5-forecasting-accuracy/sales_train_validation.csv')
print(sale_data.shape)
sale_data.head(5)
day_data = sale_data[[f'd_{day}' for day in range(1,1914)]]
total_sum = np.sum(day_data,axis=0).values
plt.plot(total_sum)
plt.hist(day_data[day_data<100].values.reshape(-1),bins=100);<jupyter_output>/usr/local/lib/python3.6/dist-packages/numpy/lib/histograms.py:839: RuntimeWarning: invalid value encountered in greater_equal
  keep = (tmp_a >= first_edge)
/usr/local/lib/python3.6/dist-packages/numpy/lib/histograms.py:840: RuntimeWarning: invalid value encountered in less_equal
  keep &= (tmp_a <= last_edge)<jupyter_text>--- 2、特征工程选定机器学习的建模方案，核心思想是对时间序列抽取窗口特征。抽取窗口特征：- 前7天- 前28天- 前7天均值- 前28天均值关联其他维度信息- 日期- 价格<jupyter_code># helper functions to reduce memory
def reduce_mem_usage(df, verbose=True):
    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']
    start_mem = df.memory_usage().sum() / 1024**2    
    for col in df.columns:
        col_type = df[col].dtypes
        if col_type in numerics:
            c_min = df[col].min()
            c_max = df[col].max()
            if str(col_type)[:3] == 'int':
                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
                    df[col] = df[col].astype(np.int8)
                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
                    df[col] = df[col].astype(np.int16)
                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
                    df[col] = df[col].astype(np.int32)
                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
                    df[col] = df[col].astype(np.int64)  
            else:
                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
                    df[col] = df[col].astype(np.float16)
                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
                    df[col] = df[col].astype(np.float32)
                else:
                    df[col] = df[col].astype(np.float64)    
    end_mem = df.memory_usage().sum() / 1024**2
    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))
    return df
def create_train_data(train_start=750,test_start=1800,is_train=True):
    """
    traing_start: 取多长的数据进行训练, 从第750个到1931个
    test_start: 取多长的数据进行预测
    """
    # 基本参数
    PRICE_DTYPES = {"store_id": "category", "item_id": "category", "wm_yr_wk": "int16","sell_price":"float32" }
    CAL_DTYPES={"event_name_1": "category", "event_name_2": "category", "event_type_1": "category", 
            "event_type_2": "category", "weekday": "category", 'wm_yr_wk': 'int16', "wday": "int16",
            "month": "int16", "year": "int16", "snap_CA": "float32", 'snap_TX': 'float32', 'snap_WI': 'float32' }

    start_day = train_start if is_train else test_start
    numcols = [f"d_{day}" for day in range(start_day,1914)]
    catcols = ['id', 'item_id', 'dept_id','store_id', 'cat_id', 'state_id']
    SALE_DTYPES = {numcol:"float32" for numcol in numcols} 
    SALE_DTYPES.update({col: "category" for col in catcols if col != "id"})

    # 加载price数据
    price_data = pd.read_csv('/content/m5-forecasting-accuracy/sell_prices.csv',dtype=PRICE_DTYPES)
    print(f"price data shape is {price_data.shape}")
    # 加载cal数据
    cal_data = pd.read_csv('/content/m5-forecasting-accuracy/calendar.csv',dtype=CAL_DTYPES)
    print(f"calender data shape is {cal_data.shape}")
    # 加载sale数据
    sale_data = pd.read_csv('/content/m5-forecasting-accuracy/sales_train_validation.csv',dtype=SALE_DTYPES,usecols=catcols+numcols)
    print(f"sale data shape is {sale_data.shape}")


    # 类别标签转换
    for col, col_dtype in PRICE_DTYPES.items():
        if col_dtype == "category":
            price_data[col] = price_data[col].cat.codes.astype("int16")
            price_data[col] -= price_data[col].min()

    cal_data["date"] = pd.to_datetime(cal_data["date"])
    for col, col_dtype in CAL_DTYPES.items():
        if col_dtype == "category":
            cal_data[col] = cal_data[col].cat.codes.astype("int16")
            cal_data[col] -= cal_data[col].min()


    for col in catcols:
        if col != "id":
            sale_data[col] = sale_data[col].cat.codes.astype("int16")
            sale_data[col] -= sale_data[col].min()

    # 注意提交格式里有一部分为空
    if not is_train:
        for day in range(1913+1, 1913+ 2*28 +1):
            sale_data[f"d_{day}"] = np.nan

    sale_data = pd.melt(sale_data,
            id_vars = catcols,
            value_vars = [col for col in sale_data.columns if col.startswith("d_")],
            var_name = "d",
            value_name = "sales")
    print(f'sale data shape is {sale_data.shape} after melt')
    sale_data = sale_data.merge(cal_data, on= "d", copy = False)
    print(f'sale data shape is {sale_data.shape} after merging with cal_data')
    sale_data = sale_data.merge(price_data, on = ["store_id", "item_id", "wm_yr_wk"], copy = False)
    print(f'sale data shape is {sale_data.shape} after merging with price_data')
    return sale_data
sale_train_data = create_train_data(train_start=1500,is_train=True)
sale_test_data = create_train_data(is_train=False)
sale_train_data.tail()
sale_test_data.head()
sale_test_data.tail()
sale_data = pd.concat([sale_train_data,sale_test_data])

print ("Full Sale Data set created.")
print(sale_data.shape)
sale_data.columns
sale_data.head()
sale_data.tail()
# 可以在这里加入更多的特征抽取方法
# 获取7天前的数据，28天前的数据
# lags = [7, 28], 因为验证集28天现在sale为空，我们对验证集不能取lag=7
lags = [28]
lag_cols = [f"lag_{lag}" for lag in lags]

# 如果是测试集只需要计算一天的特征，减少计算量
# 注意训练集和测试集特征生成要一致

for lag, lag_col in zip(lags, lag_cols):
    sale_data[lag_col] = sale_data[["id","sales"]].groupby("id")["sales"].shift(lag)

# 将获取7天前的数据，28天前的数据做移动平均
# wins = [7, 28]
# use the rolling windows from 'simple GroupKFold CV'
wins = [7, 15, 30, 90]


for win in wins :
    for lag,lag_col in zip(lags, lag_cols):
        sale_data[f"rmean_{lag}_{win}"] = sale_data[["id", lag_col]].groupby("id")[lag_col].transform(lambda x : x.rolling(win).mean()) 
        sale_data[f'rstd_{lag}_{win}'] = sale_data[['id', lag_col]].groupby('id')[lag_col].transform(lambda x: x.rolling(win).std())

sale_data = reduce_mem_usage(sale_data)

# 处理时间特征
# 有的时间特征没有，通过datetime的方法自动生成
date_features = {
        "wday": "weekday",
        "week": "weekofyear",
        "month": "month",
        "quarter": "quarter",
        "year": "year",
        "mday": "day",
        "dayofweek": "dayofweek",
        "dayofyear": "dayofyear"
    }

for date_feat_name, date_feat_func in date_features.items():
    if date_feat_name in sale_data.columns:
        sale_data[date_feat_name] = sale_data[date_feat_name].astype("int16")
    else:
        sale_data[date_feat_name] = getattr(sale_data["date"].dt, date_feat_func).astype("int16")
sale_data.shape
sale_data[(sale_data['date'] >= '2016-01-25') & (sale_data['date'] <= '2016-05-22')]
sale_data.columns
# 清洗数据，选择需要训练的数据
# sale_data.dropna(inplace=True)
cat_feats = ['item_id', 'dept_id','store_id', 'cat_id', 'state_id'] + ["event_name_1", "event_name_2", "event_type_1", "event_type_2"]
useless_cols = ["id", "date", "sales","d", "wm_yr_wk"]
print('-'*50)
print('\n')
print(f'Training model with {len(sale_data.columns)} features...')
X_train = sale_data[(sale_data['date'] >= '2016-01-25') & (sale_data['date'] <= '2016-04-24')]
y_train = X_train["sales"]
X_test = sale_data[(sale_data['date'] >= '2016-04-25') & (sale_data['date'] <= '2016-05-22')]
# X_test.drop(['sales'], axis=1, inplace=True)

train_cols = X_train.columns[~X_train.columns.isin(useless_cols)]
X_train = X_train[train_cols]
X_test = X_test[train_cols]
X_train.head()
y_train.head()
print(X_train.shape)
print(y_train.shape)
X_train.reset_index(inplace=True, drop=True)
X_train.head()
y_train.reset_index(inplace=True, drop=True)
y_train.head()
X_test.head()
X_test.tail()
test = sale_data[(sale_data['date'] >= '2016-04-25') & (sale_data['date'] <= '2016-05-22')][['id', 'date', 'sales']]
test.head()
test.tail()
del(sale_train_data,sale_test_data)<jupyter_output><empty_output><jupyter_text>--- 3、模型调参  LightGBM Manual tuning<jupyter_code># define custom loss function
def custom_asymmetric_train(y_pred, y_true):
    y_true = y_true.get_label()
    residual = (y_true - y_pred).astype("float")
    grad = np.where(residual < 0, -2 * residual, -2 * residual * 1.15)
    hess = np.where(residual < 0, 2, 2 * 1.15)
    return grad, hess

# define custom evaluation metric
def custom_asymmetric_valid(y_pred, y_true):
    y_true = y_true.get_label()
    residual = (y_true - y_pred).astype("float")
    loss = np.where(residual < 0, (residual ** 2) , (residual ** 2) * 1.15) 
    return "custom_asymmetric_eval", np.mean(loss), False
import copy
greater_is_better = False

    # params = {
    #     "objective" : "tweedie",
    #     "metric" :"rmse",
    #     "force_row_wise" : True,
    #     "learning_rate" : 0.075,
    #     "sub_feature" : 0.8,
    #     "sub_row" : 0.75,
    #     "bagging_freq" : 1,
    #     "lambda_l2" : 0.1,
    #     "metric": ["rmse"],
    #     "nthread": 8,
    #     "tweedie_variance_power":1.1, # 老师发现1.2比较好。
    # 'verbosity': 1,
    # 'num_iterations' : 1500,
    # 'num_leaves': 128,
    # "min_data_in_leaf": 104,
    # }

# lgb_metric = 'rmse'
lgb_metric = "custom_asymmetric_eval"
default_lgb_params = {}
# default_lgb_params["objective"] = "tweedie"
# default_lgb_params["tweedie_variance_power"] = 1.1
default_lgb_params["learning_rate"] = 0.1 #0.05
default_lgb_params["metric"] = lgb_metric
default_lgb_params["bagging_freq"] = 1
# default_lgb_params['colsample_bytree'] = 0.85
# default_lgb_params['colsample_bynode'] = 0.85
# default_lgb_params['min_data_per_leaf'] = 25
default_lgb_params["seed"] = 1234

params_lgb_space = {}
params_lgb_space['feature_fraction'] = [0.1, 0.3, 0.5, 0.7, 0.9]
params_lgb_space['num_leaves'] = [3, 7, 15, 31, 63, 127]
params_lgb_space['max_depth'] = [3, 7, 10, 15, 31, -1]
params_lgb_space['min_gain_to_split'] = [0, 0.1, 0.3, 1, 1.5, 2, 3]
params_lgb_space['bagging_fraction'] = [0.2, 0.4, 0.6, 0.8, 1]
params_lgb_space['min_sum_hessian_in_leaf'] = [0, 0.0001, 0.001, 0.1, 1, 3, 10]
params_lgb_space['lambda_l2'] = [0, 0.01, 0.1, 0.5, 1, 10]
params_lgb_space['lambda_l1'] = [0, 0.01, 0.1, 0.5, 1, 10] 
kfolds = 5
kf = model_selection.GroupKFold(n_splits=kfolds)
group = X_train['week'].astype(str) + '_' + X_train['year'].astype(str)
kf_ids = list(kf.split(X_train, y_train, group))
kf_ids
from sklearn.model_selection import KFold,StratifiedKFold
kfolds = 5
random_state = 42

skf = StratifiedKFold(n_splits=kfolds, shuffle=True, random_state=random_state)
# folds = kf.split(X_train, y_train, group), 
#                 folds=skf.split(X_train, y_train),
best_lgb_params = copy.copy(default_lgb_params)

for p in params_lgb_space:
    print ("\n Tuning parameter %s in %s" % (p, params_lgb_space[p]))

    params = best_lgb_params
    scores = []    
    for v in params_lgb_space[p]:
        print ('\n    %s: %s' % (p, v), end="\n")
        params[p] = v
        lgb_cv = lgb.cv(params,
                lgb.Dataset(X_train,
                            label=y_train,
                            categorical_feature=cat_feats
                            ),
                num_boost_round=100000,
                fobj = custom_asymmetric_train,
                feval = custom_asymmetric_valid,
                nfold=kfolds,
                stratified=False,
                early_stopping_rounds=50,
                verbose_eval=500)
        # print(lgb_cv)
        if greater_is_better:
            best_lgb_score = min(lgb_cv['%s-mean' % (lgb_metric)])
        else:
            best_lgb_score = min(lgb_cv['%s-mean' % (lgb_metric)])
        best_lgb_iteration = len(lgb_cv['%s-mean' % (lgb_metric)])
        print (', best_score: %f, best_iteration: %d' % (best_lgb_score, best_lgb_iteration))
        scores.append([v, best_lgb_score])
    # best param value in the space
    best_param_value = sorted(scores, key=lambda x:x[1],reverse=greater_is_better)[0][0]
    best_param_score = sorted(scores, key=lambda x:x[1],reverse=greater_is_better)[0][1]
    best_lgb_params[p] = best_param_value
    print("Best %s is %s with a score of %f" %(p, best_param_value, best_param_score))

print('\n Best manually tuned parameters:', best_lgb_params)<jupyter_output>Tuning parameter feature_fraction in [0.1, 0.3, 0.5, 0.7, 0.9]

    feature_fraction: 0.1<jupyter_text>Best feature_fraction is 0.9 with a score of 4.608439Best num_leaves is 63 with a score of 4.593394Best max_depth is 31 with a score of 4.593394Best min_gain_to_split is 0.1 with a score of 4.609980Best bagging_fraction is 1 with a score of 4.611675Best min_sum_hessian_in_leaf is 0 with a score of 4.611248Best lambda_l2 is 10 with a score of 4.605351<jupyter_code>best_lgb_params = {'learning_rate': 0.05, 
                   'metric': 'custom_asymmetric_eval', 
                   'bagging_freq': 1, 
                   'colsample_bytree': 0.85, 
                   'colsample_bynode': 0.85, 
                   'min_data_per_leaf': 25, 
                   'seed': 1234, 
                   'lambda_l1': 0.01,
                   'feature_fraction': 0.9,
                   'num_leaves': 63,
                   'max_depth': 31,
                   'min_gain_to_split': 0.1,
                   'bagging_fraction': 1,
                   'min_sum_hessian_in_leaf': 0,
                   'lambda_l2': 10
                   }
# groupkfold cv
# Best manually tuned parameters: {'learning_rate': 0.05, 'metric': 'custom_asymmetric_eval', 
#'bagging_freq': 1, 'colsample_bytree': 0.85, 'colsample_bynode': 0.85, 'min_data_per_leaf': 25,
# 'seed': 1234, 'lambda_l1': 0.01, 
#'feature_fraction': 0.9, 'num_leaves': 63, 'max_depth': 31, 'min_gain_to_split': 0.1, 
#'bagging_fraction': 1, 'min_sum_hessian_in_leaf': 0, 'lambda_l2': 10}
print ('\n Best manually tuned parameters:', best_lgb_params)<jupyter_output>Best manually tuned parameters: {'learning_rate': 0.05, 'metric': 'custom_asymmetric_eval', 'bagging_freq': 1, 'colsample_bytree': 0.85, 'colsample_bynode': 0.85, 'min_data_per_leaf': 25, 'seed': 1234, 'lambda_l1': 0.01, 'feature_fraction': 0.9, 'num_leaves': 63, 'max_depth': 31, 'min_gain_to_split': 0.1, 'bagging_fraction': 1, 'min_sum_hessian_in_leaf': 0, 'lambda_l2': 10}<jupyter_text>Automated tuning with Bayesian Optimization<jupyter_code>!pip install bayesian-optimization
from bayes_opt import BayesianOptimization
def lgb_evaluate(
                 num_leaves,
                 max_depth,
                 min_sum_hessian_in_leaf,
                 min_gain_to_split,
                 feature_fraction,
                 bagging_fraction,
                 lambda_l2,
                 lambda_l1
                 ):
    params = dict()
    params["learning_rate"] = 0.05
    params["metric"] = lgb_metric
    params["bagging_freq"] = 1
    params['colsample_bytree'] = 0.85
    params['colsample_bynode'] = 0.85
    params['min_data_per_leaf'] = 25
    params["seed"] = 1234
    params['num_leaves'] = int(num_leaves)  
    params['max_depth'] = int(max_depth) 
    params['min_sum_hessian_in_leaf'] = min_sum_hessian_in_leaf
    params['min_gain_to_split'] = min_gain_to_split    
    params['feature_fraction'] = feature_fraction
    params['bagging_fraction'] = bagging_fraction
    params['bagging_freq'] = 1
    params['lambda_l2'] = lambda_l2
    params['lambda_l1'] = lambda_l1   
    params["metric"] = lgb_metric

    lgb_cv = lgb.cv(params,
                    lgb.Dataset(X_train,
                                label=y_train
                            ),
                    num_boost_round=100000,
                    fobj = custom_asymmetric_train,
                    feval = custom_asymmetric_valid,
                    folds = kf.split(X_train, y_train, group), 
                    nfold=5,
                    stratified=False,
                    early_stopping_rounds=50,
                    verbose_eval=False)
    if greater_is_better:
        best_lgb_score = min(lgb_cv['%s-mean' % (lgb_metric)])
    else:
        best_lgb_score = min(lgb_cv['%s-mean' % (lgb_metric)])
    best_lgb_iteration = len(lgb_cv['%s-mean' % (lgb_metric)])
    print (', best_score: %f, best_iteration: %d' % (best_lgb_score, best_lgb_iteration))

    return -best_lgb_score


lgb_BO = BayesianOptimization(lgb_evaluate, 
                             {
                              'num_leaves': (10, 20),
                              'max_depth': (2, 20),
                              'min_sum_hessian_in_leaf': (5, 15),
                              'min_gain_to_split': (0,0),
                              'feature_fraction': (0.2, 0.4),
                              'bagging_fraction': (0.8,1),
                              'lambda_l2': (5, 15),
                              'lambda_l1': (0.1, 5)
                             }
                            )
## I use 5, 20 to save time but you may want to change it to larger numbers,e.g. 8, 30 
lgb_BO.maximize(init_points=5, n_iter=20) 


lgb_BO_scores = pd.DataFrame([p['params'] for p in lgb_BO.res])
lgb_BO_scores['score'] = [p['target'] for p in lgb_BO.res]
lgb_BO_scores = lgb_BO_scores.sort_values(by='score',ascending=False)

lgb_BO_scores.head()
lgb_BO_scores.to_csv("/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/lgb_BO_scores.csv", index=False)
lgb_BO_scores = pd.read_csv("/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/lgb_BO_scores.csv")
lgb_best_params = lgb_BO_scores.T.to_dict().get(lgb_BO_scores.index.values[0])
lgb_best_params['objective'] = "tweedie"
lgb_best_params["tweedie_variance_power"] = 1.1
lgb_best_params['learning_rate'] = 0.01 ## from 0.05 to 0.01
lgb_best_params['seed'] = 1234
lgb_best_params['metric'] = lgb_metric
lgb_best_params['bagging_freq'] = 1

lgb_best_params['num_leaves'] = int(lgb_best_params['num_leaves'])
lgb_best_params['max_depth'] = int(lgb_best_params['max_depth'])

print(lgb_best_params)

lgb_cv = lgb.cv(lgb_best_params,
                lgb.Dataset(X_train,
                            label=y_train
                            ),
                num_boost_round=100000,
                nfold=5,
                stratified=True,
                early_stopping_rounds=50,
                verbose_eval=100)
if greater_is_better:
    best_lgb_score = min(lgb_cv['%s-mean' % (lgb_metric)])
else:
    best_lgb_score = min(lgb_cv['%s-mean' % (lgb_metric)])
best_lgb_iteration = len(lgb_cv['%s-mean' % (lgb_metric)])
print (', best_score: %f, best_iteration: %d' % (best_lgb_score, best_lgb_iteration))
lgb_best_params
lgb_best_params = {'bagging_fraction': 0.9233320946989939,
 'bagging_freq': 1,
 'feature_fraction': 0.3788606347629591,
 'lambda_l1': 4.888559071941389,
 'lambda_l2': 5.073224584824621,
 'learning_rate': 0.01,
 'max_depth': 9,
 'metric': 'rmse',
 'min_gain_to_split': 0.0,
 'min_sum_hessian_in_leaf': 11.07704576482789,
 'num_leaves': 18,
 'objective': 'tweedie',
 'score': -1.9418723740899249,
 'seed': 1234,
 'tweedie_variance_power': 1.1}

best_lgb_iteration = 14125<jupyter_output><empty_output><jupyter_text>single model submissionNow let's retrain the model with the learning rate 0.01 as well as the tuned iterations(num_boost_round) and generate the submission.<jupyter_code>model = lgb.train(lgb_best_params,
                  lgb.Dataset(X_train,
                              label=y_train
                              ),
                  num_boost_round=best_lgb_iteration)
test_data = create_train_data(is_train=False)
date = datetime(2016,4, 25)
day = date + timedelta(days=0)
test_data = test_data[test_data.date >= day - timedelta(days=57)]
# test_data.head()
# test_data.tail()
test_data = create_feature(test_data)
test_data.head()
test_data.tail()
test_data.shape
X_test = test_data[(test_data['date'] >= '2016-04-25') & (test_data['date'] <= '2016-05-22')]
X_test = X_test[train_cols]
X_test.head()
X_test.shape
preds = model.predict(X_test)
test = test_data[(test_data['date'] >= '2016-04-25') & (test_data['date'] <= '2016-05-22')]
test = test[['id', 'date', 'sales']]
test['sales'] = preds
test.to_csv('/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/test.csv', index=False)
predictions = pd.pivot(test, index = 'id', columns = 'date', values = 'sales').reset_index()
predictions.head(10)
predictions.shape
predictions.columns = ['id'] + ['F' + str(i + 1) for i in range(28)]

# read submission
submission = pd.read_csv('/content/m5-forecasting-accuracy/sample_submission.csv')



evaluation_rows = [row for row in submission['id'] if 'evaluation' in row] 
evaluation = submission[submission['id'].isin(evaluation_rows)]

validation = submission[['id']].merge(predictions, on = 'id')
final = pd.concat([validation, evaluation])
final.to_csv('/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/submission_finetune_random5fold.csv', index = False)<jupyter_output><empty_output><jupyter_text>4、 Stacking   Level 1: 6 LightGBM modelsWe will create out-out-fold predictions for 6 LightGBM models (5 automatically tuned models + 1 manually tuned model) and<jupyter_code>from sklearn.metrics import mean_squared_error
kf = model_selection.GroupKFold(n_splits=5)
group = X_train['week'].astype(str) + '_' + X_train['year'].astype(str)
kf_ids = list(kf.split(X_train, y_train, group))
kf_ids[0]
def lgb_stack(params_list, X_train, y_train, X_test,
              kfolds, num_boost_round=10000000, early_stopping_rounds=0, fold='groupkfold', shuffle=True, randome_state=42):
    if fold == 'stratified':
        kf = model_selection.StratifiedKFold(
            n_splits=kfolds, random_state=randome_state, shuffle=shuffle)
        kf_ids = list(kf.split(train_x, train_y))
    elif fold == 'groupkfold':
        kf = model_selection.GroupKFold(n_splits=kfolds)
        group = X_train['week'].astype(str) + '_' + X_train['year'].astype(str)
        kf_ids = list(kf.split(X_train, y_train, group))
    else:
        kf = model_selection.KFold(n_splits=kfolds, random_state=randome_state)
        kf_ids = list(kf.split(train_x, train_y))

    train_x = X_train.values 
    test_x = X_test.values 
    train_y = y_train.values 

    train_blend_x = np.zeros((train_x.shape[0], len(params_list)))
    test_blend_x = np.zeros((test_x.shape[0], len(params_list)))
    blend_scores = np.zeros((kfolds, len(params_list)))

    print("Start blending.")
    for j, params in enumerate(params_list):
        print("Blending model", j+1, params)
        test_blend_x_j = np.zeros((test_x.shape[0]))
        for i, (train_ids, val_ids) in enumerate(kf_ids):
            start = time.time()
            print("Model %d fold %d" % (j+1, i+1))
            train_x_fold = train_x[train_ids]
            train_y_fold = train_y[train_ids]
            val_x_fold = train_x[val_ids]
            val_y_fold = train_y[val_ids]
            print(i, params)

            # Set n_estimators to a large number for early_stopping
            if early_stopping_rounds == 0:
                model = lgb.train(params,
                                  lgb.Dataset(train_x_fold, train_y_fold),
                                  num_boost_round=num_boost_round,
                                  verbose_eval=500
                                  )
                val_y_predict_fold = model.predict(val_x_fold)
                score = np.sqrt(mean_squared_error(val_y_fold, val_y_predict_fold))
                print("Score for Model %d fold %d: %f " % (j+1,i+1,score))
                blend_scores[i, j] = score
                train_blend_x[val_ids, j] = val_y_predict_fold
                test_blend_x_j = test_blend_x_j + model.predict(test_x)
                print("Model %d fold %d finished in %d seconds." % (j+1,i+1, time.time()-start))
            else:
                model = lgb.train(params,
                                  lgb.Dataset(train_x_fold, train_y_fold),
                                  valid_sets=[lgb.Dataset(val_x_fold,
                                                          val_y_fold,
                                                          )],
                                  valid_names=['valid'],
                                  num_boost_round=num_boost_round,
                                  early_stopping_rounds=early_stopping_rounds,
                                  verbose_eval=500
                                  )
                best_iteration = model.best_iteration
                print(model.best_score['valid']['rmse'])
                val_y_predict_fold = model.predict(val_x_fold, num_iteration=best_iteration)
                score = np.sqrt(mean_squared_error(val_y_fold, val_y_predict_fold))
                print("Score for Model %d fold %d: %f " % (j+1,i+1,score))
                blend_scores[i, j] = score
                train_blend_x[val_ids, j] = val_y_predict_fold
                test_blend_x_j = test_blend_x_j + \
                    model.predict(test_x, num_iteration=best_iteration)
                print("Model %d fold %d finished in %d seconds." % (j+1,i+1, time.time()-start))  
            print(time.time()-start)
        test_blend_x[:,j] = test_blend_x_j/kfolds
        print("Score for model %d is %f" % (j+1, np.mean(blend_scores[:, j])))
    return train_blend_x, test_blend_x, blend_scores
lgb_params_list = []
lgb_metric = 'rmse'
for i in range(5):
    lgb_param = lgb_BO_scores.iloc[i].to_dict()
    lgb_param.pop('score')
    lgb_param['max_depth'] = int(lgb_param['max_depth'])
    lgb_param['num_leaves'] = int(lgb_param['num_leaves'])
    lgb_param["objective"] = "tweedie"
    lgb_param["tweedie_variance_power"] = 1.1
    lgb_param['learning_rate'] = 0.01
    lgb_param['metric'] = lgb_metric
    lgb_param['seed'] = 42
    lgb_params_list.append(lgb_param)

## Best manual-tuned parameters
lgb_params_list.append(best_lgb_params)    
print(lgb_params_list)
import pickle
train_blend_x_lgb_01, test_blend_x_lgb_01, blend_scores_lgb_01 = lgb_stack(lgb_params_list,
                                                                           X_train,
                                                                           y_train,
                                                                           X_test,
                                                                           early_stopping_rounds=200,
                                                                           kfolds=5)


pickle.dump(train_blend_x_lgb_01, open(
    '/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/train_blend_x_lgb_01.pkl', 'wb'))

pickle.dump(test_blend_x_lgb_01, open(
    '/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/test_blend_x_lgb_01.pkl', 'wb'))
np.sqrt(mean_squared_error(y_train,train_blend_x_lgb_01.mean(axis=1)))
[np.sqrt(mean_squared_error(y_train,train_blend_x_lgb_01[:,n])) for n in range(6)]
train_blend_x_lgb_01.shape<jupyter_output><empty_output><jupyter_text>Level 2 stacking I : ridge regressionWe will be using the predictions (Level 1 features) to train another model as Level 2 which will then be used to make the actual predictions.<jupyter_code># Stacke level 1 features
train_blend_x = train_blend_x_lgb_01
test_blend_x = test_blend_x_lgb_01
from sklearn.metrics import make_scorer
def my_rmse(y_true, y_preds):
    return np.sqrt(mean_squared_error(y_true, y_preds))

my_rmse_score = make_scorer(my_rmse, greater_is_better=False)
from sklearn.linear_model import Ridge,ElasticNet, SGDRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import GridSearchCV
def search_model(train_x, train_y, est, param_grid, n_jobs, cv, refit=False):
    ##Grid Search for the best model
    model = GridSearchCV(estimator  = est,
                                     param_grid = param_grid,
                                     scoring    = my_rmse_score,
                                     verbose    = 10,
                                     n_jobs  = n_jobs,
                                     iid        = True,
                                     refit    = refit,
                                     cv      = cv)
    # Fit Grid Search Model
    model.fit(train_x, train_y)
    print("Best score: %0.3f" % model.best_score_)
    print("Best parameters set:", model.best_params_)
    print("Scores:", model.cv_results_)
    return model
param_grid = {
              "alpha":[0.001,0.01,0.1,1,10,30,100]
              }
model = search_model(train_blend_x,
                     y_train, 
                     Ridge(), 
                     param_grid, 
                     n_jobs=1, 
                     cv=kf.split(X_train, y_train, group), 
                     refit=True)   

print ("best alpha:", model.best_params_)


preds_ridge = model.predict(test_blend_x)
preds_ridge
ridge_test = test_data[(test_data['date'] >= '2016-04-25') & (test_data['date'] <= '2016-05-22')]
ridge_test = ridge_test[['id', 'date', 'sales']]
ridge_test ['sales'] = preds_ridge
ridge_test.head()
ridge_predictions = pd.pivot(ridge_test, index = 'id', columns = 'date', values = 'sales').reset_index()
ridge_predictions.shape
ridge_predictions.columns = ['id'] + ['F' + str(i + 1) for i in range(28)]

# read submission
submission = pd.read_csv('/content/m5-forecasting-accuracy/sample_submission.csv')



evaluation_rows = [row for row in submission['id'] if 'evaluation' in row] 
evaluation = submission[submission['id'].isin(evaluation_rows)]

validation = submission[['id']].merge(ridge_predictions, on = 'id')
ridge_final = pd.concat([validation, evaluation])
ridge_final.to_csv('/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/sub_stacking_ridge.csv', index = False)<jupyter_output><empty_output><jupyter_text>Level 2 stacking II: LightGBM regressionStacking can also be done by with both level 1 predictions and original features.<jupyter_code># len(X_train.values[0])
# train_blend_x_lgb_01
# sparse.coo_matrix(train_blend_x_lgb_01)
# sparse.hstack([X_train.values, sparse.coo_matrix(train_blend_x_lgb_01)]).tocsr()
train_blend_x = sparse.hstack([X_train.values, sparse.coo_matrix(train_blend_x_lgb_01)]).tocsr()
test_blend_x = sparse.hstack([X_test.values, sparse.coo_matrix(test_blend_x_lgb_01)]).tocsr()

print (train_blend_x.shape, train_blend_x.shape)
lgb_params = {'learning_rate': 0.05, 'metric': 'rmse',
              'bagging_freq': 1, 'seed': 1234, 'objective': 'regression',
              'num_leaves': 7, 'verbose': 1,
              'max_depth': 5, 'min_gain_to_split': 0,
              'feature_fraction': 0.1,
              'bagging_fraction': 0.9,
              'min_sum_hessian_in_leaf': 1,
              'lambda_l2': 0, 'lambda_l1': 0
              }


lgb_cv = lgb.cv(lgb_params,
                lgb.Dataset(train_blend_x,
                            label=y_train
                            ),
                num_boost_round=100000,
                nfold=5,
                folds=kf.split(X_train, y_train, group), 
                early_stopping_rounds=50,
                verbose_eval=500)

best_lgb_score = min(lgb_cv['rmse-mean'])
best_lgb_iteration = len(lgb_cv['rmse-mean'])
print(', best_score: %f, best_iteration: %d' %
      (best_lgb_score, best_lgb_iteration))
model = lgb.train(lgb_params,
            lgb.Dataset(train_blend_x,
                        label=y_train
                        ),
            num_boost_round=best_lgb_iteration)
preds_lgb = model.predict(test_blend_x)
lgb_test = test_data[(test_data['date'] >= '2016-04-25') & (test_data['date'] <= '2016-05-22')]
lgb_test = lgb_test[['id', 'date', 'sales']]
lgb_test ['sales'] = preds_lgb
lgb_test.head()
lgb_predictions = pd.pivot(lgb_test, index = 'id', columns = 'date', values = 'sales').reset_index()
lgb_predictions.shape
lgb_predictions.columns = ['id'] + ['F' + str(i + 1) for i in range(28)]

# read submission
submission = pd.read_csv('/content/m5-forecasting-accuracy/sample_submission.csv')



evaluation_rows = [row for row in submission['id'] if 'evaluation' in row] 
evaluation = submission[submission['id'].isin(evaluation_rows)]

validation = submission[['id']].merge(lgb_predictions, on = 'id')
lgb_final = pd.concat([validation, evaluation])
lgb_final.to_csv('/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/sub_stacking_lgb.csv', index = False)<jupyter_output><empty_output><jupyter_text>Averaging level 2 predictionsWe can also create a quick ensemble by averaging level 2 predictions:<jupyter_code>preds_l2 = preds_ridge*0.5 + preds_lgb*0.5
l2_test = test_data[(test_data['date'] >= '2016-04-25') & (test_data['date'] <= '2016-05-22')]
l2_test = l2_test[['id', 'date', 'sales']]
l2_test ['sales'] = preds_l2
l2_test.head()
l2_predictions = pd.pivot(l2_test, index = 'id', columns = 'date', values = 'sales').reset_index()
l2_predictions.columns = ['id'] + ['F' + str(i + 1) for i in range(28)]

# read submission
submission = pd.read_csv('/content/m5-forecasting-accuracy/sample_submission.csv')



evaluation_rows = [row for row in submission['id'] if 'evaluation' in row] 
evaluation = submission[submission['id'].isin(evaluation_rows)]

validation = submission[['id']].merge(l2_predictions, on = 'id')
l2_final = pd.concat([validation, evaluation])
l2_final.to_csv('/content/drive/My Drive/colab/kaggle_M5/dataapplication/output/sub_stacking_l2.csv', index = False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13241970021413277, 'fraction_numerical': 0.03283369022127052, 'mean_word_length': 2.9974891577265463, 'pattern_counts': {'":': 29, '<': 58, '<?xml version=': 0, '>': 56, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 39, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '508814', 'license': 'MIT', 'n_tokens_mistral': 13841, 'n_tokens_neox': 12847, 'n_words': 2204, 'path': 'M5_baseline_groupKfold_custom_lossfunction_20200525.ipynb', 'repo_name': 'xfsm1912/kaggle_M5-Forecasting_Accuracy', 'size': '35025'}