張家晨，左興權，黃 海，韓 靜，張百勝

1.北京郵電大學 計算機學院（國家示范性軟件學院），北京 100876

2.中興通信股份有限公司，上海 201203

隨著5G網(wǎng)絡以及移動端設備的逐漸普及，人們對于移動數(shù)據(jù)流量的需求不斷增加。隨之而來的數(shù)據(jù)流量將會以指數(shù)形式增長，對于網(wǎng)絡運營商而言，提供不受阻礙、無處不在的高質(zhì)量服務意義重大。其中，如何建立準確的長周期預測基站小區(qū)網(wǎng)絡流量的模型，以指導運營商擴容基站小區(qū)無線設備將是一個挑戰(zhàn)。

基站小區(qū)網(wǎng)絡流量通常以時間序列的形式呈現(xiàn)，預測周期可分為短周期時間序列預測（預測時長按分鐘、小時和天）和長周期時間序列預測（預測時長按天、周和月）。短周期的基站小區(qū)網(wǎng)絡流量預測相對成熟，主要包括三類方法，即傳統(tǒng)統(tǒng)計模型建模方法、神經(jīng)網(wǎng)絡頇測方法、基于序列分解結合以上兩種方法的混合模型方法。

傳統(tǒng)統(tǒng)計模型建模方法主要包括自回歸模型（autoregressive model，AR）[1]、滑動平均模型（moving average model，MA）[2]、回歸滑動平均模型（autoregressive moving average model，ARMA）[3]、差分整合移動平均自回歸模型（autoregressive integrated moving average model，ARIMA）[4]等模型及其改進。文獻[5]使用線性回歸和非線性回歸方法，對4G基站流量進行預測，實驗結果表明非線性預測算法性能優(yōu)于線性機器學習算法對基站流量進行直接建模。文獻[6]使用自回歸滑動平均模型與三指數(shù)平滑模型，預測上傳和下載熱點的流量，結果表明三指數(shù)平滑模型具有更好的預測性能。文獻[7]提出了基于SARIMA（seasonal auto regressive integrated moving average）模型對網(wǎng)絡流量進行預測，實驗表明當預測步長小于10時，預測精度更高。

用于網(wǎng)絡流量預測的神經(jīng)網(wǎng)絡模型包括循環(huán)神經(jīng)網(wǎng)絡（recurrent neural network，RNN）[8]、BP神經(jīng)網(wǎng)絡（back propagation neural network，BPNN）[9]及其變體[10]。文獻[11]利用RNN為網(wǎng)絡流量建立預測模型，可準確預測各種網(wǎng)絡應用場景下的流量變化趨勢。文獻[12]利用螢火蟲群算法優(yōu)化BP神經(jīng)網(wǎng)絡，對網(wǎng)絡流量進行預測，仿真結果表明該算法具有較高的預測精度。

當前，以“分解和集成”方法為基礎的混合模型預測方法在網(wǎng)絡流量預測中表現(xiàn)良好。文獻[13]提出基于EEMD與RBF神經(jīng)網(wǎng)絡組合模型對網(wǎng)絡流量進行預測。文獻[14]為了解決自相似網(wǎng)絡流量預測，提出了集合EMD和ARMA的組合模型。文獻[15]基于小波分解聯(lián)合使用長短期記憶網(wǎng)絡（long short-term memory，LSTM）模型，該組合預測模型在預測流量突發(fā)性方面的性能優(yōu)于普通的LSTM網(wǎng)絡模型。文獻[16]提出了一種基于小波變換和SARIMA模型的方法，該方法比單一預測模型具有更高的預測精度。文獻[17]提出基于經(jīng)驗模態(tài)分解（empirical mode decomposition，EMD）和支持向量機（support vector machine，SVM）模型對網(wǎng)絡流量進行預測，比單獨使用SVM模型具有較高的準確度。文獻[18]提出小波分解聯(lián)合BP神經(jīng)網(wǎng)絡的方法對網(wǎng)絡流量進行預測，通過使用一種自適應學習率的方法對BP進行優(yōu)化以提高預測準確度。

可以看出，以上研究大多對基站網(wǎng)絡流量進行短期預測，缺乏長期預測研究，而只有對基站網(wǎng)絡流量進行長期預測，才能更有效地指導基站的擴容。此外，基站網(wǎng)絡流量數(shù)據(jù)存在較大的噪聲以及數(shù)據(jù)缺失現(xiàn)象，當前方法普遍未全面考慮這些情況，導致難以保證實際預測效果。

