前言

之前，我們都在影像、語言等基礎應用上打轉，這次我們要來探討一個可應用在企業運作上的實例，銷售預測主要是希望藉由過去的銷售記錄預測下一個週期的銷售量，在統計上，我們會使用簡單迴歸，乃至複雜的『時間序列分析』(Time Series Analysis)來預測銷售趨勢，因為，當期的銷售量通常會與前期的銷售量有緊密的關係，除非公司發生重大事件，否則，應該會循著規律變化。還記得嗎? 在『自然語言處理』時，我們會使用LSTM考慮上下文的關係，這個模型恰好與前面講的銷售量預測不謀而合，所以，本篇就以 LSTM 模型來預測銷售量。

銷售量預測的樣態很多種，包括營收、利潤、來客數、遊園人數、銷售產品數/金額、...等等，都屬於同一範疇，本篇會以航空公司的每月乘客人數為例，使用 LSTM 模型預測下個月的乘客數。

時間序列(Time Series Analysis)概念淺介

簡單迴歸(Regression) 公式 y=ax+b，是基於 y(i) 與 y(j) 是相互獨立，沒有任何關聯，但在銷售量的表現上，這個假設並不合理，公司銷售業績通常不會暴漲暴跌，而是『逐步』上升或下跌，也就是與前期的表現有緊密的關聯，另外，大部分的公司也會有淡、旺季，即所謂的『季節效應』(Seasonal Effect)，因此，使用更複雜的『時間序列分析』(Time Series Analysis)預測會更貼近事實，時間序列分析的模型因應問題的型態不同也有很多種，我們以ARIMA(Autoregressive Integrated Moving Average)為例，做簡單的說明，如果你對以下的說明很頭痛，可跳過，直接看 Neural Network 怎麼作。

參考Wiki的簡介，ARIMA 是 ARMA 的擴充模型，而 ARMA 就等於 AR + MA，AR 談的就是與前期的關係，有一重要的參數 p，代表期數，就是假設當期與前面p期有關係，MA 是移動平均數(Moving Average)就是避免異常值影響過大，故利用前幾期的平均數來代表當期量，所以，也有一個重要的參數 q，代表計算移動平均的期數，另外，ARIMA 多一個參數 d，稱為差分階數，是要使時間序列維持『穩態』(stationarity)，所以，總共有三個參數 p、q、d，相關的理論是一門學科，筆者只能點到為止。

至於p、q如何決定是多少呢? 可利用『自我相關函數』(Autocorrelation Function，ACF)或『偏自我相關函數』(Partial Autocorrelation Function, PACF)來協助判斷p及q值，至於 d 值的判斷可參考『Identifying the order of differencing in an ARIMA model』，經過反覆的數據觀察、分析與實驗，才可以得到較接近事實的預測模型。

綜合以上說明，我們可以了解以統計模型預測銷售量，是須具備一定的知識以及行業別知識(Domain Know How)，才能建立較正確的預測模型，接下來，我們來看看 Neural Network 的方法。

相關的實作可使用 StatsModels 套件，它支援『時間序列分析』與繪圖。

Neural Network 作法

以下說明及實作主要參考『Time Series Prediction with LSTM Recurrent Neural Networks in Python with Keras』，資料集採用『航空公司每月乘客人數』，可至這裡下載，我們利用以下程式畫出乘客人數的折線圖：

import pandas
import matplotlib.pyplot as plt
dataset = pandas.read_csv('international-airline-passengers.csv', usecols=[1], engine='python', skipfooter=3)
plt.plot(dataset)
plt.show()

圖. 航空公司每月乘客人數

我們可以觀察到數據有一定的趨勢與波動性，暑假(7~9月)乘客人數比較多，表示有季節效應，我們可以用ACF或PACF確認一下：

import numpy as np
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.api import qqplot
# 畫出 ACF 12 期的效應
sm.graphics.tsa.plot_acf(dataset, lags=12)
plt.show()
# 畫出 PACF 12 期的效應
sm.graphics.tsa.plot_pacf(dataset, lags=12)
plt.show()

** 注意，StatsModels 套件，必須以以下指令安裝，用 pip install 會有錯 **
conda install -c conda-forge statsmodels

圖. ACF 12 期的效應

圖. PACF 12 期的效應

接著，我們就以 LSTM 模型實作，程式碼來自『Time Series Prediction with LSTM Recurrent Neural Networks in Python with Keras』，我加了一些註解，也可至這裡下載，範例在 TimeSeries 資料夾，如下：

# LSTM for international airline passengers problem with regression framing
import numpy
import matplotlib.pyplot as plt
from pandas import read_csv
import math
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error


