一、需求背景
我們福祿網路致力於為廣大使用者提供智慧化充值服務，包括各類通訊充值卡（比如移動、聯通、電信的話費及流量充值）、遊戲類充值卡（比如王者榮耀、吃雞類點券、AppleStore充值、Q幣、鬥魚幣等）、生活服務類（比如肯德基、小鹿茶等），網娛類（比如QQ各類鑽等），作為一個服務提供商，商品質量的穩定、持續及充值過程的便捷一直是我們在業內的口碑。
在整個商品流通過程中，如何做好庫存的管理，以充分提高庫存運轉週期和資金使用效率，一直是個難題。基於此，我們提出了智慧化的庫存管理服務，根據訂單資料及商品資料，來預測不同商品隨著時間推移的日常消耗情況。
二、演算法選擇
目前成熟的時間序列預測演算法很多，但商業領域效能優越的卻不多，經過多種嘗試，給大家推薦2種時間序列演算法：facebook開源的Prophet演算法和LSTM深度學習演算法。
現將個人理解的2種演算法特性予以簡要說明：
- (1)、在訓練時間上，prophet幾十秒就能出結果，而lstm往往需要1個半小時，更是隨著網路層數和特徵數量的增加而增加。
- (2)、Prophet是一個為商業預測而生的時間序列預測模型，因此在很多方便都有針對性的優化，而lstm的初衷是nlp。
- (3)、Prophet無需特徵處理即可使用，引數調優也明確簡單。而lstm則需要先進行必要的特徵處理，其次要進行正確的網路結構設計，因此lstm相對prophet更為複雜。
- (4)、Lstm需要更多的資料進行學習，否則無法消除欠擬合的情形。而prophet不同，prophet基於統計學，有完整的數學理論支撐，因此更容易從少量的資料中完成學習。
- (5)、傳統的時間序列預測演算法只支援單緯度，但LSTM能支援多緯度，也就是說LSTM能考慮促銷活動，目標使用者特性，產品特性等
三、資料來源
- (1)、訂單資料
- (2)、產品分類資料
四、資料形式
time,product,cnt
2019-10-01 00,**充值,6
2019-10-01 00,***遊戲,368
2019-10-01 00,***,1
2019-10-01 00,***,11
2019-10-01 00,***遊戲,17
2019-10-01 00
,三網***,39
2019-10-01 00,**網,6
2019-10-01 00,***,2
欄位說明：
- Time：小時級時間
- Product：產品名稱或產品的分類名稱，目前使用的是產品2級分類，名稱
- Cnt：成功訂單數量
目前的時間序列是由以上time和cnt組成，product是用於區分不同時間序列的欄位。
五、特徵處理
時間序列一般不進行特徵處理，當然可以根據具體情況進行歸一化處理或是取對數處理等。
六、演算法選擇
目前待選的演算法主要有2種：
- (1)、Prophet
Facebook開源的時間序列預測演算法，考慮了節假日因素。 - (2)、LSTM
優化後的RNN深度學習演算法。
七、演算法說明
7.1 prophet
7.1.1Prophet的核心是調參，步驟如下：
- 1、首先我們去除資料中的異常點（outlier），直接賦值為none就可以，因為Prophet的設計中可以通過插值處理缺失值，但是對異常值比較敏感。
- 2、選擇趨勢模型，預設使用分段線性的趨勢，但是如果認為模型的趨勢是按照log函式方式增長的，可設定growth='logistic'從而使用分段log的增長方式
- 3、 設定趨勢轉折點（changepoint），如果我們知道時間序列的趨勢會在某些位置發現轉變，可以進行人工設定，比如某一天有新產品上線會影響我們的走勢，我們可以將這個時刻設定為轉折點。如果自己不設定，演算法會自己總結changepoint。
- 4、 設定週期性，模型預設是帶有年和星期以及天的週期性，其他月、小時的週期性需要自己根據資料的特徵進行設定，或者設定將年和星期等週期關閉。
設定節假日特徵，如果我們的資料存在節假日的突增或者突降，我們可以設定holiday引數來進行調節，可以設定不同的holiday，例如五一一種，國慶一種，影響大小不一樣，時間段也不一樣。 - 5、 此時可以簡單的進行作圖觀察，然後可以根據經驗繼續調節上述模型引數，同時根據模型是否過擬合以及對什麼成分過擬合，我們可以對應調節seasonality_prior_scale、holidays_prior_scale、changepoint_prior_scale引數。
以上是理論上的調參步驟，但我們在實際情況下在建議使用grid_search（網格尋參）方式，直接簡單效果好。當機器效能不佳時網格調參配合理論調參方法可以加快調參速度。建議初學者使用手動調參方式以理解每個引數對模型效果的影響。
holiday.csv
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from fbprophet import Prophet
data = pd.read_csv('../data/data2.csv', parse_dates=['time'], index_col='time')
def get_product_data(name, rule=None):
product = data[data['product'] == name][['cnt']]
product.plot()
if rule is not None:
product = product.resample(rule).sum()
product.reset_index(inplace=True)
product.columns = ['ds', 'y']
return product
holidays = pd.read_csv('holiday.csv', parse_dates=['ds'])
holidays['lower_window'] = -1
holidays = holidays.append(pd.DataFrame({
'holiday': '雙11',
'ds': pd.to_datetime(['2019-11-11', '2020-11-11']),
'lower_window': -1,
'upper_window': 1,
})).append(pd.DataFrame({
'holiday': '雙12',
'ds': pd.to_datetime(['2019-12-12', '2020-12-12']),
'lower_window': -1,
'upper_window': 1,
})
)
def predict(name, rule='1d', freq='d', periods=1, show=False):
ds = get_product_data(name, rule=rule)
if ds.shape[0] < 7:
return None
m = Prophet(holidays=holidays)
m.fit(ds)
future = m.make_future_dataframe(freq=freq, periods=periods) # 建立資料預測框架，資料粒度為天，預測步長為一年
forecast = m.predict(future)
if show:
m.plot(forecast).show() # 繪製預測效果圖
m.plot_components(forecast).show() # 繪製成分趨勢圖
mse = forecast['yhat'].iloc[ds.shape[0]] - ds['y'].values
mse = np.abs(mse) / (ds['y'].values + 1)
return [name, mse.mean(), mse.max(), mse.min(), np.quantile(mse, 0.9), np.quantile(mse, 0.8), mse[-7:].mean(),
ds['y'].iloc[-7:].mean()]
if __name__ == '__main__':
products = set(data['product'])
p = []
for i in products:
y = predict(i)
if y is not None:
p.append(y)
df = pd.DataFrame(p, columns=['product', 'total_mean', 'total_max', 'total_min', '0.9', '0.8', '7_mean',
'7_real_value_mean'])
df.set_index('product', inplace=True)
product_sum: pd.DataFrame = data.groupby('product').sum()
df = df.join(product_sum)
df.sort_values('cnt', ascending=False, inplace=True)
df.to_csv('result.csv', index=False)
結果如下：由於行數較多這裡只展示前1行
根據結果，對比原生資料，可以得出如下結論：
就演算法與產品的匹配性可分為3個型別：
- （1）與演算法較為匹配，演算法的歷史誤差8分為數<=0.2的
- （2）與演算法不太匹配的，演算法的歷史誤差8分為數>0.2的
- （3）資料過少的，無法正常預測的。目前僅top10就能佔到整體訂單數的90%以上。
7.1.2 部分成果展示
A. 因素分解圖
上圖中主要分為3個部分，分別對應prophet 3大要素，趨勢、節假日或特殊日期、週期性（包括年週期、月週期、week週期、天週期以及使用者自定義的週期）
下面依照上面因素分解圖的順序依次對圖進行說明：
- （1)、Trend:
即趨勢因素圖。描述時間序列的趨勢。Prophet支援線性趨勢和logist趨勢。通過growth引數設定，當然模型能自己根據時間序列的走勢判斷growth型別。這也是prophet實現的比較智慧的一點。 - （2)、Holidays
即節假日及特殊日期因素圖。描述了節假日及使用者自定義的特殊日期對時間序列的影響。正值為正影響，負值為負影響。從圖中可以看出這個商品對節假日比較敏感。節假日是根據holidays引數設定的。 - （3)、weekly
星期週期性因素圖。正常情況下，如果是小時級別資料將會有天週期圖。有1年以上完整資料並且時間序列有典型的年週期性會有年週期圖。如果你覺得這個有年週期，但模型並不這麼認為，你可以通過設定yearly_seasonality設定一個具體的數值。這個數值預設情況下為10（weekly_seasonality預設為3），這個值代表的是傅立葉級數的項數，越大模型越容易過擬合，過小則會導致欠擬合，一般配合seasonality_prior_scale使用。
B.預測曲線與實際值對比
7.2 lstm
LSTM（長短記憶網路）主要用於有先後順序的序列型別的資料的深度學習網路。是RNN的優化版本。一般用於自然語言處理，也可用於時間序列的預測。
簡單來說就是，LSTM一共有三個門，輸入門，遺忘門，輸出門， i 、o、 f 分別為三個門的程度引數， g 是對輸入的常規RNN操作。公式裡可以看到LSTM的輸出有兩個，細胞狀態c 和隱狀態 h，c是經輸入、遺忘門的產物，也就是當前cell本身的內容，經過輸出門得到h，就是想輸出什麼內容給下一單元。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import torch
from torch import nn
from sklearn.preprocessing import MinMaxScaler
ts_data = pd.read_csv('../data/data2.csv', parse_dates=['time'], index_col='time')
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
n_vars = 1 if type(data) is list else data.shape[1]
df = pd.DataFrame(data)
cols, names = list(), list()
# input sequence (t-n, ... t-1)
for i in range(n_in, 0, -1):
cols.append(df.shift(i))
names += [('var%d(t-%d)' % (j + 1, i)) for j in range(n_vars)]
# forecast sequence (t, t+1, ... t+n)
for i in range(0, n_out):
cols.append(df.shift(-i))
if i == 0:
names += [('var%d(t)' % (j + 1)) for j in range(n_vars)]
else:
names += [('var%d(t+%d)' % (j + 1, i)) for j in range(n_vars)]
# put it all together
agg = pd.concat(cols, axis=1)
agg.columns = names
# drop rows with NaN values
if dropnan:
agg.dropna(inplace=True)
return agg
def transform_data(feature_cnt=2):
yd = ts_data[ts_data['product'] == '移動話費'][['cnt']]
scaler = MinMaxScaler(feature_range=(0, 1))
yd_scaled = scaler.fit_transform(yd.values)
yd_renamed = series_to_supervised(yd_scaled
, n_in=feature_cnt).values.astype('float32')
n_row = yd_renamed.shape[0]
n_train = int(n_row * 0.7)
train_X, train_y = yd_renamed[:n_train, :-1], yd_renamed[:n_train, -1]
test_X, test_y = yd_renamed[n_train:, :-1], yd_renamed[n_train:, -1]
# 最後，我們需要將資料改變一下形狀，因為 RNN 讀入的資料維度是 (seq, batch, feature)，所以要重新改變一下資料的維度，這裡只有一個序列，所以 batch 是 1，而輸入的 feature 就是我們希望依據的幾天，這裡我們定的是兩個天，所以 feature 就是 2.
train_X = train_X.reshape((-1, 1, feature_cnt))
test_X = test_X.reshape((-1, 1, feature_cnt))
print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)
# 轉化成torch 的張量
train_x = torch.from_numpy(train_X)
train_y = torch.from_numpy(train_y)
test_x = torch.from_numpy(test_X)
test_y = torch.from_numpy(test_y)
return scaler, train_x, train_y, test_x, test_y
scaler, train_x, train_y, test_x, test_y = transform_data(24)
# lstm 網路
class lstm_reg(nn.Module): # 括號中的是python的類繼承語法，父類是nn.Module類 不是引數的意思
def __init__(self, input_size, hidden_size, output_size=1, num_layers=2): # 建構函式
# inpu_size 是輸入的樣本的特徵維度， hidden_size 是LSTM層的神經元個數，
# output_size是輸出的特徵維度
super(lstm_reg, self).__init__() # super用於多層繼承使用，必須要有的操作
self.rnn = nn.LSTM(input_size, hidden_size, num_layers) # 兩層LSTM網路，
self.reg = nn.Linear(hidden_size, output_size) # 把上一層總共hidden_size個的神經元的輸出向量作為輸入向量，然後迴歸到output_size維度的輸出向量中
def forward(self, x): # x是輸入的資料
x, _ = self.rnn(x) # 單個下劃線表示不在意的變數，這裡是LSTM網路輸出的兩個隱藏層狀態
s, b, h = x.shape
x = x.view(s * b, h)
x = self.reg(x)
x = x.view(s, b, -1) # 使用-1表示第三個維度自動根據原來的shape 和已經定了的s,b來確定
return x
def train(feature_cnt, hidden_size, round, save_path='model.pkl'):
# 我使用了GPU加速，如果不用的話需要把.cuda()給註釋掉
net = lstm_reg(feature_cnt, hidden_size)
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(net.parameters(), lr=1e-2)
for e in range(round):
# 新版本中可以不使用Variable了
# var_x = Variable(train_x).cuda()
# var_y = Variable(train_y).cuda()
# 將tensor放在GPU上面進行運算
var_x = train_x
var_y = train_y
out = net(var_x)
loss = criterion(out, var_y)
optimizer.zero_grad()
loss.backward()
optimizer.step()
if (e + 1) % 100 == 0:
print('Epoch: {}, Loss:{:.5f}'.format(e + 1, loss.item()))
# 儲存訓練好的模型引數
torch.save(net.state_dict(), save_path)
return net
if __name__ == '__main__':
net = train(24, 8, 5000)
# criterion = nn.MSELoss()
# optimizer = torch.optim.Adam(net.parameters(), lr=1e-2)
pred_test = net(test_x) # 測試集的預測結果
pred_test = pred_test.view(-1).data.numpy() # 先轉移到cpu上才能轉換為numpy
# 乘以原來歸一化的刻度放縮回到原來的值域
origin_test_Y = scaler.inverse_transform(test_y.reshape((-1,1)))
origin_pred_test = scaler.inverse_transform(pred_test.reshape((-1,1)))
# 畫圖
plt.plot(origin_pred_test, 'r', label='prediction')
plt.plot(origin_test_Y, 'b', label='real')
plt.legend(loc='best')
plt.show()
# 計算MSE
# loss = criterion(out, var_y)？
true_data = origin_test_Y
true_data = np.array(true_data)
true_data = np.squeeze(true_data) # 從二維變成一維
MSE = true_data - origin_pred_test
MSE = MSE * MSE
MSE_loss = sum(MSE) / len(MSE)
print(MSE_loss)
八、兩種演算法的比較
- (1)在訓練時間上，prophet幾十秒就能出結果，而lstm往往需要1個半小時，更是隨著網路層數和特徵數量的增加而增加。
- (2)Prophet是一個為商業預測而生的時間序列預測模型，因此在很多方便都有針對性的優化，而lstm的初衷是nlp。
- (3)Prophet無需特徵處理即可使用，引數調優也明確簡單。而lstm則需要先進行必要的特徵處理，其次要進行正確的網路結構設計，因此lstm相對prophet更為複雜。
- (4)Lstm需要更多的資料進行學習，否則無法消除欠擬合的情形。而prophet不同，prophet基於統計學，有完整的數學理論支撐，因此更容易從少量的資料中完成學習。
參考文獻：
【1】Prophet官方文件：https://facebook.github.io/prophet/
【2】Prophet論文：https://peerj.com/preprints/3190/
【3】Prophet-github:https://github.com/facebook/prophet
【4】LSTM http://colah.github.io/posts/2015-08-Understanding-LSTMs/
【5】基於LSTM的關聯時間序列預測方法研究 尹康 《北京交通大學》 2019年 cnki地址：http://cdmd.cnki.com.cn/Article/CDMD-10004-1019209125.htm
