宋 楠，劉際鑫，李 林，仇道霞

（山東省煙草專(zhuān)賣(mài)局（公司），山東 濟(jì)南 250101）

0 引言

1 Holt-Winters模型

通過(guò)分析多個(gè)地市的卷煙銷(xiāo)量數(shù)量發(fā)現(xiàn)，卷煙的銷(xiāo)量存在季節(jié)性的變化規(guī)律，卷煙銷(xiāo)量和卷煙銷(xiāo)量增長(zhǎng)率在時(shí)間序列上存在著一段周期的與下一個(gè)周期的關(guān)聯(lián)。

Holt和Winters將Holt方法進(jìn)行拓展用來(lái)捕獲季節(jié)因素。Holt-Winters季節(jié)性方法包括預(yù)測(cè)方程和三個(gè)平滑方程：一個(gè)用于水平lt，一個(gè)用于趨勢(shì)bt，另一個(gè)用于季節(jié)性分量st，相應(yīng)的平滑參數(shù)分別為2α β和γ。我們用m來(lái)表示季節(jié)頻率，即一年中包含的季節(jié)數(shù)。例如，季度數(shù)據(jù)的m=4，月度數(shù)據(jù)的m=12這種方法有兩種不同的季節(jié)性組成部分。當(dāng)季節(jié)變化在該時(shí)間序列中大致保持不變時(shí)，通常選擇加法模型；而當(dāng)季節(jié)變化與時(shí)間序列的水平成比例變化時(shí)，通常選擇乘法模型。

在加法模型中，季節(jié)性分量在觀(guān)測(cè)序列的尺度上以絕對(duì)值表示，在水平方程中，時(shí)間序列通過(guò)減去季節(jié)分量進(jìn)行季節(jié)性調(diào)整，并且每年的季節(jié)性分量加起來(lái)大約為零。在乘法模型中，季節(jié)性分量用相對(duì)數(shù)（百分比）表示，時(shí)間序列通過(guò)除以季節(jié)性分量來(lái)進(jìn)行季節(jié)性調(diào)整，并且每年的季節(jié)性分量加起來(lái)約為m。

1.1 Holt-Winters加法模型

加法模型的分量形式為：

水平方程表示在 t時(shí)刻，季節(jié)性調(diào)整的觀(guān)察值 （yt - st-m）與非季節(jié)性預(yù)測(cè)值 （ lt-1+ bt-1）之間的加權(quán)平均值。趨勢(shì)方程與Holt線(xiàn)性方法相同。季節(jié)性方程表示當(dāng)前季節(jié)性指數(shù)， （ yt- lt-1- bt-1），和去年同一季節(jié)（即m個(gè)時(shí)間段前）的季節(jié)性指數(shù)之間的加權(quán)平均值。季節(jié)性分量的方程通常表示為：

如果用平滑方程中的lt代替上述分量形式中的水平，可以得到：

2.1 兩組女性血清中AsAb抗體水平檢測(cè)結(jié)果比較 在A(yíng)sAb-T檢測(cè)水平方面，不孕不育組總陽(yáng)性率(65.8%)與健康對(duì)照組總陽(yáng)性率(5.4%)相比明顯提高(P<0.05)；其中不孕不育組女性血清IgG-AsAb、IgM- AsAb、IgA- AsAb的總陽(yáng)性率均明顯高于健康對(duì)照組(均P<0.05)。見(jiàn)表1。

與規(guī)定的季節(jié)性分量的平滑方程相同，其中γ = γ*（ 1 - α ）。通常的參數(shù)限制是 0 ≤ γ*≤ 1,將其轉(zhuǎn)換為0 ≤ γ ≤ 1 -α。

1.2 Holt-Winters乘法模型

乘法模型的分量形式為：

