?
目錄
一、題目概要
二、導入包和數據集
三、數據處理
四、描述性分析
五、探索性數據分析
六、模型一：線性回歸
七、模型2：隨機森林
一、題目概要
在Kaggle競賽中，要求我們應用時間序列預測，根據厄瓜多爾大型雜貨零售商Corporación Favorita的數據預測商店銷售情況，建立一個模型，準確地預測在不同商店銷售的商品的單位銷量。準確的預測可以減少與庫存過多相關的食物浪費，提高客戶滿意度。
在六個可用的數據文件中，我們分析了其中的三個，即訓練、測試和存儲。雖然我們在這個項目中沒有研究每日油價或假日事件的影響，但我們希望在這門課之外花更多的時間來深入學習和成長。
在我們的分析中，我們探索了兩種不同的時間序列預測模型:線性回歸和隨機森林。通過準備線性回歸的數據，我們發現了一些有趣的見解，包括周末的銷售增長，公共部門支付工資時的銷售增長，以及11月和12月的銷售增長。我們也注意到2014年和2015年的銷量大幅下降。這兩種增加和減少都可能是由于商店促銷、假期、油價或世界事件，我們無法在分析中調查。
對于線性回歸，我們刪除了沒有在該特定商店銷售的產品族，對商店進行聚類，并將銷售較低的產品分組到一個產品族中。我們調查并去除異常值，并評估訓練數據的季節性。
對于隨機森林，我們處理分類變量并去除異常值。
線性回歸被證明是低效的，表現出指數性質，而隨機森林被證明是一個更容易實現的模型。然而，如果沒有線性回歸的數據處理，我們將無法發現我們所做的洞察，因為隨機森林是一個黑盒模型。
我們還包含了一些關于特征重要性、超參數優化和殘差圖的進一步發現的最終想法。
綜上所述，隨機森林是預測商店銷售額的較好模型。
#This Python 3 environment comes with many helpful analytics libraries installed
#It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python
#For example, here's several helpful packages to load
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
#Input data files are available in the read-only "../input/" directory
#For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory
import os
for dirname, _, filenames in os.walk('/kaggle/input'):
for filename in filenames:
print(os.path.join(dirname, filename))
#You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using "Save & Run All"
#You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session
/kaggle/input/store-sales-time-series-forecasting/oil.csv
/kaggle/input/store-sales-time-series-forecasting/sample_submission.csv
/kaggle/input/store-sales-time-series-forecasting/holidays_events.csv
/kaggle/input/store-sales-time-series-forecasting/stores.csv
/kaggle/input/store-sales-time-series-forecasting/train.csv
/kaggle/input/store-sales-time-series-forecasting/test.csv
/kaggle/input/store-sales-time-series-forecasting/transactions.csv
二、導入包和數據集
#Import packages
#BASE
# ------------------------------------------------------
import numpy as np
import pandas as pd
import os
import gc
import warnings
#Machine Learning
# ------------------------------------------------------
import statsmodels.api as sm
import sklearn
#Data Visualization
# ------------------------------------------------------
#import altair as alt
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
warnings.filterwarnings('ignore')
from plotly.offline import init_notebook_mode, iplot
init_notebook_mode(connected=True)
#Import datasets
train = pd.read_csv("/kaggle/input/store-sales-time-series-forecasting/train.csv",parse_dates=['date'])
test = pd.read_csv("/kaggle/input/store-sales-time-series-forecasting/test.csv",parse_dates=['date'])
stores = pd.read_csv("/kaggle/input/store-sales-time-series-forecasting/stores.csv")
我們將檢查以下數據集中每個列的數據類型。為了使用我們導入的包執行時間序列預測，我們必須確保已將日期解析為日期。稍后，我們將把“對象”數據類型轉換為類別。我們還將查看我們的數據集，看看是否需要進行任何清理或操作。
train.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 3000888 entries, 0 to 3000887
Data columns (total 6 columns):
| Column | Dtype | |
|---|---|---|
| 0 | id | int64 |
| 1 | date | datetime64[ns] |
| 2 | store_nbr | int64 |
| 3 | family | object |
| 4 | sales | float64 |
| 5 | onpromotion | int64 |
dtypes: datetime64ns, float64(1), int64(3), object(1)
memory usage: 137.4+ MB
train.head()
id date store_nbr family sales onpromotion
0 0 2013-01-01 1 AUTOMOTIVE 0.0 0
1 1 2013-01-01 1 BABY CARE 0.0 0
2 2 2013-01-01 1 BEAUTY 0.0 0
3 3 2013-01-01 1 BEVERAGES 0.0 0
4 4 2013-01-01 1 BOOKS 0.0 0
stores.head()
store_nbr city state type cluster
0 1 Quito Pichincha D 13
1 2 Quito Pichincha D 13
2 3 Quito Pichincha D 8
3 4 Quito Pichincha D 9
4 5 Santo Domingo Santo Domingo de los Tsachilas D 4
三、數據處理
從train.csv和transactions.csv中:
date-記錄數據的日期。
store_nbr -標識銷售產品的商店。
family-標識所售產品的類型。
sales -給出給定日期某一特定商店某一產品系列的總銷售額。小數值是可能的。
onpromotion -給出在給定日期某商店促銷的產品族的總數量。
transactions -在給定日期在商店中發生的交易總數。
刪除不銷售特定系列產品的商店的銷售
在快速查看我們的訓練數據集之后，我們可以看到有很多零。有些商店可能不銷售某些產品，因為它們不是該產品的合適商店。在這種情況下，我們將刪除這些值，當預測，他們不應該有任何銷售。
zeros = train.groupby(['id', 'store_nbr', 'family']).sales.sum().reset_index().sort_values(['family','store_nbr'])
zeros = zeros[zeros.sales == 0]
zeros
id store_nbr family sales
0 0 1 AUTOMOTIVE 0.0
10692 10692 1 AUTOMOTIVE 0.0
30294 30294 1 AUTOMOTIVE 0.0
40986 40986 1 AUTOMOTIVE 0.0
53460 53460 1 AUTOMOTIVE 0.0
... ... ... ... ...
2981153 2981153 54 SEAFOOD 0.0
2984717 2984717 54 SEAFOOD 0.0
2986499 2986499 54 SEAFOOD 0.0
2993627 2993627 54 SEAFOOD 0.0
2998973 2998973 54 SEAFOOD 0.0
#full outer joining the tables and removing the rows where they match to get rid of the zeros
join = train.merge(zeros[zeros.sales == 0].drop("sales",axis = 1), how='outer', indicator=True)
train1 = join[~(join._merge == 'both')].drop(['id', '_merge'], axis = 1).reset_index()
train1 = train1.drop(['index', 'onpromotion'], axis=1)
train1
date store_nbr family sales
0 2013-01-01 25 BEAUTY 2.000
1 2013-01-01 25 BEVERAGES 810.000
2 2013-01-01 25 BREAD/BAKERY 180.589
3 2013-01-01 25 CLEANING 186.000
4 2013-01-01 25 DAIRY 143.000
... ... ... ... ...
2061753 2017-08-15 9 POULTRY 438.133
2061754 2017-08-15 9 PREPARED FOODS 154.553
2061755 2017-08-15 9 PRODUCE 2419.729
2061756 2017-08-15 9 SCHOOL AND OFFICE SUPPLIES 121.000
2061757 2017-08-15 9 SEAFOOD 16.000
集群存儲
訓練數據集就是我們用來創建模型的數據。共有 54 家商店，每家商店有 33 個產品系列。由于這是一個擁有 300 多萬行的大型數據集，將商店分組將大大減少我們后續模型的計算量。幸運的是，"商店 "數據集已經為我們將商店聚類為 17 個聚類。為了簡單起見和時間限制，我們還選擇忽略 "促銷 "變量。
def group_clusters (df) :
#left join train and stores on store number
jointr = df.merge(stores, on='store_nbr', how='left', indicator=False)
#replacing all store numbers with their cluster and grouping them by cluster
grouped = jointr.groupby(['date', 'cluster', 'family']).sum('sales').reset_index()
#removing columns id, store_nbr, and type as they are aggregated values with no significance
grouped = grouped.drop(['store_nbr'], axis=1)
return grouped
grouped = group_clusters (train1)
grouped
date cluster family sales
0 2013-01-01 1 BEAUTY 2.000
1 2013-01-01 1 BEVERAGES 810.000
2 2013-01-01 1 BREAD/BAKERY 180.589
3 2013-01-01 1 CLEANING 186.000
4 2013-01-01 1 DAIRY 143.000
... ... ... ... ...
742573 2017-08-15 17 PLAYERS AND ELECTRONICS 25.000
742574 2017-08-15 17 POULTRY 686.941
742575 2017-08-15 17 PREPARED FOODS 91.976
742576 2017-08-15 17 PRODUCE 5031.190
742577 2017-08-15 17 SEAFOOD 52.876
產品系列分組
接下來，我們將對這些系列進行探討，以幫助我們更好地了解哪些系列對商店的總銷售額貢獻最大，從而有可能對它們進行分組。
#group by 'family', then sort by sales, and create a column for aggregate sales percent
temp = grouped.groupby('family').sum('sales').reset_index().sort_values(by='sales', ascending=False)
#aggregated sales
temp = temp[['family','sales']]
temp['percent']=(temp['sales']/temp['sales'].sum())
temp['percent'] = temp['percent'].apply(lambda x: f'{x:.0%}')
temp['cumulative']=(temp['sales']/temp['sales'].sum()).cumsum()
temp['cumulative'] = temp['cumulative'].apply(lambda x: f'{x:.0%}')
temp.head()
family sales percent cumulative
12 GROCERY I 3.434627e+08 32% 32%
3 BEVERAGES 2.169545e+08 20% 52%
30 PRODUCE 1.227047e+08 11% 64%
7 CLEANING 9.752129e+07 9% 73%
8 DAIRY 6.448771e+07 6% 79%
上表顯示了對所有商店和日期的銷售額貢獻最大的五個系列：食品雜貨 I、飲料、農產品、清潔和奶制品。食品雜貨 I 在總銷售額中所占比例最大，為 32%。前五大產品系列占總銷售額的 79%。下圖更清楚地顯示了各產品系列的累計銷售百分比。
#plot ranked category sales
fig1 = px.bar(temp, x="family",y="sales",title = "Sales",text="cumulative")
fig1.show()
我們決定將其余產品系列歸入一個名為 "其他 "的系列，以減少系列類別的數量。這將使我們以后使用虛擬變量進行回歸時更加方便。
#list of the top 5 families
top5 = ['GROCERY I','BEVERAGES','PRODUCE','CLEANING','DAIRY']
#removing the top 5 families so we can get the list of remaining families
tmp = grouped[~grouped['family'].isin(top5)]
#the list of families that we want to group into 'OTHERS'
tmp['family'].unique()
#replace the above list with 'OTHERS'
trainc = grouped.copy()
trainc['family'] = grouped['family'].replace(['AUTOMOTIVE', 'BABY CARE', 'BEAUTY', 'BOOKS', 'BREAD/BAKERY',
'CELEBRATION', 'DELI', 'EGGS', 'FROZEN FOODS', 'GROCERY II',
'HARDWARE', 'HOME AND KITCHEN I', 'HOME AND KITCHEN II',
'HOME APPLIANCES', 'HOME CARE', 'LADIESWEAR', 'LAWN AND GARDEN',
'LINGERIE', 'LIQUOR,WINE,BEER', 'MAGAZINES', 'MEATS',
'PERSONAL CARE', 'PET SUPPLIES', 'PLAYERS AND ELECTRONICS',
'POULTRY', 'PREPARED FOODS', 'SCHOOL AND OFFICE SUPPLIES',
'SEAFOOD'],'OTHERS')
newtrain = trainc.groupby(['date', 'cluster', 'family']).sum('sales').reset_index()
newtrain
date cluster family sales
0 2013-01-01 1 BEVERAGES 810.000000
1 2013-01-01 1 CLEANING 186.000000
2 2013-01-01 1 DAIRY 143.000000
3 2013-01-01 1 GROCERY I 700.000000
4 2013-01-01 1 OTHERS 672.618999
... ... ... ... ...
165214 2017-08-15 17 CLEANING 1357.000000
165215 2017-08-15 17 DAIRY 1377.000000
165216 2017-08-15 17 GROCERY I 4756.000000
165217 2017-08-15 17 OTHERS 3773.369000
165218 2017-08-15 17 PRODUCE 5031.190000
總之，我們現在已經刪除了從未有過銷售額的產品系列，將商店聚類為 17 個群組，并將銷售額較低的產品系列歸入 "其他 "群組。這樣，我們就大大減少了數據集中的記錄數量，從而改進了計算，方便了稍后的回歸。
四、描述性分析
newtrain.groupby('family').describe()['sales'].applymap(lambda x: f"{x:.2f}")
count mean min 25% 50% 75% max std
family
BEVERAGES 28398.00 7639.78 242.00 2990.25 5536.50 10469.50 58848.00 6686.87
CLEANING 28399.00 3433.97 186.00 1462.00 2846.00 4639.00 22544.00 2545.52
DAIRY 28398.00 2270.85 66.00 979.25 1838.00 3022.00 14339.00 1764.45
GROCERY I 28398.00 12094.61 700.00 5012.50 9925.00 16164.00 138535.00 9418.56
OTHERS 28399.00 8046.55 596.89 3704.74 6454.25 10780.98 105366.90 6125.40
PRODUCE 23227.00 5282.85 1.00 82.00 4045.53 7953.38 31994.97 5711.49
上表顯示了產品系列的簡要說明。我們將繪制產品系列分組圖，以確定數據的形狀以及是否存在異常值。
# Create a subplot grid with 3 rows and 2 columns
fig, axes = plt.subplots(3, 2, figsize=(12, 12))
# Flatten axes for easier indexing
axes = axes.ravel()
# Unique family categories
families = newtrain['family'].unique()
# Plot sales boxplots for different family categories
for i, family in enumerate(families):
filtered_data = newtrain[newtrain['family'] == family]
sns.boxplot(data=filtered_data, x='sales', ax=axes[i])
axes[i].set_xlabel('Sales')
axes[i].set_title(f'Boxplot of Sales for {family} Family')
# Hide extra subplots
for j in range(len(families), len(axes)):
axes[j].axis('off')
plt.tight_layout() # Automatically adjust subplot layout to prevent overlap
plt.show()
根據上述方框圖，我們可以得出結論：它們高度偏斜。雖然直觀地認為雜貨店偶爾也會有一些高銷量，但為了簡單起見，我們還是決定按照方框圖作為剔除異常值的理由。
我們還注意到，許多產品系列的銷售額為零。雖然這可能是由于商店關閉造成的，但我們也不能忽視人為失誤的可能性。然而，這些零銷售額的存在可能會拉低平均值。因此，我們不會完全刪除方框圖之外的所有異常值。
刪除異常值
#function for removing outliers
def remove_outliers (df) :
# Calculate the first quartile (Q1)
q1 = df.groupby('family')['sales'].transform('quantile', 0.25)
# Calculate the third quartile (Q3)
q3 = df.groupby('family')['sales'].transform('quantile', 0.75)
# Calculate the Interquartile Range (IQR)
IQR = q3 - q1
# Define the lower and upper bounds for outliers
lbound = q1 - 1.5 * IQR
ubound = q3 + 1.5 * IQR
# Filter the dataset to remove outliers
no_outliers = df[~((df['sales'] < lbound) | (df['sales'] > ubound))]
return no_outliers
no_outliers = remove_outliers (newtrain)
no_outliers
date cluster family sales
0 2013-01-01 1 BEVERAGES 810.000000
1 2013-01-01 1 CLEANING 186.000000
2 2013-01-01 1 DAIRY 143.000000
3 2013-01-01 1 GROCERY I 700.000000
4 2013-01-01 1 OTHERS 672.618999
... ... ... ... ...
165214 2017-08-15 17 CLEANING 1357.000000
165215 2017-08-15 17 DAIRY 1377.000000
165216 2017-08-15 17 GROCERY I 4756.000000
165217 2017-08-15 17 OTHERS 3773.369000
165218 2017-08-15 17 PRODUCE 5031.190000
# plot new boxplots for no_outliers
# Create a subplot grid with 3 rows and 2 columns
fig, axes = plt.subplots(3, 2, figsize=(12, 12))
# Flatten axes for easier indexing
axes = axes.ravel()
# Unique family categories
families = no_outliers['family'].unique()
# Plot sales boxplots for different family categories
for i, f in enumerate(families):
filtered = no_outliers[no_outliers['family'] == f]
sns.boxplot(data=filtered, x='sales', ax=axes[i])
axes[i].set_xlabel('Sales')
axes[i].set_title(f'Boxplot of Sales for {f} Family')
# Hide extra subplots
for j in range(len(families), len(axes)):
axes[j].axis('off')
plt.tight_layout() # Automatically adjust subplot layout to prevent overlap
plt.show()
五、探索性數據分析
現在，我們可以繪制前五大產品系列和 "其他 "產品系列的圖表，以確定我們應該使用哪種回歸方法。
wtrain = no_outliers.set_index("date").groupby("family").resample("W").sales.sum().reset_index()
top6 = ['GROCERY I','BEVERAGES','PRODUCE','CLEANING','DAIRY', 'OTHERS']
cond = wtrain['family'].isin(top6)
px.line(wtrain[cond], x = "date", y = "sales", color = "family", title = "Weekly Total Sales by Family")
上圖顯示了明顯的周和月季節性跡象，以及所有產品系列的增長趨勢。有趣的是，我們可以看到 2014 年以及 2015 年和 2017 年 8 月的兩次銷售大跌。這些下滑可能是由于世界重大事件、節假日或油價造成的，雖然我們已經獲得了這些數據，但我們選擇暫時忽略它們。也許將來可以對它們進行進一步調查。
然而，農產品銷售卻呈現出一種奇怪的模式，似乎沒有任何銷售。下面的圖表將對此進行進一步研究。
#sns.set(rc={"figure.figsize":(17, 10)})
dtrain = no_outliers.set_index("date").groupby("family").resample("D").sales.sum().reset_index()
gfig = sns.lineplot(dtrain[dtrain['family']=='GROCERY I'],x='date',y='sales').set_title('GROCERY I SALES')
上圖顯示了雜貨 I 在所有商店和日期的銷售情況。它清楚地描述了每天的季節性和持續增長的趨勢。每年年初的銷售額似乎為零，這可能與年初的節日或活動有關，但我們不會在本報告中對此進行調查。
bfig = sns.lineplot(dtrain[dtrain['family']=='BEVERAGES'],x='date',y='sales').set_title('BEVERAGES SALES')
飲料銷售也呈現出明顯的季節性，但在 2014 年和 2015 年，季節性隨著趨勢的變化而增強，從 2015 年年中開始似乎趨于穩定。這可能是我們以后需要調查的問題。與雜貨 I 的銷售一樣，年初似乎沒有銷售。
pfig = sns.lineplot(dtrain[dtrain['family']=='PRODUCE'],x='date',y='sales').set_title('PRODUCE SALES')
plt.xticks(rotation=45)
plt.show()
農產品銷售的表現最為奇怪。2013 年的銷售量似乎很低，多次出現銷售高峰，然后再次下降。從 2015 年年中開始，銷售量繼續與季節性和趨勢保持一致。
ptrain = dtrain[(dtrain['date'].dt.year == 2013) & (dtrain['family'] == 'PRODUCE')]
pfig2013 = sns.lineplot(ptrain,x='date',y='sales').set_title('PRODUCE SALES in 2013')
plt.xticks(rotation=45)
plt.show()
在上圖中，我們可以更近距離地觀察 2013 年的農產品銷售情況。雖然銷售量較低，但仍顯示出明顯的季節性。
cfig = sns.lineplot(dtrain[dtrain['family']=='CLEANING'],x='date',y='sales').set_title('CLEANING SALES')
dfig = sns.lineplot(dtrain[dtrain['family']=='DAIRY'],x='date',y='sales').set_title('DAIRY SALES')
ofig = sns.lineplot(dtrain[dtrain['family']=='OTHERS'],x='date',y='sales').set_title('OTHERS SALES')
家庭清潔、乳制品和其他產品的銷售表現出明顯的每日季節性和一致的趨勢。
下面兩個數字簡單地說明了每個類別的總銷售額或平均銷售額。
# plot a bar chart
sales_by_family = newtrain.groupby('family')['sales'].sum().reset_index().sort_values(by='sales', ascending=False)
sns.barplot(data=sales_by_family, x='family', y='sales').set_title("Sum Sales by Category")
Text(0.5, 1.0, 'Sum Sales by Category')
# plot a bar chart
sales_by_family2 = newtrain.groupby('family')['sales'].mean().reset_index().sort_values(by='sales', ascending=False)
sns.barplot(data=sales_by_family2, x='family', y='sales').set_title("Average Sales by Category")
Text(0.5, 1.0, 'Average Sales by Category')
由于存在一些季節性因素，我們可以做一些工程，看看數據中的“天數”是否有意義.
六、模型一：線性回歸
現在我們將嘗試使用假人對以下數據集執行線性回歸。由于我們分別處理每個聚類，我們將對每個聚類運行線性回歸。
#function for running linear regression
def lin_reg (df, c):
y = df['sales']
#exclude family_6 and month_12
x = df[[
'4d_within_pay',
'weekend',
'family_1',
'family_2',
'family_3',
'family_4',
'family_5',
'month_1',
'month_2',
'month_3',
'month_4',
'month_5',
'month_6',
'month_7',
'month_8',
'month_9',
'month_10',
'month_11',
]]
x = sm.add_constant(x)
model = sm.OLS(y, x).fit()
#plotting the residual and normal probability plots side by side
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sns.residplot(x=model.fittedvalues, y=model.resid, lowess=True, line_kws={"color": "red"}, ax=axes[0])
axes[0].set_xlabel('Predicted Sales')
axes[0].set_ylabel('Residuals')
axes[0].set_title(f'Residual Plot for Cluster {c}')
sm.qqplot(model.resid, line='s', ax=axes[1])
axes[1].set_xlabel('Predicted Sales')
axes[1].set_ylabel('Actual Sales')
axes[1].set_title(f'Normal Probability Plot for Cluster {c}')
plt.tight_layout
plt.show
return model
#cycling through the clusters and running linear regression on each one
for c in range(1,18):
result = lin_reg(train_m1[train_m1['cluster'] == c], c)
print("\nRegression Results for Cluster " + str(c) + "\n")
print(result.summary())
模型1的目的是觀察銷售是否隨時間呈線性關系。由于這是最基本的回歸方法，這將幫助我們確定下一步應該做什么。如果模型表現良好，我們就可以使用線性回歸來預測銷售。然而，這個模型被證明是很差的，我們可以在數據清洗過程中看到這一點，調整r平方，殘差和正態概率圖。
首先，數據處理花費了大量的時間。如上所述，有54家商店和33個產品系列，這已經在我們的數據中產生了太多不同的變量，迫使我們進行聚類和分組。還有季節性因素需要處理，為了讓我們用多元回歸來做這個，我們必須創造假人，這又創造了太多的變量。為了不違反簡約原則，擁有太多變量并不是一個好主意，而且僅僅為了準備線性回歸的數據集而需要進行的處理量是低效的。由于季節性似乎是加性的，我們可以探索指數平滑方法，如Holt-Winters，但我們選擇不在這里繼續研究。
由于我們的模型包括許多預測變量，我們將查看調整后的r平方值，這對于使用線性回歸運行的每個集群來說都很差。雖然一些集群輸出的調整后的r平方值高于0.7，這被認為是好的，但大多數集群的r平方值在0.5左右。這意味著大約50%的時間，我們的模型沒有捕捉到銷售數量的變化。然而，所有運行中預測變量的p值均小于0.05，表明它們具有顯著性。
首先看一下集群1的殘差圖，這些點是隨機的，正態概率圖幾乎是一條直線。這些都是一個好模型的代表。然而，當我們繼續觀察殘差圖時，我們注意到它們之間確實存在一種模式。雖然我們不能根據我們對不同類型回歸的知識來定義這個模式，但這個模式確實證明了我們使用的線性回歸模型不是一個很好的擬合。一些正態概率圖也表現出指數增長，最值得注意的是集群12、15和16的圖。在這種情況下，我們可以嘗試通過對銷售額取對數來使用指數趨勢模型，但由于時間限制以及我們知道回歸模型在這里不是很適合的事實，我們選擇在此時不探索指數趨勢。
我們已經得出結論，由于變量的數量，調整后的r平方和殘差圖，模型1不是預測該數據集中商店銷售的好模型。然而，我們確實注意到，這個模型是在我們從分析中省略了促銷、假期和油價等其他變量后創建的。如果我們選擇包括這些，我們可能會得到更好的結果，但在這個時候，我們不會用它來做任何預測，而是將轉向另一種可能給我們更好的結果的技術。
七、模型2：隨機森林
我們在隨機森林中創建預測商店銷售模型的第二種方法。總之，隨機森林是一種監督學習和集成算法，是決策樹的擴展。它既可以用于分類模型，也可以用于回歸模型，具有較高的精度，并且能夠有效地處理大型數據集。然而，這是一個黑盒模型，因為我們不能真正解釋隨機森林中的單個決策樹。
Random Forest能夠處理數據中的許多商店、產品系列和季節性，因此我們不需要擔心分組、創建假人和對銷售額取對數。因此，我們可以回頭查看未處理的數據。值得慶幸的是，隨機森林有一個包，我們可以導入它，而不必分解決策樹的細節。
from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
#from sklearn.linear_model import LinearRegression
#from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import GridSearchCV,RandomizedSearchCV
id date store_nbr family sales onpromotion
0 0 2013-01-01 1 AUTOMOTIVE 0.000 0
1 1 2013-01-01 1 BABY CARE 0.000 0
2 2 2013-01-01 1 BEAUTY 0.000 0
3 3 2013-01-01 1 BEVERAGES 0.000 0
4 4 2013-01-01 1 BOOKS 0.000 0
... ... ... ... ... ... ...
3000883 3000883 2017-08-15 9 POULTRY 438.133 0
3000884 3000884 2017-08-15 9 PREPARED FOODS 154.553 1
3000885 3000885 2017-08-15 9 PRODUCE 2419.729 148
3000886 3000886 2017-08-15 9 SCHOOL AND OFFICE SUPPLIES 121.000 8
3000887 3000887 2017-08-15 9 SEAFOOD 16.000 0
# use original train dataset and drop all the zero sales value
train_rf = join[~(join._merge == 'both')].drop(['id', '_merge'], axis = 1).reset_index()
train_rf = train_rf.drop(['index'], axis=1)
train_rf
date store_nbr family sales onpromotion
0 2013-01-01 25 BEAUTY 2.000 0
1 2013-01-01 25 BEVERAGES 810.000 0
2 2013-01-01 25 BREAD/BAKERY 180.589 0
3 2013-01-01 25 CLEANING 186.000 0
4 2013-01-01 25 DAIRY 143.000 0
... ... ... ... ... ...
2061753 2017-08-15 9 POULTRY 438.133 0
2061754 2017-08-15 9 PREPARED FOODS 154.553 1
2061755 2017-08-15 9 PRODUCE 2419.729 148
2061756 2017-08-15 9 SCHOOL AND OFFICE SUPPLIES 121.000 8
2061757 2017-08-15 9 SEAFOOD 16.000 0
rotrain = remove_outliers(train_rf)
rotrain['codetime'] = (rotrain['date'] - rotrain['date'].min()).dt.days
rotrain['month'] = rotrain['date'].dt.month
rotrain['day'] = rotrain['date'].dt.day
rotrain['day_of_week'] = rotrain['date'].dt.dayofweek
train2 = rotrain.copy()
分類變量預處理與異常值去除
date store_nbr family sales onpromotion codetime month day day_of_week
0 2013-01-01 25 BEAUTY 2.000 0 0 1 1 1
1 2013-01-01 25 BEVERAGES 810.000 0 0 1 1 1
2 2013-01-01 25 BREAD/BAKERY 180.589 0 0 1 1 1
3 2013-01-01 25 CLEANING 186.000 0 0 1 1 1
4 2013-01-01 25 DAIRY 143.000 0 0 1 1 1
... ... ... ... ... ... ... ... ... ...
2061752 2017-08-15 9 PLAYERS AND ELECTRONICS 6.000 0 1687 8 15 1
2061753 2017-08-15 9 POULTRY 438.133 0 1687 8 15 1
2061754 2017-08-15 9 PREPARED FOODS 154.553 1 1687 8 15 1
2061755 2017-08-15 9 PRODUCE 2419.729 148 1687 8 15 1
2061757 2017-08-15 9 SEAFOOD 16.000 0 1687 8 15 1
在模型1中，我們在刪除異常值之前將產品族分組在一起，但是現在我們將刪除異常值，同時保持所有產品族的原樣。
隨機森林包需要將所有的分類變量轉換為數值變量，因此我們將為每個家庭添加假人。
train2.info()
train2_dummy = pd.get_dummies(train2, columns=['family'])
train2_dummy
<class 'pandas.core.frame.DataFrame'>
Index: 1931233 entries, 0 to 2061757
Data columns (total 9 columns):
Column Dtype
0 date datetime64[ns]
1 store_nbr int64
2 family object
3 sales float64
4 onpromotion int64
5 codetime int64
6 month int32
7 day int32
8 day_of_week int32
dtypes: datetime64[ns](1), float64(1), int32(3), int64(3), object(1)
memory usage: 125.2+ MB
date store_nbr sales onpromotion codetime month day day_of_week family_AUTOMOTIVE family_BABY CARE ... family_MAGAZINES family_MEATS family_PERSONAL CARE family_PET SUPPLIES family_PLAYERS AND ELECTRONICS family_POULTRY family_PREPARED FOODS family_PRODUCE family_SCHOOL AND OFFICE SUPPLIES family_SEAFOOD
0 2013-01-01 25 2.000 0 0 1 1 1 False False ... False False False False False False False False False False
1 2013-01-01 25 810.000 0 0 1 1 1 False False ... False False False False False False False False False False
2 2013-01-01 25 180.589 0 0 1 1 1 False False ... False False False False False False False False False False
3 2013-01-01 25 186.000 0 0 1 1 1 False False ... False False False False False False False False False False
4 2013-01-01 25 143.000 0 0 1 1 1 False False ... False False False False False False False False False False
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2061752 2017-08-15 9 6.000 0 1687 8 15 1 False False ... False False False False True False False False False False
2061753 2017-08-15 9 438.133 0 1687 8 15 1 False False ... False False False False False True False False False False
2061754 2017-08-15 9 154.553 1 1687 8 15 1 False False ... False False False False False False True False False False
2061755 2017-08-15 9 2419.729 148 1687 8 15 1 False False ... False False False False False False False True False False
2061757 2017-08-15 9 16.000 0 1687 8 15 1 False False ... False False False False False False False False False True
# model training with RandomForest
Xrf = train2_dummy.drop(['sales','date'], axis=1)
yrf = train2_dummy['sales']
# splitting data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(Xrf, yrf, test_size=0.3)
# creating the RandomForest model (hyperparameter temporarily set to: n_estimators=50, max_depth=10, n_jobs=-1, random_state=42)
rfmodel = RandomForestRegressor(n_estimators=50, max_depth=10, n_jobs=-1, random_state=42)
# training the model
rfresult = rfmodel.fit(X_train, y_train)
# test the random forest model for the test part of dataset train2_dummy
y_pred = rfmodel.predict(X_test)
print('The model score is: ',rfmodel.score(X_test,y_test))
The model score is: 0.9150032211616177
# RMSE
#y_pred = rfmodel.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
rmse = mse**.5
print(mse)
print(rmse)
80046.52756956528
282.92495041894995
# RMSLE for train2_dummy
log_actual = np.log1p(y_test)
log_pred = np.log1p(y_pred)
rmsle = np.sqrt(np.mean((log_pred - log_actual) ** 2))
rmsle
print("RMSLE:", rmsle)
RMSLE: 1.4590769436890072
mse = mean_squared_error(y_test, y_pred)
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print("Mean Squared Error:", mse)
print("Mean Absolute Error:", mae)
print("R-squared:", r2)
Mean Squared Error: 80046.52756956528
Mean Absolute Error: 134.8164037890721
R-squared: 0.9150032211616177
模型2的討論
模型2比模型1表現得好得多，我們將評估數據處理、性能度量和殘差圖，以了解為什么它表現得更好。
在查看模型2的調整后的r平方時，它遠高于我們從模型1得到的調整后的r平方值0.9。這意味著90%的銷售差異可以用商店數量和產品系列來解釋。雖然我們無法看到任何p值或系數，但我們能夠計算RMSE，如果我們決定在將來創建它們，可以將其與其他模型進行比較。
隨機森林的殘差圖顯示了一個明確的模式，盡管我們不能確定它是什么。這些模式可能源于其他預測變量，如促銷、假期和石油銷售，我們無法將其包含在本次分析中。
我們無法進一步研究隨機森林模型中顯示的趨勢，但我們可以得出結論，模型2在預測商店銷售方面比模型表現得更好。
?
浙公網安備 33010602011771號