本文提出了一種結合EEMD與Prophet模型的基站小區(qū)網(wǎng)絡流量長期預測方法（E-Prophet）。EEMD是對EMD分解方法的一種改進，解決了EMD因極值點分布不均勻導致的“模式混合”現(xiàn)象，通過EEMD分解能夠對信號進行平穩(wěn)化處理，將時間序列分解為具有不同時間尺度的分量序列。相比小波分解，EEMD不依賴基函數(shù)且具有很好的自適應性。Prophet模型[19]能把時間序列分解成不同的成分，其中趨勢項成分對于具有長期趨勢性的時間序列預測具有較好的效果，季節(jié)項成分對于具有周期性（年、月和星期）的時間序列具有較好的效果。Prophet模型已經(jīng)應用于時間序列預測領域，例如，2019新冠肺炎的分析和預測[20]，以及對LTE（long time evolution）小區(qū)流量增長和擁塞預測[21]等。通過觀察，發(fā)現(xiàn)基站小區(qū)網(wǎng)絡流量數(shù)據(jù)具有一定周期性，與此同時由于小波分解采用不同的小波基和分解層數(shù)對結果影響很大，而且分解效果不如EMD[22]理想，因此本文提出基于EMD/EEMD分解結合Prophet模型對基站小區(qū)網(wǎng)絡流量進行預測，其主要過程包括兩步：首先通過EMD/EEMD分解模型對基站小區(qū)網(wǎng)絡流量進行分解；其次使用Prophet模型對各個分量進行預測，并累加求和各個分量最終預測值得出預測結果。結合EEMD對網(wǎng)絡流量良好的分解特性，以及Prohpet模型對各分量時間序列準確的預測特性，E-Prophet方法能更好地適應于基站小區(qū)網(wǎng)絡流量的長期預測。

本文的主要貢獻為引入Prophet模型到小區(qū)流量預測，Prophet模型雖然已被應用于多種時間序列預測中，但還沒有研究將Prophet用于基站小區(qū)網(wǎng)絡流量預測。通過真實數(shù)據(jù)的實驗證明E-Prophet比其他傳統(tǒng)模型具有更好的預測效果。

1 基站小區(qū)網(wǎng)絡流量的長周期序列特性

Paxson和Floyd的研究發(fā)現(xiàn)[23]，小時級以上粗時間粒度的流量行為由于主要受外界因素的影響，不具有自相似性，而是一種非線性復雜的過程。本文研究的基站小區(qū)網(wǎng)絡流量數(shù)據(jù)是按天統(tǒng)計，每天的流量數(shù)據(jù)為自忙時的流量，所謂自忙時是指當天流量最大的整點時刻。因此基站小區(qū)流量是以小時粒度采集，具有非線性特點，而且有時候因為設備或者人為因素，網(wǎng)絡流量在采集過程中可能出現(xiàn)部分缺失值，這也為基站小區(qū)網(wǎng)絡流量長期預測帶來很大挑戰(zhàn)。同時，用戶在使用網(wǎng)絡服務過程中受外界因素影響，網(wǎng)絡行為具有一定的規(guī)律性和偶然性，因而產(chǎn)生的非線性宏觀網(wǎng)絡流量具有一定的周期性、突發(fā)性和偶然性。根據(jù)基站小區(qū)網(wǎng)絡流量這些特性，網(wǎng)絡流量X(t)可分解為趨勢部分A(t)、周期部分P(t)、突發(fā)部分B(t)和隨機部分R(t),如式（1）所示[24]：

Prophet是一個類似GAM（generative additive model）的時間序列模型，模型本身采用加法模型，分別由趨勢項擬合函數(shù)、季節(jié)項擬合函數(shù)、節(jié)假日擬合函數(shù)和隨機項組成，非常適合具有上述分布特點的網(wǎng)絡流量進行長期預測，具體表現(xiàn)為：Prophet模型趨勢項和季節(jié)項能夠分別擬合網(wǎng)絡流量非周期性和周期性的長期變化規(guī)律；Prophet模型節(jié)假日項能夠對因特殊節(jié)日而導致的有規(guī)律的網(wǎng)絡流量暴跌或暴增進行擬合；同時如果數(shù)據(jù)中存在部分缺失片段，傳統(tǒng)自回歸模型如ARIMA則需要人為插入缺失值后才能進行預測，而Prophet能夠根據(jù)已有的數(shù)據(jù)進行擬合，忽略缺失值對模型預測帶來的影響。此外，通過EEMD分解后的流量序列成分較簡單，使用Prophet模型后能達到更好的預測效果。