以一個(gè)卷煙單品為例，從日銷(xiāo)量圖可以看出，單品的每日銷(xiāo)量有周期性的規(guī)律，日銷(xiāo)量數(shù)據(jù)也可以保持一段時(shí)間的趨勢(shì)性。holt-winters算法主要考慮的是趨勢(shì)性和周期性，因此可以使用holt-winters算法對(duì)單品的銷(xiāo)量數(shù)據(jù)進(jìn)行推斷。但是我們從日銷(xiāo)量圖中也可以看出，銷(xiāo)量數(shù)據(jù)也會(huì)出現(xiàn)較大的波動(dòng)，這些波動(dòng)很難捕捉到，而且會(huì)對(duì)基于holt-winters的銷(xiāo)量推斷產(chǎn)生影響，由于這些銷(xiāo)量數(shù)據(jù)是真實(shí)的，不是由于數(shù)據(jù)填充錯(cuò)誤導(dǎo)致的，因此后面我們引入prophet算法。

圖1 卷煙銷(xiāo)量波動(dòng)圖Fig.1 Fluctuation chart of cigarette sales

2 prophet模型

prophet算法模型不僅可以處理時(shí)間序列存在一些異常值的情況，也可以處理部分缺失值的情形，還能夠幾乎全自動(dòng)地預(yù)測(cè)時(shí)間序列未來(lái)的走勢(shì)。prophet算法模型是基于時(shí)間序列分解和機(jī)器學(xué)習(xí)的擬合來(lái)做的，其中在擬合模型的時(shí)候使用了pyStan這個(gè)開(kāi)源工具，因此能夠在較快的時(shí)間內(nèi)得到需要預(yù)測(cè)的結(jié)果。使用prophet模型可以構(gòu)建趨勢(shì)項(xiàng)模型、構(gòu)建季節(jié)性趨勢(shì)、對(duì)節(jié)假日效應(yīng)進(jìn)行預(yù)估。

2.1 趨勢(shì)項(xiàng)模型

在Prophet算法里面，趨勢(shì)項(xiàng)有兩個(gè)重要的函數(shù)，一個(gè)是基于邏輯回歸函數(shù)（logistic function）的，另一個(gè)是基于分段線(xiàn)性函數(shù)（piecewise linear function）的。

如果回顧邏輯回歸函數(shù)的話(huà)，一般都會(huì)想起這樣的形式：

式（6）中C稱(chēng)為曲線(xiàn)的最大漸近值，k表示曲線(xiàn)的增長(zhǎng)率，m表示曲線(xiàn)的中點(diǎn)。當(dāng)C=1，k=1，m=0時(shí)，恰好是常見(jiàn)的 sigmoid函數(shù)的形式。從sigmoid的函數(shù)表達(dá)式來(lái)看，它滿(mǎn)足的微分方程為：

那么，使用分離變量法來(lái)求解微分方程y′=y（1-y），可以得到：

在Prophet里面，是需要設(shè)置變點(diǎn)的位置的，而每一段的趨勢(shì)和走勢(shì)也是會(huì)根據(jù)變點(diǎn)的情況而改變的。在程序里面有兩種方法，一種是通過(guò)人工指定的方式指定變點(diǎn)的位置；另外一種是通過(guò)算法來(lái)自動(dòng)選擇。在默認(rèn)的函數(shù)里面，Prophet 會(huì)選擇n_changepoints = 25個(gè)變點(diǎn)，然后設(shè)置變點(diǎn)的范圍是前80%（changepoint_range），也就是在時(shí)間序列的前 80%的區(qū)間內(nèi)會(huì)設(shè)置變點(diǎn)。通過(guò)forecaster.py里面的set_changepoints函數(shù)可以知道，首先要看一些邊界條件是否合理，例如時(shí)間序列的點(diǎn)數(shù)是否少于n_changepoints 等內(nèi)容；其次如果邊界條件符合，那變點(diǎn)的位置就是均勻分布的，這一點(diǎn)可以通過(guò)np.linspace這個(gè)函數(shù)看出來(lái)。

