Deepar
- 2021 年 6 月 22 日
- AI
deepar這篇論文的亮點不僅在於提出了新的模型結構，許多關於時序預測的見解也很具有啟發性，下面是對論文的閱讀筆記和思考：
1、
在現實世界的預測問題中，試圖共同學習多個時間序列時經常遇到的一個挑戰是，時間序列的數量級差異很大，而且數量級的分布具有很強的傾斜性。這個問題如圖1所示，圖中顯示了亞馬遜銷售的數百萬件商品的銷售速度(即平均每周銷售一件商品)的分布情況。分布在幾個數量級上，近似冪律。
據我們所知，這一發現是新的(儘管可能並不令人驚訝)，並且對試圖從這些數據集學習全局模型的預測方法具有基本意義。由於該分布的無標度特性，很難將數據集劃分為具有一定速度帶的時間序列子組，並為它們學習單獨的模型，因為每個速度子組都會有相似的偏態。此外，基於分組的正則化方案，比如Chapados[4]提出的方案，可能會失敗，因為每個分組的速度將有很大的不同。最後，這種偏態分布使得使用某些常用的標準化技術，如輸入標準化或批標準化[14]，效果較差。
這裡說的是最常見的多序列預測問題，web traffic，rossman store forecasting等比賽以及電商的商品銷量預測的都是典型的多序列預測問題，不同商品的base level（時序數據四分量之一）差異很大，整體也呈現冪律分布，常規標準化處理無效，因為常規的標準化基於全局的統計特徵無法改變不同商品之間的相對大小關係，使得不同商品的base level差異仍舊很大，冷門商品的銷量和熱門商品的銷量相差成千上萬倍，這導致輸入數據的量綱差異極大（一般基本是做分組標準化），即不同的商品自己內部做標準化，從而使得不同商品的序列數據具有可比性。
2、模型說明部分
為了便於說明，這裡我們以電商銷量預測的問題為例：
問題定義：
其中：
表示已知的，從時刻1~時刻t0 -1的商品的銷量數據
表示從時刻1到時刻T的所有的協變數，這裡需要說明，我們假設[1, t0− 1]為已知銷量的時間區間，[t0, T]為我們需要預測的時間區間，而Xi,1:T是橫跨了[1, t0− 1]和[t0, T]兩個時間區間的，這意味著X 必然是 靜態特徵（例如商品的大類，商品所在的省，城市或商品的歷史銷量的統計值）或者是可以推斷的動態特徵，time varying known features，例如月份，年份，隨著時間變化，但是他們未來的值是可以直接推算出來的，而銷量這種數據屬於不可推斷的動態特徵（同時也是標籤），即time varying unknown features。
模型結構圖如上，在理解上圖之前，需要了解兩個背景知識：
1、時間序列預測中的自回歸模型：時間序列預測中進行自回歸式的預測本質上是把多步預測轉化為多個遞歸式的單步預測，一個典型的例子就是基於RNN的seq2seq結構。seq2seq的結構有兩種形式：
簡單的seq2seq的結構是無法使用預測的數據的特徵，僅僅是將t時刻的隱層輸出作為decoder部分的共同輸入。
實際上上圖的seq2seq結構才算是在做真正意義上的自回歸，t+1時刻的預測結果被用於幫助預測t+2時刻的結果，t+2時刻的預測結果被用於幫助預測t+3時刻的結果。
當然，在這個過程中，外生變數是可以放進來的，假設我們的encoder部分的輸入的序列數據是[5,4,3,2,1,0]對應t-5到t，外生變數是[a5,b5,c5]….[a0,b0,c0]，則encoder部分的輸入數據的形式就是[5,a5,b5,c5]…..[0,a0,b0,c0]，decoder部分的初始輸入是t時刻的特徵即[0,a0,b0,c0]，我們在decoder的t時刻部分產生了t+1時刻的預測 t+1_pred，然後我們將t+1_pred和t+1時刻的可推斷特徵（例如城市，省份，店鋪，品類，月份，星期等）合併得到了[t+1_pred,a-1,b-1,c-1]，作為t+2時刻的輸入特徵幫助預測t+2，依次類推。
由於這種遞歸式的訓練很容易產生誤差累積的問題，即如果t+1的預測誤差大，則後續的decoder都是在錯誤的預測結果上繼續訓練的，所以更常見的訓練方式是使用teacher forcing：
decoder的時候直接使用真實的值幫助模型訓練。
而常規的LSTM是直接做的向量輸出：
seq2seq式的encoder-decoder結構相對於LSTM這樣單一的encoder的結構有兩個好處：
1、』seq2seq可以支援不定長的輸入和輸出，而單一的encoder結構只能支援不定長輸入無法支援不定長的輸出，因為單一的encoder結構最終是通過一個固定的dense層完成預測任務的，dense層的size是固定不變的；
2、seq2seq的訓練過程可以使用到未來的特徵，當然這些特徵必須是可以推斷的（前面描述過過了），這可能可以幫助模型更好的學習。
deepar的模型結構介於常規的RNN和基於RNN的seq2seq之間，這在原文中提到了：
即deepar中，encoder和decoder使用的是同一個RNN模型，這個RNN模型就是常規意義上的RNN模型，可以是單層的LSTM或多層的LSTM。
那麼deepar的模型是怎麼訓練的？和使用teacher forcing的seq2seq一樣，deepar的訓練和預測過程是不一樣的，基於teacher forcing訓練策略的seq2seq訓練的時候用到了真實的標籤，但是預測的時候我們是沒有真實標籤可用的，因此這類seq2seq預測的行為退化為了普通的recursive seq2seq：
在程式碼實現上也比較麻煩，我們需要手動寫decoder部分的循環（不知道有沒有什麼方便的開源可以實現這個功能，沒有找過）
deepar的訓練過程如上，其實去掉中間的
部分，deepar看起來就是一個LSTM，為了方便描述，中間部分暫時先不考慮，看一下deepar的輸入和輸出是什麼
輸入部分：Xi(t-1)，即t-1時刻的外生變數，即可推斷的特徵，比如前面說過的城市，省份，店鋪，品類，月份，星期等，Zi（t-2），即t-2時刻的標籤，比如商品的銷量（也就是我們想要預測的目標），通過一個RNN cell得到 t-1時刻的表徵，然後接一個dense（1）的層得到Zi(t-2)，依次類推。其中t-1，t和t+1都是我們需要預測的標籤，和teacher forcing seq2seq一樣，訓練的過程中，待預測的標籤被放進來用於模型訓練了。
然後在預測階段：
和recursive seq2seq一樣，進行迭代預測，這種預測方式可以預測無限長的序列，只不過當預測的序列變長之後精度往往很差，這也是時序預測問題的一個挑戰，長時序的預測如何保證精度的問題，和本文無關，不廢話。
2、
當然deepar沒這麼簡單，deepar的設計本身不是針對點預測，而是針對概率預測的。主要問題在於，怎麼用nn來對標籤的概率分布直接進行建模。
如果熟悉貝葉斯深度學習，對這一套操作應該會比較熟悉了：
1、首先，人工指定一個標籤的概率分布，這裡我們使用常見的一元高斯分布（因為這裡我們僅僅預測一維的數據，而不像天氣預報那樣要預測溫度，濕度，風力等的多個維度的數據）：
一元高斯分布的概率密度函數如上，也就是論文中提到的：
2、問題轉換，現在我們的問題不是直接擬合標籤的值z了，而是要去擬合我們假定的標籤z服從的高斯分布的參數mu和sigma，那麼怎麼設計呢？思路非常的簡單，簡單的令人髮指，以下圖為例：
我們要做的事情就是，把network的輸出，以LSTM為例，假設network使用的是LSTM，我們令LSTM return sequences，然後LSTM的每一個時間步的hidden state粉筆接兩個Dense(1)，我們假設為Dense A(1) 和 Dense B(1)，隱層的神經元只有一個，源程式碼中的實現對應的就是：
Python"> mu = K.dot(inputs, self.W_mu)
mu = K.bias_add(mu, self.b_mu, data_format='channels_last')
sigma = K.dot(inputs, self.W_sigma)
sigma = K.bias_add(sigma, self.b_sigma, data_format='channels_last')
sigma = K.softplus(sigma) + K.epsilon()這裡的mu表示我們要求解的高斯分布的 參數之一，即U，均值，我們直接將LSTM的hidden state 接DenseA(1)得到1個1維的向量作為對高斯分布的均值U的估計，然後如法炮製我們將LSTM的hidden state又接了一個 DenseB(1)得到高斯分布的標準差 sigma，因為標準差必然是正數，所以為了避免參數訓練的過程中sigma變成負數，這裡加入了
對sigma的值進行調整。
這樣，deepar的整個前向傳播的過程就定義完畢了，對應論文中的：
總結一下，對於一個時序樣本來說，這個時序樣本包含了n個time step的時間切片，進入了network（這個network是可以處理序列數據的網路結構例如LSTM，CNN，seq2seq的網路機構或者Transformer）處理之後，吐出了n個新的序列，然後這個序列分別接兩個dense層，兩個dence層的神經元數量均為1，一個dense層用於擬合人為指定的高斯分布的u均值，一個用於擬合認為指定的高斯分布的標準差sigma，然後為了避免sigma小於0，使用softplus，即log(1+e**x)來保證sigma大於0。這裡的兩個dense層被封裝到一個layer，我們稱這個layer為gaussian layer，專門用於擬合高斯分布的特殊網路層。
3、模型訓練：
原文的實驗設計有幾個地方需要注意：
第一：缺失值用0填充
第二：在nlp中，seq2seq使用teacher forcing有一個問題，就是模型訓練和預測行為的不一致會導致模型的泛化性能下降，但是deepar的論文里提到他們沒有遇到這個問題。
第三：關於縮放的問題，考慮到模型的輸入是多個時間序列，這些時間序列的量級可能並不一樣，因此我們需要對它們做放縮，對於每一個時間序列i，都對應有一個放縮因子（from 段易通：概率自回歸預測——DeepAR模型淺析）
可以看到，這裡的處理其實就是每個時間序列做分組標準化，每個時序數據除以自己的均值，為了避免均值為0，所以對均值+1避免除0的問題。
原文中還提到了關於負二項似然的分布假設，沒怎麼用過，所以不想看了。。到時候有需要再看好了，問題不大。
第四：關於取樣的問題
這個工程技巧很有意思，原文中的意思是，我們對於銷量數據比較全，也就是比較熱門，上了比較久的商品的預測更看重，而對於冷門商品，銷量數據很少或者是新商品則其重要性相對低一些，所以選擇有權取樣的方式進行取樣，關於時序數據的取樣，文中的意思是這樣的，比如假設我們的training window是100，prediction window是20，那麼實際上每一次入模的樣本總的長度是120，正常的操作是我們要從某個商品的歷史數據中對數據進行切分，比如某個商品的時序數據的長度一共是1000，則這個商品的序列數據可以拆分為多條樣本，按照文中加權取樣的意思，對於熱門的商品取樣的權重更大，即這種商品的序列數據會更頻繁的出現在訓練數據中。
第五：特徵工程
特徵工程這塊兒就沒什麼特別的，日期特徵，一些靜態特徵的處理之類的和常規的時序數據的特徵工程沒有太大區別，最後都進行了standardscaler標準化
第六：損失函數是啥？
我們前向傳播，將每一個樣本映射為了 u和sigma兩個值，那麼怎麼和實際的標籤y關聯起來，顯然我們這個時候沒有辦法使用常規的用於回歸的loss，比如mae，mape，mse之類的，因為這類損失函數無法處理 u，sigma，y_true 三個輸入的關係。這塊兒的loss設計也非常直觀，我們高中的時候，高斯分布的均值和方差是有公式的，其推導過程如下
就是最大似然估計，如果從nn的角度來看，我們就是希望能夠讓：
這個式子越小越好，那好辦，我們直接帶入這個公式就可以了，於是就有
def gaussian_loss(y_true, mu, sigma):
loss = (
tf.math.log(sigma) # 上述公式第一項
+ 0.5 * tf.math.log(2 * np.pi) #上述公式第二項，這項是常數項可以省去
+ 0.5 * tf.square(tf.math.divide(y_true - mu, sigma)) #這裡和公式有兩個區別
#一個是 y_true - mu 這項應該平方，這裡直接用1次方，個人認為主要是擔心y_true和mu差異太大，
#平方之後的差異更大，導致loss的變動範圍太大所以直接省去了二次項，對公式中的第三項用了
# tf.sqaure進行放縮，同樣也是為了避免 除以sigma之後的結果太大的問題
#
)
return tf.reduce_mean(loss)（這裡的設計是把sigma當作方差而不是標準差）
最後，前向傳播和損失函數都已經完事兒了，來看下模型是怎麼預測的。
現在我們直接進行predict一個樣本a1，得到的不是y1，而是y1所在的高斯分布的u和sigma（對於每一個y，其所在的高斯分布都是不同的，這一點非常重要，我們是認為每一個觀測點服從的高斯分布都是不同的，（學完deepar，對貝葉斯統計和貝葉斯深度學習的理解會有一定的幫助），如果某個y的sigma很大，則這個y的預測的不確定性很高，我們可以很容易畫出不確定性的預測結果如下圖：
）
一般來說，我們得到了高斯分布的均值和方差，則直觀的思路是用高斯分布的期望值，
但是實際上，我們的所謂的高斯分布的估計都僅僅是建立在一個樣本上的而已，無論是訓練還是預測都一樣，一個樣本就對應一個mu和一個sigma，即一個樣本就對應一個特定的基於該樣本的高斯分布，直接使用均值會一個問題：
這個均值不準，根據期望公式可以知道，高斯分布的期望的計算其實就是取樣一個x，計算這個x的概率密度，然後相乘，然後累加進來，其實就是加權求和，因為x的取值是無限的在正負無窮之間，所以通過積分公式來計算，但是實際上我們不必高斯分布的取值範圍里所有的點，我們只需要使用置信度高的點就可以了，那麼實現的方法也很簡單：
def gaussian_sample(mu, sigma):
mu = tf.expand_dims(mu, axis=-1)
sigma = tf.expand_dims(sigma, axis=-1)
samples = tf.random.normal((300,), mean=mu, stddev=sigma)
return tf.reduce_mean(samples, axis=-1)既然有了均值和方差，那麼就從這個分布里隨機取樣出一些點直接做平均。這裡的程式碼設計是取樣300個點。實際使用的話，可以調整。
個人使用的一些感受:
deepar的設計是非常靈活的，但是局限的地方在於，比如上面，我們優化的是gaussian_loss，並不是直接優化mape，rmse這類的損失函數，所以相對來說，沒有辦法直接去優化我們的業務指標。除此 之外，如果假定的其它分布，則最終的輸出層需要做修改，比如possion分布，tweetie 分布，需要實現了解這些分布的參數，然後相應地去設計輸出層。
最後是git上的一些源程式碼實現。
deepar的實現git上不太多，
gluonts里的核心就使用了rnn的encoder層而沒有decoder層，
//github.com/arrigonialberto86/deepar/blob/master/deepar/model/lstm.py
tf的實現和論文的圖比較契合，不過這個實現寫的不是很好，封裝的方式比較不舒服，而且也沒有實現原文中scaling的部分，
//www.youtube.com/watch?v=OP4wa3cFSG8&list=PLmFBTSyQ2rzUy1K929-6RRfFcfDZPUq6U&index=4
最後在youtube的一個公開的講座上找到了相應的實現，整體的模型設計比較keras原生態，移植起來比較方便，不過用到了encoder+decoder兩個結構。
其實就論文原圖來說，中間network的設計可以很靈活，lstm，tnn，transformer其實都可以作為這裡的network，只要加入的網路組件是一個 seq2seq式的並且具有記憶功能的模型結構就可以，即輸入n個序列，吐出n個序列的模型結構。
直接來看下這段程式碼：
class GaussianLayer(keras.layers.Layer):
def __init__(self, **kwargs):
super(GaussianLayer, self).__init__(**kwargs)
self.W_mu = None
self.b_mu = None
self.W_sigma = None
self.b_sigma = None
def build(self, input_shape):
super(GaussianLayer, self).build(input_shape)
dim = input_shape[-1]
self.W_mu = self.add_weight(
name='W_mu',
shape=(dim, 1),
initializer='glorot_normal',
trainable=True,
)
self.b_mu = self.add_weight(
name='b_mu',
shape=(1,),
initializer='zeros',
trainable=True,
)
self.W_sigma = self.add_weight(
name='W_sigma',
shape=(dim, 1),
initializer='glorot_normal',
trainable=True,
)
self.b_sigma = self.add_weight(
name='b_sigma',
shape=(1,),
initializer='zeros',
trainable=True,
)
def call(self, inputs):
mu = K.dot(inputs, self.W_mu)
mu = K.bias_add(mu, self.b_mu, data_format='channels_last')
sigma = K.dot(inputs, self.W_sigma)
sigma = K.bias_add(sigma, self.b_sigma, data_format='channels_last')
sigma = K.softplus(sigma) + K.epsilon()
return tf.squeeze(mu, axis=-1), tf.squeeze(sigma, axis=-1)
def gaussian_loss(y_true, mu, sigma):
loss = (
tf.math.log(sigma)
+ 0.5 * tf.math.log(2 * np.pi)
+ 0.5 * tf.square(tf.math.divide(y_true - mu, sigma))
)
return tf.reduce_mean(loss)
def gaussian_sample(mu, sigma):
mu = tf.expand_dims(mu, axis=-1)
sigma = tf.expand_dims(sigma, axis=-1)
samples = tf.random.normal((300,), mean=mu, stddev=sigma)
return tf.reduce_mean(samples, axis=-1)
class DeepAR(keras.Model):
def __init__(self, fw, rnn_units):
super(DeepAR, self).__init__()
self.encoder_layer = keras.layers.RNN(create_rnn_cells(rnn_units))
self.repeat_layer = keras.layers.RepeatVector(fw)
self.decoder_layer = keras.layers.RNN(create_rnn_cells(rnn_units), return_sequences=True)
self.gaussian_layer = GaussianLayer()
def call(self, inputs):
outputs = self.encoder_layer(inputs)
outputs = self.repeat_layer(outputs)
outputs = self.decoder_layer(outputs)
return self.gaussian_layer(outputs)
def train_step(self, inputs):
x, y = inputs
with tf.GradientTape() as tape:
mu, sigma = self(x)
loss_val = self.loss(y, mu, sigma)
grads = tape.gradient(loss_val, self.trainable_weights)
grads = [tf.clip_by_value(grad, -1., 1.) for grad in grads]
self.optimizer.apply_gradients(zip(grads, self.trainable_weights))
return {'loss': loss_val}
def test_step(self, inputs):
x, y = inputs
mu, sigma = self(x, training=False)
loss_val = self.loss(y, mu, sigma)
return {'loss': loss_val}
def predict_step(self, inputs):
x = inputs
mu, sigma = self(x, training=False)
return gaussian_sample(mu, sigma)補充一個假設y服從泊松分布的：
class PoissonLayer(keras.layers.Layer):
def __init__(self, **kwargs):
super(PoissonLayer, self).__init__(**kwargs)
self.W_rate = None
self.b_rate = None
def build(self, input_shape):
super(PoissonLayer, self).build(input_shape)
dim = input_shape[0][-1]
self.W_rate = self.add_weight(
name='W_rate',
shape=(dim, 1),
initializer='glorot_uniform',
trainable=True,
)
self.b_rate = self.add_weight(
name='b_rate',
shape=(1,),
initializer='zeros',
trainable=True,
)
def call(self, inputs):
x, scale = inputs
rate = K.dot(x, self.W_rate)
rate = K.bias_add(rate, self.b_rate, data_format='channels_last')
rate = K.softplus(rate) + K.epsilon()
rate = rate * scale
return tf.squeeze(rate, axis=-1)
class MeanScalerLayer(keras.layers.Layer):
def call(self, inputs):
scale = tf.reduce_mean(tf.abs(inputs), axis=1, keepdims=True) + 1.0
outputs = inputs / scale ##除以均值
return outputs, scale
def poisson_loss(y_true, rate):
loss = -1.0 * (
y_true * tf.math.log(rate)
- tf.math.lgamma(y_true + 1.0)
- rate
)
return tf.reduce_mean(loss)
def poisson_sample(rate):
sample = tf.random.poisson((300,), lam=rate)
sample = tf.reduce_mean(sample, axis=0)
return tf.expand_dims(sample, axis=-1)
class DeepARPos(keras.Model):
def __init__(self, fw, rnn_units):
super(DeepARPos, self).__init__()
self.mean_scaler_layer = MeanScalerLayer()
self.encoder_layer = keras.layers.RNN(create_rnn_cells(rnn_units))
self.repeat_layer = keras.layers.RepeatVector(fw)
self.decoder_layer = keras.layers.RNN(create_rnn_cells(rnn_units), return_sequences=True)
self.poisson_layer = PoissonLayer()
def call(self, inputs):
outputs, scale = self.mean_scaler_layer(inputs)
outputs = self.encoder_layer(outputs)
outputs = self.repeat_layer(outputs)
outputs = self.decoder_layer(outputs)
return self.poisson_layer([outputs, scale])
def train_step(self, inputs):
x, y = inputs
with tf.GradientTape() as tape:
rate = self(x)
loss_val = self.loss(y, rate)
grads = tape.gradient(loss_val, self.trainable_weights)
grads = [tf.clip_by_value(grad, -1., 1.) for grad in grads]
self.optimizer.apply_gradients(zip(grads, self.trainable_weights))
return {'loss': loss_val}
def test_step(self, inputs):
x, y = inputs
rate = self(x, training=False)
loss_val = self.loss(y, rate)
return {'loss': loss_val}
def predict_step(self, inputs):
x = inputs
rate = self(x, training=False)
return poisson_sample(rate)感覺tensorflow probability裡面應該會有現成的實現吧，還沒看過，再說
其它有意思的資源：