2 預測流程及方法

2.1 E-Prophet預測流程

基于E-Prophet預測流程方法如圖1所示。首先，對每個基站原始小區(qū)網(wǎng)絡流量進行EEMD分解；然后，對各個分量使用Prophet模型進行預測；最后，再對各個分量預測結果進行求和，得到該基站小區(qū)網(wǎng)絡流量的最終預測結果。

圖1 E-Prophet的網(wǎng)絡流量預測流程Fig.1 Network traffic forecasting process of E-Prophet

2.2 EMD和EEMD分解

經(jīng)驗模態(tài)分解（EMD）[25]是一種數(shù)據(jù)處理方法，用于處理非線性和非平穩(wěn)數(shù)據(jù)。EMD通過希爾伯特變換將非線性和非平穩(wěn)時間序列分解為分量子序列，包括多個IMF分量和一個殘差項，而無需任何基函數(shù)或濾波函數(shù)。

IMF分量需滿足兩個判定條件：第一，在整個數(shù)據(jù)序列中，極值點的數(shù)量與過零點的數(shù)量相等或最多相差一個，所謂過零點是曲線經(jīng)過縱坐標等于零的點。第二，在任一時間點上，由局部極大值點確定的極大值包絡線，和由局部極小值點確定的極小值包絡線的均值為零，即信號需滿足對時間軸的局部對稱性。EMD的步驟如下[26]：

（1）令{x(t)∈X:t=1,2,…,N}表示基站小區(qū)網(wǎng)絡流量訓練集時間序列，其中t表示天數(shù)。根據(jù)時間序列確定所有局部極值，并將所有的極大值和極小值通過三次樣條線（cubic spline line）連接[27]起來構造上包絡線xmax(t)和下包絡線xmin(t)。

（2）上包絡線和下包絡線的均值m(t)由式（2）得出：

（3）從原始時間序列中減去m(t)以獲得分量l(t)：

（4）檢查l(t)序列是否滿足上述中作為IMF分量的兩個條件。如果滿足這兩個條件，則用l(t)替換x(t)；如果不滿足，則重復步驟（1）～（3），直到l(t)滿足作為IMF分量的兩個判定條件。此時，l(t)為一個IMF分量，令l(t)等于Imf1(t)。通過式（4）從原始流量序列x(t)減去Imf1(t)以得到殘差項r1(t)：

然后，將r1(t)作為更新后的新原始流量序列。重復步驟（1）～（3），獲得Imf2(t),Imf3(t),…,Imfn(t)分量，最后得到殘差項序列rn(t)。原始流量序列與各分量之間的關系見式（5）：

通過上述過程，原始的基站小區(qū)網(wǎng)絡流量序列可以分解為若干個IMF分量和一個殘差項分量。基站小區(qū)網(wǎng)絡流量波動性較大，導致在使用EMD分解時容易發(fā)生“模式混合”現(xiàn)象，為了解決EMD存在的“模式混合”問題[28]，Wu和Huang[29]提出了集合經(jīng)驗模式分解方法。EEMD的思想是利用白噪聲具有頻率均勻分布的特性，把高斯白噪聲添加到信號中，使信號在不同尺度上具有連續(xù)性，從而改變了信號極值點特性，促進抗混分解，有效避免模式混合現(xiàn)象。為此，本文采用EEMD取代EMD來分解流量信號。EEMD分解步驟如下[30]：

（1）初始化EEMD中的參數(shù)，如白噪聲幅度。

（2）在初始基站小區(qū)網(wǎng)絡流量x(t)中添加M次白噪聲：

其中，wj(t)表示第j次添加白噪聲序列，xj(t)表示第j次添加白噪聲后網(wǎng)絡流量。

（3）通過EMD將添加了白噪聲的網(wǎng)絡流量序列xj(t)分解為n個IMF分量Imfj i(i=1,2,…,n)和殘差項分量之和。