假設(shè)已經(jīng)放置了S個(gè)變點(diǎn)了，并且變點(diǎn)的位置是在時(shí)間戳 sj,1≤j≤s上，那么在這些時(shí)間戳上，需要給出增長(zhǎng)率的變化，也就是在時(shí)間戳sj上發(fā)生的 change in rate。可以假設(shè)有這樣一個(gè)向量： δ ∈Rs，其中δj表示在時(shí)間戳sj上的增長(zhǎng)率的變化量。如果一開(kāi)始的增長(zhǎng)率使用k來(lái)代替，則在時(shí)間戳t上的增長(zhǎng)率是：

以實(shí)際的煙草銷(xiāo)量數(shù)據(jù)為例，通過(guò)銷(xiāo)量數(shù)據(jù)可以得出，銷(xiāo)量數(shù)據(jù)是存在分段線(xiàn)性銷(xiāo)量趨勢(shì)的，因此可以通過(guò)設(shè)置趨勢(shì)來(lái)捕獲銷(xiāo)量的走勢(shì)。

圖2 卷煙銷(xiāo)量趨勢(shì)項(xiàng)模型分析圖Fig.2 Analysis diagram of cigarette sales trend item model

2.2 季節(jié)性趨勢(shì)

幾乎所有的時(shí)間序列預(yù)測(cè)模型都會(huì)考慮這個(gè)因素，因?yàn)闀r(shí)間序列通常會(huì)隨著天，周，月，年等季節(jié)性的變化而呈現(xiàn)季節(jié)性的變化，也稱(chēng)為周期性的變化。卷煙的銷(xiāo)售同樣會(huì)受季節(jié)的影響。對(duì)于周期函數(shù)而言，大家能夠馬上聯(lián)想到的就是正弦余弦函數(shù)。而在數(shù)學(xué)分析中，區(qū)間內(nèi)的周期性函數(shù)是可以通過(guò)正弦和余弦的函數(shù)來(lái)表示的：假設(shè)f（x）是以2π為周期的函數(shù)，那么它的傅立葉級(jí)數(shù)為：

使用傅立葉級(jí)數(shù)來(lái)模擬時(shí)間序列的周期性。假設(shè) P表示時(shí)間序列的周期，p=365.25，表示以年為周期，p=7，表示以周為周期。它的傅立葉級(jí)數(shù)的形式為：

對(duì)于以年為周期的序列（P=365.25）而言，N=10；對(duì)于以周為周期的序列（p=7）而言，N=3。參數(shù)可以形成列向量：

我們以一個(gè)單品為例，拆分prophet模型的趨勢(shì)項(xiàng)與周期項(xiàng)得到如下圖所示，改單品整體的銷(xiāo)量為下降趨勢(shì)，有周期性的銷(xiāo)量波動(dòng)。

圖3 卷煙銷(xiāo)量季節(jié)性趨勢(shì)分析圖Fig.3 Seasonal trend analysis of cigarette sales

2.3 節(jié)假日效應(yīng)

通過(guò)對(duì)時(shí)間序列的深度分析，我們發(fā)現(xiàn)卷煙的銷(xiāo)量除了季節(jié)性的變化趨勢(shì)外，還會(huì)受節(jié)假日的影響，如春節(jié)、中秋等節(jié)假日會(huì)使卷煙的銷(xiāo)量發(fā)生比較大的波動(dòng)，所以我們又結(jié)合節(jié)假日的特征對(duì)模型進(jìn)行了進(jìn)一步訓(xùn)練。

