### 一、需求背景 我們福祿網路致力於為廣大使用者提供智慧化充值服務，包括各類通訊充值卡（比如移動、聯通、電信的話費及流量充值）、遊戲類充值卡（比如王者榮耀、吃雞類點券、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的關聯時間序列預測方法研究 [尹康](http://search.cnki.com.cn/Search/Result?author=%E5%B0%B9%E5%BA%B7) 《北京交通大學》 2019年 cnki地址：http://cdmd.cnki.com.cn/Article/CDMD-10004-1019209125.htm