（4）計算添加M次白噪聲后得到流量的IMF平均值，作為最終IMF分量。

2.3 Prophet模型

Prophet使用pyStan開源工具中的L-BFGS（limited memory Broyden-Fletcher-Goldfarb-Shanno）方法[31]進行函數(shù)的擬合。L-BFGS是解無約束非線性規(guī)劃問題最常用的方法。Prophet模型通過式（10）對時間序列進行擬合：

其中，y(t)∈{Imf1,Imf2,…,Imfn,rn(t)}；g(t)為趨勢項函數(shù)，適合非周期性序列建模；s(t)為季節(jié)項（周、月、季度和年等）函數(shù)，用于對周期性變化明顯的序列進行建模；h(t)為某固定時間段（例如節(jié)假日等）引起的變化函數(shù)；?t表示噪聲項。

本文使用簡單的分段常數(shù)函數(shù)對趨勢項g(t)進行擬合：

其中，k～Normal(0,5)為流量增長率參數(shù)δ(δ∈?S)為調(diào)整增長速率的增長率向量，m～Normal(0,5)為偏移量參數(shù)。

季節(jié)項s(t)通過式（12）來擬合，如果流量數(shù)據(jù)是以天為單位，那么數(shù)據(jù)可能存在年周期性和星期周期性。

由于本文中基站小區(qū)網(wǎng)絡流量數(shù)據(jù)的時間跨度小于一年，因此不具有年周期性而具有星期周期性。式（12）中P表示周期，當設定時間序列按照幾天為周期，則設定P等于幾，在此設定按星期為周期，所以將P設為7。N作為低通濾波器參數(shù)，對于具有星期周期性的序列，當N設置為3時擬合效果較好[19]。

關于節(jié)假日項h(t)的擬合，假如有L個節(jié)假日，對于第i個節(jié)假日，使用Di表示該節(jié)假日所在的一段時間區(qū)間，參數(shù)κi來表示節(jié)假日時間區(qū)間Di的影響程度，則h(t)可表示為：

噪聲項?t是一個服從正太分布的一個隨機值。

3 實驗過程

3.1 實驗數(shù)據(jù)說明

實驗中基站小區(qū)網(wǎng)絡流量數(shù)據(jù)來源于某區(qū)域1 000個基站小區(qū)網(wǎng)絡流量，單位為KB。流量數(shù)據(jù)按天統(tǒng)計，每天的流量數(shù)據(jù)為自忙時的流量。流量數(shù)據(jù)的時間范圍為2017年9月1日到2018年3月31日，選取2017/9/1—2018/2/28共175天的有效流量數(shù)據(jù)作為模型的訓練集，選取2018/3/1—2018/3/31共31天的流量數(shù)據(jù)作為模型的測試集。

3.2 預測性能評估指標

鑒于網(wǎng)絡流量數(shù)據(jù)數(shù)量級較大，本文選擇均方對數(shù)誤差根（RMSLE）指標來評價模型預測性能，該指標越小表示預測越準確。

其中，N表示要預測流量的天數(shù)，表示網(wǎng)絡流量預測值，yi表示網(wǎng)絡流量真實值；當評價多個小區(qū)預測結果時，N為要預測流量的天數(shù)與基站小區(qū)數(shù)量的乘積。

3.3 實驗設置

本文將E-Prophet模型與SARIMA、LSTM、Prophet以及EMD-Prophet四個模型進行對比。Prophet模型如前所述，其他三個模型說明如下：

SARIMA模型：季節(jié)性差分自回歸滑動平均模型（SARIMA）通過對時間序列進行變換和擬合來建模。其中ARIMA模型在擬合平穩(wěn)過程時表現(xiàn)出很強的適用性，同時對于數(shù)據(jù)的時間趨勢可以進行差分變換，使其轉化為平穩(wěn)的時間序列數(shù)據(jù)，相應地得到預測結果后也要進行差分還原。

LSTM模型[32]：LSTM是循環(huán)神經(jīng)網(wǎng)絡的一種類型，它通過記憶長時間的信息來解決長期依賴性問題。LSTM網(wǎng)絡由LSTM單元組成，每個LSTM單元由以下幾個部分組成：（1）一個細胞，（2）一個輸入門，（3）一個輸出門，（4）一個遺忘門。通過對神經(jīng)網(wǎng)絡設定合理的輸入輸出步長來訓練網(wǎng)絡參數(shù)。