在 Prophet算法模型里面，通過(guò)對(duì)節(jié)假日的特征分析，訓(xùn)練得到節(jié)假日對(duì)卷煙銷(xiāo)售的影響和關(guān)系。由于每個(gè)節(jié)假日對(duì)時(shí)間序列的影響程度不同，例如春節(jié)、國(guó)慶節(jié)是七天的假期，勞動(dòng)節(jié)等假期則較短。因此，不同的節(jié)假日可以看成相互獨(dú)立的模型，并且可以為不同的節(jié)假日設(shè)置不同的前后窗口值，表示該節(jié)假日會(huì)影響前后一段時(shí)間的時(shí)間序列。用數(shù)學(xué)語(yǔ)言來(lái)說(shuō)，對(duì)于第 i個(gè)節(jié)假日來(lái)說(shuō)，Di表示該節(jié)假日的前后一段時(shí)間。為了表示節(jié)假日效應(yīng)，需要一個(gè)相應(yīng)的指示函數(shù)（indicator function），同時(shí)需要一個(gè)參數(shù)ki來(lái)表示節(jié)假日的影響范圍。假設(shè)有L個(gè)節(jié)假日，那么：

2.4 模型擬合

按照以上的解釋?zhuān)瑫r(shí)間序列已經(jīng)可以通過(guò)增長(zhǎng)項(xiàng)、季節(jié)項(xiàng)、節(jié)假日項(xiàng)來(lái)構(gòu)建：

下一步則需要擬合函數(shù)，以下以30天卷煙銷(xiāo)量預(yù)測(cè)為例。

根據(jù)卷煙日銷(xiāo)量數(shù)據(jù)匯聚成 30天的銷(xiāo)量數(shù)據(jù)，當(dāng)前的銷(xiāo)量為30天的歷史銷(xiāo)量的總和，使用prophet算法的預(yù)測(cè)共得到三條銷(xiāo)量的曲線(xiàn)，最上面一條為銷(xiāo)量預(yù)測(cè)的上界即最大的銷(xiāo)量，中間一條曲線(xiàn)為預(yù)測(cè)的銷(xiāo)量，最下面一條曲線(xiàn)為銷(xiāo)量預(yù)測(cè)的下界，散點(diǎn)為實(shí)際的銷(xiāo)量，通過(guò)圖4可以看出，雖然30天總的銷(xiāo)量相比于日銷(xiāo)量來(lái)說(shuō)數(shù)據(jù)的分布更穩(wěn)定，但是依然會(huì)數(shù)據(jù)波動(dòng)非常大，而且只用銷(xiāo)量的數(shù)據(jù)，算法在預(yù)測(cè)結(jié)果的穩(wěn)定性方面沒(méi)有考慮更多的因素，因此引入機(jī)器學(xué)習(xí)與時(shí)間序列模型融合的方法，在機(jī)器學(xué)習(xí)特征設(shè)計(jì)的時(shí)候使用了時(shí)間序列預(yù)測(cè)的特征，同時(shí)加入了非時(shí)間序列的特征。

圖4 卷煙銷(xiāo)量模型擬合分析圖Fig.4 Fitting analysis diagram of cigarette sales model

3 xgboost模型

XGBoost是 boosting算法的其中一種。Boosting算法的思想是將許多弱分類(lèi)器集成在一起形成一個(gè)強(qiáng)分類(lèi)器。因?yàn)閄GBoost是一種提升樹(shù)模型，所以它是將許多樹(shù)模型集成在一起，形成一個(gè)很強(qiáng)的分類(lèi)器，而所用到的樹(shù)模型則是CART回歸樹(shù)模型。該算法思想就是不斷地添加樹(shù)，不斷地進(jìn)行特征分裂來(lái)生長(zhǎng)一棵樹(shù)，每次添加一個(gè)樹(shù)，其實(shí)是學(xué)習(xí)一個(gè)新函數(shù)，去擬合上次預(yù)測(cè)的殘差。當(dāng)訓(xùn)練完成得到k棵樹(shù)，則需預(yù)測(cè)一個(gè)樣本的分?jǐn)?shù)，根據(jù)樣本的特征，在每棵樹(shù)中會(huì)落到對(duì)應(yīng)的一個(gè)葉子節(jié)點(diǎn)，每個(gè)葉子節(jié)點(diǎn)就對(duì)應(yīng)一個(gè)分?jǐn)?shù)，最后只需要將每棵樹(shù)對(duì)應(yīng)的分?jǐn)?shù)加起來(lái)就是該樣本的預(yù)測(cè)值：