# 產生 (X, Y) 資料集, Y 是下一期的乘客數
def create_dataset(dataset, look_back=1):
	dataX, dataY = [], []
	for i in range(len(dataset)-look_back-1):
		a = dataset[i:(i+look_back), 0]
		dataX.append(a)
		dataY.append(dataset[i + look_back, 0])
	return numpy.array(dataX), numpy.array(dataY)

# 載入訓練資料
dataframe = read_csv('international-airline-passengers.csv', usecols=[1], engine='python', skipfooter=3)
dataset = dataframe.values
dataset = dataset.astype('float32')
# 正規化(normalize) 資料，使資料值介於[0, 1]
scaler = MinMaxScaler(feature_range=(0, 1))
dataset = scaler.fit_transform(dataset)

# 2/3 資料為訓練資料， 1/3 資料為測試資料
train_size = int(len(dataset) * 0.67)
test_size = len(dataset) - train_size
train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:]

# 產生 (X, Y) 資料集, Y 是下一期的乘客數(reshape into X=t and Y=t+1)
look_back = 1
trainX, trainY = create_dataset(train, look_back)
testX, testY = create_dataset(test, look_back)
# reshape input to be [samples, time steps, features]
trainX = numpy.reshape(trainX, (trainX.shape[0], 1, trainX.shape[1]))
testX = numpy.reshape(testX, (testX.shape[0], 1, testX.shape[1]))

# 建立及訓練 LSTM 模型
model = Sequential()
model.add(LSTM(4, input_shape=(1, look_back)))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
model.fit(trainX, trainY, epochs=100, batch_size=1, verbose=2)

# 預測
trainPredict = model.predict(trainX)
testPredict = model.predict(testX)

# 回復預測資料值為原始數據的規模
trainPredict = scaler.inverse_transform(trainPredict)
trainY = scaler.inverse_transform([trainY])
testPredict = scaler.inverse_transform(testPredict)
testY = scaler.inverse_transform([testY])

# calculate 均方根誤差(root mean squared error)
trainScore = math.sqrt(mean_squared_error(trainY[0], trainPredict[:,0]))
print('Train Score: %.2f RMSE' % (trainScore))
testScore = math.sqrt(mean_squared_error(testY[0], testPredict[:,0]))
print('Test Score: %.2f RMSE' % (testScore))

# 畫訓練資料趨勢圖
# shift train predictions for plotting
trainPredictPlot = numpy.empty_like(dataset)
trainPredictPlot[:, :] = numpy.nan
trainPredictPlot[look_back:len(trainPredict)+look_back, :] = trainPredict

# 畫測試資料趨勢圖
# shift test predictions for plotting
testPredictPlot = numpy.empty_like(dataset)
testPredictPlot[:, :] = numpy.nan
testPredictPlot[len(trainPredict)+(look_back*2)+1:len(dataset)-1, :] = testPredict

# 畫原始資料趨勢圖
# plot baseline and predictions
plt.plot(scaler.inverse_transform(dataset))
plt.plot(trainPredictPlot)
plt.plot(testPredictPlot)
plt.show()

程式執行

將下載的資料檔international-airline-passengers.csv與程式SimpleLSTM.py放在同一目錄，在DOS內執行以下指令：
python SimpleLSTM.py

執行結果如下：

圖. LSTM 預測結果，藍線為實際值，綠線為訓練結果，紅色為測試結果

程式說明

程式不複雜，可參閱中文註解，處理流程如下：

1. 資料檔只有兩欄，日期及乘客數，正規化(normalize) 資料，使資料值介於[0, 1]。
2. 我們要用前期預測當期，故將資料轉為(前期乘客數， 當期乘客數)，當作(X, Y)。
3. 建立及訓練 LSTM 模型，模型很簡單，就只有一個LSTM層及output層。
4. 訓練 LSTM 模型並進行預測。
5. 針對實際值、預測值進行繪圖。

可以看到很神奇的現象，我們沒有特別處理季節效應，但 LSTM 經過 100 輪的訓練後，幫我們考慮到了，模型只使用前期資料為X變數。另外，我們也可以使用多期的資料作為X變數，或者，針對季節效應對 LSTM 作『記憶重置』(reset state)、以期數(Time Steps)當input處理、...等等的方法，針對問題特性，多實驗看看，也許會有意想不到的效果。

結論

其實，我們周遭不乏其他的例子可以實驗，例如雨量、氣候、上市公司營收、股價預測、房價預測...等等，都屬於時間序列，尤其，政府目前提供許多 Open Data，取得資料並不困難，重要的是，我們怎麼看待資料，是否有創新的想法吧。