EMD-Prophet模型：與E-Prophet模型類似，差別在于利用EMD對基站小區(qū)網(wǎng)絡流量進行分解。

以上模型的參數(shù)設置為：（1）對于SARIMA模型，通過網(wǎng)格搜索的方式[33]來確定其參數(shù)，預測步長設為31天。（2）對于LSTM模型，輸入步長設為18，輸出步長設為3。由于預測時長為31天，通過滑動窗口的方式把預測的3天的數(shù)據(jù)作為輸入，以此類推，迭代11次預測未來33天數(shù)據(jù)，并取前31天數(shù)據(jù)值作為最終預測結果。（3）對于Prophet模型，使用星期周期和月周期參數(shù)，預測步長參數(shù)為31天。

通過設置以下實驗場景對各種模型進行比較和驗證：場景一，使用EMD和EEMD方法對基站流量進行分解，以此比較兩種方法的流量分解效果；場景二，在單基站長期預測中對各個模型預測性能進行比較，以此觀察不同預測模型的擬合效果；場景三，在多基站長期預測中對各個模型預測性能進行比較，以此來對比不同模型長期預性能；場景四，在缺失數(shù)據(jù)情況下對不同模型預測性能進行比較，以此分析模型的魯棒性。

實驗環(huán)境的硬件配置：Intel?CoreTMi7-6700 CPU@3.40 GHz，22核，64 GB內(nèi)存。軟件配置為：Python3.7.3，tensorflow1.15，fbprophet0.5，EMD-signal[34]0.2.10。

3.4 EMD/EEMD流量分解效果對比

為了比較EMD和EEMD對于流量分解的效果，隨機從1 000個基站小區(qū)選取一個基站作為案例研究。原始基站小區(qū)網(wǎng)絡流量數(shù)據(jù)如圖2所示，圖3和圖4分別表示EMD和EEMD對原始數(shù)據(jù)流量進行分解的效果。

圖2 某基站原始網(wǎng)絡流量訓練集Fig.2 Raw network traffic training set for base station

圖3 某基站原始網(wǎng)絡流量的EMD分解Fig.3 EMD decomposition of base station’s raw network traffic training set

通過圖5可以看出，原始的網(wǎng)絡流量按星期進行劃分之后，周六、日的采集流量值普遍比其他星期采集的流量值小。通過圖6可得，Imf1分量不同星期流量都比原始分量更加平穩(wěn)。通過圖7可得，Imf2分量不同星期流量值分布區(qū)間有顯著區(qū)別。

圖5 某基站原始網(wǎng)絡流量按星期劃分Fig.5 Raw network traffic training set of base station divided by week

通過圖8可得，Imf3分量不同星期對應的網(wǎng)絡流量具有相似的周期性。通過圖9～11可得，其他分量不同星期所對應的網(wǎng)絡流量具有相似的趨勢性，而且不同星期取值呈現(xiàn)一定的按星期周期交換的規(guī)律，即不同星期對應的流量大小順序，會隨著星期間隔變化而呈現(xiàn)出一定的周期性。這是因為原始數(shù)據(jù)流量經(jīng)過EEMD分解后的各個分量具有一定的平穩(wěn)性、趨勢性以及按星期周期變化的規(guī)律。

圖6 Imf1分量按星期劃分Fig.6 Imf1 components by day of week

圖7 Imf2分量按星期劃分Fig.7 Imf2 components by day of week

圖8 Imf3分量按星期劃分Fig.8 Imf3 components by day of week

綜上可看出，相對于EMD，EEMD能使分解的序列表現(xiàn)出更好的平滑性、平穩(wěn)性以及趨勢性和周期性，更加適合Prophet模型的擬合。

3.5 單基站不同模型預測結果對比分析

本節(jié)對3.4節(jié)案例研究中的基站小區(qū)使用不同模型進行預測，預測步長為31天，結果對比如圖12所示。

圖9 Imf4分量按星期劃分Fig.9 Imf4 components by day of week

圖10 Imf5分量按星期劃分Fig.10 Imf5 components by day of week

圖11 殘差項分量按星期劃分Fig.11 Residual components by week

圖12 某基站不同模型預測結果對比Fig.12 Comparison of prediction results of different models of base station