式（21）中，wq（x）為葉子節(jié)點(diǎn)的分?jǐn)?shù)，f（x）為其中一棵回歸樹(shù)。

4 銷(xiāo)量預(yù)測(cè)具體流程

由于時(shí)間序列預(yù)測(cè)銷(xiāo)量只能使用歷史的銷(xiāo)量數(shù)據(jù)，不能使用其它數(shù)據(jù)，因此我們采用時(shí)間序列模型與機(jī)器學(xué)習(xí)模型融合的方法，模型的設(shè)計(jì)既使用來(lái)時(shí)間序列銷(xiāo)量特征，又考慮了用戶(hù)特征等非時(shí)間序列銷(xiāo)量特征。prophet和 holt-winters時(shí)間序列方法作為 xgboost預(yù)測(cè)的高級(jí)特征,目標(biāo)預(yù)測(cè)不是直接預(yù)測(cè)未來(lái)的銷(xiāo)量，而是預(yù)測(cè)未來(lái)的增長(zhǎng)率,這樣目標(biāo)變量的波動(dòng)范圍會(huì)小，而且更穩(wěn)定。通過(guò)增長(zhǎng)率與當(dāng)前30天歷史的銷(xiāo)量就可以計(jì)算出未來(lái)的 30天銷(xiāo)量，7天銷(xiāo)量的方法是一致的，只是特征構(gòu)建不同。

圖5 卷煙銷(xiāo)量預(yù)測(cè)流程圖Fig.5 Flow chart of cigarette sales forecast

4.1 30天銷(xiāo)量預(yù)測(cè)的特征

表1 30天銷(xiāo)量預(yù)測(cè)的特征表Tab.1 Feature table for 30-day sales forecast

4.2 7天銷(xiāo)量的預(yù)測(cè)特征

表2 7天銷(xiāo)量的預(yù)測(cè)特征表Tab.2 7-day sales forecast feature table

4.3 模型驗(yàn)證

使用2017年1月1日到2020年7月1日，某省三個(gè)地市的六種規(guī)格的卷煙銷(xiāo)量數(shù)據(jù)，根據(jù)MAPE進(jìn)行驗(yàn)證銷(xiāo)量預(yù)測(cè)的準(zhǔn)確率，結(jié)果如表3。

表3 卷煙銷(xiāo)量預(yù)測(cè)表Tab.3 Cigarette sales forecast table

6901028315012規(guī)格編號(hào)的卷煙在使用 2020年5月29號(hào)之前的歷史數(shù)據(jù)預(yù)測(cè)2020年5月30號(hào)到6月28號(hào)周期30天的銷(xiāo)量數(shù)據(jù)時(shí)，從銷(xiāo)量數(shù)據(jù)來(lái)看，5月30號(hào)到6月28號(hào)的銷(xiāo)量數(shù)據(jù)忽然下降很多，導(dǎo)致預(yù)測(cè)結(jié)果有一定的偏差， 與銷(xiāo)售部門(mén)確認(rèn)其原因是貨源不足等因素干擾。

圖6 卷煙銷(xiāo)量異常波動(dòng)圖Fig.6 Abnormal fluctuations in cigarette sales

5 結(jié)論

為科學(xué)的制定市級(jí)煙草商業(yè)企業(yè)的卷煙投放策略，采用基于prophet、Holt-Winters和xgboost模型構(gòu)成的混合模型，對(duì)基于地市鄉(xiāng)鎮(zhèn)為單位的銷(xiāo)量前十的單品卷煙銷(xiāo)量進(jìn)行預(yù)測(cè)，提高了單品卷煙預(yù)測(cè)的精確度，在卷煙單品銷(xiāo)量預(yù)測(cè)中取得了較好的結(jié)果，但文中的混合模型僅以歷史銷(xiāo)售數(shù)據(jù)為輸入，相關(guān)結(jié)論在數(shù)據(jù)存在政策因素影響時(shí)并不適用。