由圖12可看出，真實流量值有一定的周期性，而Prophet、EMD-Prophet和E-Prophet能識別這種周期變化的特性和趨勢。SARIMA長期預測結果優(yōu)于Prophet，但從圖中預測結果看，在長期預測中SARIMA擬合趨于一條直線。由于使用EMD/EEMD分解后的數(shù)據(jù)平穩(wěn)性更好，因此結合EMD/EEMD和Prophet模型的預測結果優(yōu)于Prophet模型。此外，與EMD相比，EEMD避免了因極值點分布不均勻導致的“模式混合”現(xiàn)象，而且分解后的殘差項趨勢性更顯著。因此使用E-Prophet（EEMD+Prophet）預測效果優(yōu)于EMD-Prophet模型。分別計算上述不同模型預測結果的RMSLE值，如表1所示。由表1可看出E-Prophet模型的預測性能最好。

表1 不同模型預測某基站小區(qū)流量的RMSLE值Table 1 RMSLE value of base station cell traffic predicted by different models

3.6 多基站不同模型長期預測結果分析

上節(jié)對單基站的預測性能并不能代表模型的整體預測性能，本節(jié)對多基站網(wǎng)絡流量進行預測。為此，對875個基站小區(qū)流量進行了長期預測。需要說明的是，在1 000個小區(qū)中，有些小區(qū)具有較多的零值網(wǎng)絡流量，SARIMA無法利用網(wǎng)格搜索找到最佳參數(shù)導致無法預測。為方便對比，選取SARIMA能夠預測的所有875個基站小區(qū)流量。預測步長范圍為16～31 d。預測結果對比和統(tǒng)計分別如圖13和表2所示。

圖13 多基站不同模型長期預測結果對比Fig.13 Comparison of long-term prediction results of different models for multiple base stations

綜合圖13和表2可以看出，EMD/EEMD結合Prophet的模型在長期預測中優(yōu)于Prophet、SARIMA和LSTM模型。這是因為EMD和EEMD方法可以把時間序列分解成平穩(wěn)性和趨勢性更顯著的分量，同時保持原有時間序列的周期特征，最后通過Prophet模型對這些分量進行更準確的預測。

表2 多基站長期預測的RMSLE值對比Table 2 Statistical comparison of long-term forecast RMSLE data for multiple base stations

3.7 魯棒性實驗結果分析

為了驗證E-Prophet模型的魯棒性，對不同模型在缺失數(shù)據(jù)環(huán)境下的預測性能進行對比。選取125個零值數(shù)據(jù)較多的基站進行實驗，結果見圖14和表3。

圖14 多基站不同模型魯棒性結果對比Fig.14 Comparison of robustness results of different models of multiple base stations

從中可看出，結合EMD/EEMD和Prophet的模型在長期預測中優(yōu)于Prophet模型、LSTM模型、EMD-Prophet模型。這是因為Prophet模型能通過計算訓練集數(shù)據(jù)概率分布和周期性來確定相應參數(shù)，從而在缺失值較多的小區(qū)依然能較好預測。SARIMA對于以上零值較多的基站小區(qū)，無法得出預測結果。

4 結語

由于網(wǎng)絡流量具有非線性、非平穩(wěn)的特點，因此難以準確地進行長期預測。本文提出基于Prophet混合模型的基站網(wǎng)絡流量長期預測方法，通過實驗得出以下結論：（1）EEMD可將流量序列分解為一組平穩(wěn)或較平滑的分量，適合處理基站小區(qū)網(wǎng)絡流量序列；（2）EEMD使得每個IMF分量和殘差項具有清晰的物理含義、IMF分量更平穩(wěn)，從而使得Prophet模型對每個IMF分量和殘差項有更好的預測效果，并且采用EEMD比采用EMD的預測準確性更高；（3）EMD/EEMD結合Prophet的組合模型在長期預測中比Prophet、LSTM、SARIMA模型的準確性、魯棒性更好，并且E-Prophet優(yōu)于EMD-Prophet模型。下一步研究工作將探索將概率模型算法和深度學習模型相結合的方法來解決多區(qū)域多基站聯(lián)合預測的問題。

表3 多基站魯棒性實驗的RMSLE值對比Table 3 Statistical comparison of RMSLE data for multisite robustness experiments