關於局部加權迴歸（Locally Weighted Scatterplot Smoothing，LOWESS）、STL（Seasonal-Trend decomposition procedure b

時間 2019-11-05

標籤關於局部加權迴歸 locally weighted scatterplot smoothing lowess stl seasonal trend decomposition procedure 简体版

原文原文鏈接

1. LOWESS（Locally Weighted Scatterplot Smoothing，局部加權迴歸）

0x1：lowess算法主要解決什麼問題

1. 非線性迴歸擬合問題

LOWESS 經過取必定比例的局部數據，在這部分子集中擬合多項式迴歸曲線，這樣咱們即可以觀察到數據在局部展示出來的局部規律和局部趨勢（局部均值μ迴歸）。html

同時，由於 LOWESS 的統計窗口縮小爲局部窗口，所以擬合的迴歸曲線能夠包含週期性，波動性的信息。git

和傳統的多元線性迴歸最大的不一樣是，一般的迴歸分析每每是根據全體數據建模，這樣能夠描述總體趨勢，但現實生活中規律不老是（或者不多是）一條直線。github

LOWESS將局部範圍從左往右依次推動，最終一條連續的曲線就被計算出來了。顯然，曲線的光滑程度與咱們選取數據比例（局部窗口）有關：web

局部統計窗口比例越少，擬合越不光滑，由於過於看重局部性質。但相對的，擬合的精確度更高
局部統計窗口比例越大，擬合越光滑，由於更看重總體性質。但擬合的精確度也相應下降

下面經過一個例子來形象說明"LOWESS"和"Linear regression"在局部信息捕獲上的主要區別算法

數據文件是關於「物種數目與海拔高度」（中科院植物所-賴江山博士），數據集能夠在這裏下載。api

用迴歸模型擬合的圖示以下，本文采用R語言：安全

# 從counts.txt文件直接將數據讀入R
x = read.csv("/Users/zhenghan/Downloads/counts.txt")
print(x)

par(las = 1, mar = c(4, 4, 0.1, 0.1), mgp = c(2.5,
                                              1, 0))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
abline(lm(counts ~ altitude, x), lwd = 2, col = "red")

若是僅僅從迴歸直線來看，彷佛是海拔越高，則物種數目越多。如此推斷下去，恐怕月球或火星上該物種最多。這顯然是一個不是那麼完美的線性迴歸推斷。服務器

改用lowess進行局部加權迴歸擬合，圖示以下：網絡

# 從counts.txt文件直接將數據讀入R
x = read.csv("/Users/zhenghan/Downloads/counts.txt")
print(x)

par(las = 1, mar = c(4, 4, 0.1, 0.1))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
# 取不一樣的f參數值
for (i in seq(0.01, 1, length = 100)) {
  lines(lowess(x$altitude, x$counts, f = i), col = gray(i),
        lwd = 1.5)
  Sys.sleep(0.15)
}

上圖中，曲線顏色越淺表示所取數據比例越大，即數據點越集中。不難看出白色的曲線幾乎已呈直線狀，而黑色的線則波動較大。架構

整體看來，圖中大體有四處海拔上的物種數目偏離迴歸直線較嚴重：450 米、550 米、650 米和 700 米附近。這種偏離迴歸的現象剛好說明了實際問題的機率分佈並非嚴格遵循線性迴歸模型，而是其餘的更復雜的非線性模型。

基於LOWESS的擬合結果，若研究者的問題是，多高海拔處的物種數最多？那麼答案應該是在 650 米附近。

爲了確保咱們用 LOWESS 方法獲得的趨勢是穩定的，咱們能夠進一步用 Bootstrap 的方法驗證。由於 Bootstrap 方法是對原樣本進行重抽樣，根據抽得的不一樣樣本能夠獲得不一樣的 LOWESS 曲線，最後咱們把全部的曲線添加到圖中，看所取樣本不一樣是否會使得 LOWESS 有顯著變化。

# 從counts.txt文件直接將數據讀入R
x = read.csv("/Users/zhenghan/Downloads/counts.txt")
print(x)

set.seed(711) # 設定隨機數種子，保證本圖形能夠重製
par(las = 1, mar = c(4, 4, 0.1, 0.1), mgp = c(2.5,
                                              1, 0))
plot(x, pch = 20, col = rgb(0, 0, 0, 0.5))
for (i in 1:400) {
  idx = sample(nrow(x), 300, TRUE) # 有放回抽取300個樣本序號
  lines(lowess(x$altitude[idx], x$counts[idx]), col = rgb(0,
                                                          0, 0, 0.05), lwd = 1.5) # 用半透明顏色，避免線條重疊使得圖形看不清
  Sys.sleep(0.05)
}

能夠看出，雖然 700 米海拔附近物種數目減少的趨勢並不明顯了，這是由於這個海拔附近的觀測樣本量較少，在重抽樣的時候不容易被抽到，所以在圖中表明性不足，最後獲得的擬合曲線分佈稀疏。

可是通過 400 次重抽樣並計算 LOWESS 曲線，剛纔在第一幅圖中觀察到的趨勢大體都還存在（由於默認取數據比例爲 2/3，所以擬合曲線都比較光滑），因此總體迴歸擬合的結果是可信的。

2. 局部平滑問題

因爲樣本數據不可避免的隨機波動問題（方差存在），咱們得到的樣本數據老是存在偏離均值的狀況，爲了更好地體現出數據中的統計規律性，咱們須要對數據進行平滑處理。

舉一個具體例子來講，在作數據平滑的時候，會有遇到有趨勢或者季節性的數據，對於這樣的數據，咱們不能使用簡單的均值正負3倍標準差之外作異常值剔除，須要考慮到趨勢性等條件。使用局部加權迴歸，能夠擬合一條趨勢線，將該線做爲基線，偏離基線距離較遠的則是真正的異常值點。

0x2：從KNN線性迴歸逐漸演進爲LOWESS

上一章節討論了線性迴歸在具體工程問題中面臨的2個主要挑戰，即：

週期性、波動性非線性數據迴歸問題
數據系統性噪點（非數據自己包含的方差，而是系統自己固有的方差）的平滑化問題

面對以上問題，學術界的陸續發表了不少研究paper，很好地解決了這些問題，筆者在這裏儘可能列舉和梳理出整個理論體系的發展脈絡。

1. 多元非線性迴歸模型

人們最先提出，用多元線性迴歸模型來擬合一元線性模型沒法擬合的複雜樣本數據。可是遇到了幾個比較大的問題：

多元模型的模型函數須要人工設定，例如根據數據集繪出的圖像，半觀察半猜想要不要加個自變量的平方項，或者自變量的三角函數等等，做爲新的迴歸算子加入，直到一些評價指標上擬合效果較好爲止，這個過程就是很是繁重的特徵工程工做
多元線性迴歸很容易陷入overfit或者underfit問題。

2. KNN線性迴歸

另外一個解決問題的思路就是對數據進行局部最近鄰處理，也即KNN線性迴歸。

其中：

Nk(x)：爲距離點x最近k個點xi組成的鄰域集合（neighborhood set）
yi：爲KNN均值化後的迴歸結果
：爲統計窗口內的迴歸結果，該統計窗口的全部點都用該點代替

可是，KNN線性迴歸有幾個最大的缺點：

沒有考慮到不一樣距離的鄰近點應有不一樣的權重，無論是從統計上仍是從領域經驗角度來看，偏離均值越遠的點重要性理應越低；
擬合的曲線不連續（discontinuous），擬合曲線的毛躁點仍是不少，例以下圖

關於KNN線性迴歸的問題的其餘討論，能夠參閱這篇文章。

3. 引入kernel加權平滑的KNN線性迴歸

解決傳統KNN線性迴歸問題的一個思考方向是：避免直接求均值這種「hard regression」的方式，改成尋找一種「soft regression」的方式。

咱們在原始的KNN公式之上引入加權平滑：

權重因子函數K能夠由不少種形式，一個大的原則就是：偏離均值越遠的點，其貢獻度相對越低。

好比，Epanechnikov二次kernel：

，

其中，λ爲kernel的參數，稱之爲window width。對於KNN，只考慮最近的k個點影響；基於此，

$λ$

經kernel加權平滑後，原始的KNN迴歸擬合的曲線爲連續的了，以下圖：

可是，這種kernel迴歸一樣存在着邊界（boundary）問題，以下圖：

4. 基於最大似然估計的kernel加權平滑KNN迴歸

上一小節咱們看到，傳統kernel加權平滑的KNN迴歸，對於x序列的開始與結束區段的點，其左右鄰域是不對稱的，致使了平滑後的值偏大或偏小。

所以，須要對權值作再修正，修正的方法就是再次基於窗口內的數據點極大似然估計，使迴歸值x0進一步向均值線μ靠攏。

假定對統計窗口的迴歸量 $x_{0}$

$x_{0}$

定義目標函數：

令：

，

那麼，目標函數可改寫爲：

求偏導，可獲得：

利用矩陣運算，可直接得極大似然估計值爲：

其中，

上述迴歸方法稱之爲LOWESS (Local Weighted regression)。

5. Robust LOWESS

咱們沿着LOWESS的思路脈絡繼續討論，截止目前爲止，LOWESS彷佛圓滿地完成了全部任務，可是還有一個問題，就是樣本數據中的異常點（outlier）問題。

咱們知道，樣本數據中的方差波動來自兩個方面：

系統性偏差：因爲採樣和測量手段引發的偏差，不可避免
樣本自身偏差：因爲樣本數據自身存在的固有方差波動，致使的離羣點

樣本數據集的整體分佈由以上兩種方差波動共同做用產生，LOWESS經過kernel加權和極大似然迴歸的方式已經起到了較好的均值迴歸的做用，可是在極端的離羣點（outlier）狀況下，LOWESS的迴歸點x0依然可能偏離真值均值點μ過遠。所以，人們又提出了Robust LOWESS方法。

Robust LOWESS是Cleveland在LOWESS基礎上提出來的robust迴歸方法，能避免outlier對迴歸的影響。

在計算完估計值後，繼續計算殘差：

根據殘差計算robustnest weight：

其中，s爲殘差絕對值序列的中位值（median），B函數爲bisquare函數：

$s$

而後，用robustness weight乘以kernel weight做爲的新weight。如此，便剔除了殘差較大的異常點對於迴歸的影響，這本質上是數理統計中隨機變量標準化的一種思想應用。

0x3：LOWESS對sin正弦週期數據的擬合效果測試

爲了驗證LOWESS的實際性能，咱們來構造一組數據集驗證其效果，該數據集的特性知足：

週期性：數據集不是單純的線性迴歸關係，而是包含了週期波動性，咱們基於sin函數生成y值
隨機性：數據在不一樣的局部窗口中的密度分佈是不一樣的，所以在不一樣局部窗口之間，均值並非穩定連續變化的

lowess.py庫函數

"""
lowess: Locally linear regression
==================================

Implementation of the LOWESS algorithm in n dimensions.

References
=========
[HTF] Hastie, Tibshirani and Friedman (2008). The Elements of Statistical
Learning; Chapter 6

[Cleveland79] Cleveland (1979). Robust Locally Weighted Regression and Smoothing
Scatterplots. J American Statistical Association, 74: 829-836.

"""
import numpy as np
import numpy.linalg as la


# Kernel functions:
def epanechnikov(xx, **kwargs):
    """
    The Epanechnikov kernel estimated for xx values at indices idx (zero
    elsewhere)

    Parameters
    ----------
    xx: float array
        Values of the function on which the kernel is computed. Typically,
        these are Euclidean distances from some point x0 (see do_kernel)

    Notes
    -----
    This is equation 6.4 in HTF chapter 6

    """
    l = kwargs.get('l', 1.0)
    ans = np.zeros(xx.shape)
    xx_norm = xx / l
    idx = np.where(xx_norm <= 1)
    ans[idx] = 0.75 * (1 - xx_norm[idx] ** 2)
    return ans


def tri_cube(xx, **kwargs):
    """
    The tri-cube kernel estimated for xx values at indices idx (zero
    elsewhere)

    Parameters
    ----------
    xx: float array
        Values of the function on which the kernel is computed. Typically,
        these are Euclidean distances from some point x0 (see do_kernel)

    idx: tuple
        An indexing tuple pointing to the coordinates in xx for which the
        kernel value is estimated. Default: None (all points are used!)

    Notes
    -----
    This is equation 6.6 in HTF chapter 6
    """
    ans = np.zeros(xx.shape)
    idx = np.where(xx <= 1)
    ans[idx] = (1 - np.abs(xx[idx]) ** 3) ** 3
    return ans


def bi_square(xx, **kwargs):
    """
    The bi-square weight function calculated over values of xx

    Parameters
    ----------
    xx: float array

    Notes
    -----
    This is the first equation on page 831 of [Cleveland79].
    """
    ans = np.zeros(xx.shape)
    idx = np.where(xx < 1)
    ans[idx] = (1 - xx[idx] ** 2) ** 2
    return ans


def do_kernel(x0, x, l=1.0, kernel=epanechnikov):
    """
    Calculate a kernel function on x in the neighborhood of x0

    Parameters
    ----------
    x: float array
       All values of x
    x0: float
       The value of x around which we evaluate the kernel
    l: float or float array (with shape = x.shape)
       Width parameter (metric window size)
    kernel: callable
        A kernel function. Any function with signature: `func(xx)`
    """
    # xx is the norm of x-x0. Note that we broadcast on the second axis for the
    # nd case and then sum on the first to get the norm in each value of x:
    xx = np.sum(np.sqrt(np.power(x - x0[:, np.newaxis], 2)), 0)
    return kernel(xx, l=l)


def lowess(x, y, x0, deg=1, kernel=epanechnikov, l=1, robust=False, ):
    """
    Locally smoothed regression with the LOWESS algorithm.

    Parameters
    ----------
    x: float n-d array
       Values of x for which f(x) is known (e.g. measured). The shape of this
       is (n, j), where n is the number the dimensions and j is the
       number of distinct coordinates sampled.

    y: float array
       The known values of f(x) at these points. This has shape (j,)

    x0: float or float array.
        Values of x for which we estimate the value of f(x). This is either a
        single scalar value (only possible for the 1d case, in which case f(x0)
        is estimated for just that one value of x), or an array of shape (n, k).

    deg: int
        The degree of smoothing functions. 0 is locally constant, 1 is locally
        linear, etc. Default: 1.

    kernel: callable
        A kernel function. {'epanechnikov', 'tri_cube', 'bi_square'}

    l: float or float array with shape = x.shape
       The metric window size for the kernel

    robust: bool
        Whether to apply the robustification procedure from [Cleveland79], page
        831


    Returns
    -------
    The function estimated at x0.

    Notes
    -----
    The solution to this problem is given by equation 6.8 in Hastie
    Tibshirani and Friedman (2008). The Elements of Statistical Learning
    (Chapter 6).

    Example
    -------
    >>> import lowess as lo
    >>> import numpy as np

    # For the 1D case:
    >>> x = np.random.randn(100)
    >>> f = np.cos(x) + 0.2 * np.random.randn(100)
    >>> x0 = np.linspace(-1, 1, 10)
    >>> f_hat = lo.lowess(x, f, x0)
    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(1)
    >>> ax.scatter(x, f)
    >>> ax.plot(x0, f_hat, 'ro')
    >>> plt.show()

    # 2D case (and more...)
    >>> x = np.random.randn(2, 100)
    >>> f = -1 * np.sin(x[0]) + 0.5 * np.cos(x[1]) + 0.2*np.random.randn(100)
    >>> x0 = np.mgrid[-1:1:.1, -1:1:.1]
    >>> x0 = np.vstack([x0[0].ravel(), x0[1].ravel()])
    >>> f_hat = lo.lowess(x, f, x0, kernel=lo.tri_cube)
    >>> from mpl_toolkits.mplot3d import Axes3D
    >>> fig = plt.figure()
    >>> ax = fig.add_subplot(111, projection='3d')
    >>> ax.scatter(x[0], x[1], f)
    >>> ax.scatter(x0[0], x0[1], f_hat, color='r')
    >>> plt.show()
    """
    if robust:
        # We use the procedure described in Cleveland1979
        # Start by calling this function with robust set to false and the x0
        # input being equal to the x input:
        y_est = lowess(x, y, x, kernel=epanechnikov, l=l, robust=False)
        resid = y_est - y
        median_resid = np.nanmedian(np.abs(resid))
        # Calculate the bi-cube function on the residuals for robustness
        # weights:
        robustness_weights = bi_square(resid / (6 * median_resid))

    # For the case where x0 is provided as a scalar:
    if not np.iterable(x0):
        x0 = np.asarray([x0])

    #print "x0: ", x0
    ans = np.zeros(x0.shape[-1])
    # We only need one design matrix for fitting:
    B = [np.ones(x.shape[-1])]
    for d in range(1, deg + 1):
        B.append(x ** deg)

    #print "B: ", B

    B = np.vstack(B).T
    for idx, this_x0 in enumerate(x0.T):
        # This is necessary in the 1d case (?):
        if not np.iterable(this_x0):
            this_x0 = np.asarray([this_x0])
        # Different weighting kernel for each x0:
        W = np.diag(do_kernel(this_x0, x, l=l, kernel=kernel))
        # XXX It should be possible to calculate W outside the loop, if x0 and
        # x are both sampled in some regular fashion (that is, if W is the same
        # matrix in each iteration). That should save time.

        print "W: ", W
        if robust:
            # We apply the robustness weights to the weighted least-squares
            # procedure:
            robustness_weights[np.isnan(robustness_weights)] = 0
            W = np.dot(W, np.diag(robustness_weights))
        # try:
        # Equation 6.8 in HTF:
        BtWB = np.dot(np.dot(B.T, W), B)
        BtW = np.dot(B.T, W)
        # Get the params:
        beta = np.dot(np.dot(la.pinv(BtWB), BtW), y.T)
        # We create a design matrix for this coordinat for back-predicting:
        B0 = [1]
        for d in range(1, deg + 1):
            B0 = np.hstack([B0, this_x0 ** deg])
        B0 = np.vstack(B0).T
        # Estimate the answer based on the parameters:
        ans[idx] += np.dot(B0, beta)
    # If we are trying to sample far away from where the function is
    # defined, we will be trying to invert a singular matrix. In that case,
    # the regression should not work for you and you should get a nan:
    # except la.LinAlgError :
    #    ans[idx] += np.nan
    return ans.T

test_lowess.py

# -*- coding: utf-8 -*-

import numpy as np
import numpy.testing as npt
import lowess as lo
import matplotlib.pyplot as plt

def test_lowess():
    """
    Test 1-d local linear regression with lowess
    """
    np.random.seed(1984)
    for kernel in [lo.epanechnikov, lo.tri_cube]:
        for robust in [True, False]:
            # constructing a batch of periodic and volatile data
            x = np.random.randn(100)
            f = np.sin(x)
            x0 = np.linspace(-1, 1, 10)     # x0 means fixed length window width
            f_hat = lo.lowess(x, f, x0, kernel=kernel, l=1.0, robust=robust)
            #print f_hat
            f_real = np.sin(x0)

            plt.figure()
            plt.plot(x0, f_hat, label="lowess fitting curve", color='coral')
            plt.plot(x0, f_real, label="real sin curve", color='blue')

            plt.show()


if __name__ == '__main__':
    np.random.seed(1984)
    kernel = lo.epanechnikov
    robust = True

    # constructing a batch of periodic and volatile data
    x = np.random.randn(100)
    f = np.sin(x)
    x0 = np.linspace(-1, 1, 10)  # x0 means fixed length window width
    f_hat = lo.lowess(x, f, x0, kernel=kernel, l=1.0, robust=robust)
    # print f_hat
    f_real = np.sin(x0)

    plt.figure()
    plt.plot(x0, f_hat, label="lowess fitting curve", color='coral')
    plt.plot(x0, f_real, label="real sin curve", color='blue')

    plt.show()

能夠看到，lowess的擬合迴歸線（橙紅色）相對比較好地擬合了原始數據集的。

0x4：LOWESS和深度神經網絡在思想上的同源與聯繫

這個章節，筆者將嘗試把LOWESS和深度神經網絡進行類比，經過類比，咱們能夠更加深刻理解LOWESS的思想本質。

1. 如何解決非線性問題

在這點上，深度神經網絡用多元神經元疊加的方式，並經過分配不一樣的權重比來獲得一個複雜非線性擬合；而LOWESS的解決思路是縮小統計窗口，在小範圍內的窗口內進行迴歸分析，這也是微積分思想的一種體現。

從模型能力上來講，深度神經網絡更適合擬合複雜非線性大網絡，而LOWESS更適合擬合週期性波動數據。

2. 如何解決離羣異常點問題

在解決這個問題上，LOWESS和深度神經網絡異曲同工，都採用了非線性函數進行權重重分配與調整。在深度神經網絡中，這被稱之爲激活函數，在LOWESS中則被稱之爲kernel加權函數。

實際上，若是咱們將「Epanechnikov二次kernel」和「sigmoid函數的導數」圖像分別打印出來，會發現其函數形態都相似於下圖：

因此它們發揮的做用也是相似的，即：對遠離均值μ迴歸線的離羣點，賦值更小的權重；對靠近均值μ的數據點，賦值更大的權重。

從機率論的角度來講，這種非線性函數的加入本質上也是一種先驗假設，先驗假設的加入，使得數據驅動的後驗估計更加準確，這也是貝葉斯估計思想的核心。

另外一方面來講，深度神經網絡和深度學習的發展，使得人們能夠用深度神經網絡構建出一個超複雜的多元非線性網絡，對目標數據進行擬合，並經過剪枝/dropout等方式緩解過擬合的發生。

Relevant Link:

https://cosx.org/2008/11/lowess-to-explore-bivariate-correlation-by-yihui/
https://cosx.org/tags/lowess
https://github.com/arokem/lowess
https://blog.csdn.net/longgb123/article/details/79520982
https://www.jianshu.com/p/0002f27e80f9?utm_campaign=maleskine&utm_content=note&utm_medium=seo_notes&utm_source=recommendation
https://cloud.tencent.com/developer/article/1111722 
https://www.cnblogs.com/en-heng/p/7382979.html

0x5：LOWESS用於預測模型

首先要注意的是，和KNN同樣，LOWESS是一種非參數（nonparametric）迴歸模型。模型訓練結束後並不能獲得一組肯定的模型參數（coefficient）。

究其本質緣由是由於KNN/LOWESS在訓練並無預先假定一種固定的模型函數公式，相反，KNN/LOWESS直接基於原始數據獲得一條多項式擬合曲線，因此不少時候，咱們也能夠稱其爲平滑算法。

非參數迴歸具有不少良好的特性，例如：

關於兩個變量（X，Y）的關係探索是開放式的，不套用任何現成的數學函數，即不作任何先驗假設
所擬合的曲線能夠很好地描述變量之間的細微變化、局部變化、週期性特徵等信息
非參數迴歸提供的是萬能的擬合曲線，無論多麼複雜的曲線關係都能進行成功的擬合。

基於這個理論前提，非參數模型是沒法對未知輸入進行預測的，可是能夠對訓練集範圍內的x值進行預測。

咱們經過分析KNN的預測能力來解釋這句話，KNN模型本質上是一個非參數模型，但其又能夠進行預測，這並不矛盾，KNN model存儲的本質上是一個很是龐大的參數，也能夠說是將整個訓練集做爲參數存儲起來了。在預測時，KNN.predict要求輸入值x必須在訓練集範圍以內，而KNN所謂的最近鄰搜索其本質上和插值法在思想上是一致的。

回到LOWESS的預測任務上，LOWESS.predict要求輸入的x值必須在訓練數據範圍內。

Relevant Link:

http://webcache.googleusercontent.com/search?q=cache:EgbwJKBCol8J:landcareweb.com/questions/44622/shi-yong-ju-bu-jia-quan-hui-gui-loess-lowess-yu-ce-xin-shu-ju+&cd=2&hl=zh-CN&ct=clnk

2. 時間序列分解算法（STL，Seasonal-Trend decomposition procedure based on Loess）

0x1：時域數據中週期性信號的分解本質是信號從時域到頻域的轉換

傅立葉原理代表：任何連續測量的時序或信號，均可以表示爲不一樣頻率的正弦波信號的無限疊加。而根據該原理創立的傅立葉變換算法利用直接測量到的原始信號，以累加方式來計算該信號中不一樣正弦波信號的頻率、振幅和相位。

以下圖所示，即爲時域信號與不一樣頻率的正弦波信號的關係：

上圖中最右側展現的是時域中的一個信號，這是一個近似於矩形的波。而圖的正中間則是組成該信號的各個頻率的正弦波。

從圖中咱們能夠看出，即便角度幾乎爲直角的正弦波，其實也是由衆多的弧度圓滑的正弦波來組成的。

在時域圖像中，咱們看到的只有一個矩形波，但當咱們經過傅里葉變換將該矩形波轉換到頻域以後，咱們可以很清楚的看到許多脈衝，其中頻域圖中的橫軸爲頻率，縱軸爲振幅。所以能夠經過這個頻域圖像得知，時域中的矩形波是由這麼多頻率的正弦波疊加而成的。

時序信號本質上也是由不少衝擊信號疊加而成的，從時域上，咱們看到的是一條不規整的，上下起伏的曲線，可是經過STL分解，咱們能夠將其分解爲不一樣的信號組成部分。STL的本質就是信號的成分分解，這是STL的核心。

0x2：STL主要解決什麼問題

STL（Seasonal and Trend decomposition using Loess）是一種分解時間序列信號的方法，由 Cleveland et al. (1990) 提出。

咱們下面解釋STL算法主要解決的兩個問題，也是STL的主要應用場景。

1. 信號分解與提取

時間序列分解算法主要被用來測定及提取時間序列信號中的趨勢性和季節性信號，其分解式爲：

，其中

Tt：表示趨勢份量（trend component）
St：表示週期份量（seasonal component）
Rt：表示餘項（remainder component）

因此本質上，STL是一個信號分解算法。

2. 時序信號趨勢性和季節性統計量的定義

基於STL信號分解公式，咱們可使用機率論與數理統計裏的方法論對時序信號創建統計量，分別是趨勢性強度和季節性強度，它們都以隨機變量機率的形式定義。

這兩個度量方法很是實用，當咱們須要在不少時間序列中找到趨勢性或季節性最強的序列，就能夠基於這兩個統計量進行篩選。

1）趨勢性強度

對於趨勢性很強的數據，經季節調整後的數據應比殘差項的變更幅度更大，所以會相對較小。可是，對於沒有趨勢或是趨勢很弱的時間序列，兩個方差應大體相同。所以，咱們將趨勢強度定義爲：

這能夠給趨勢強度的衡量標準，其值在0-1之間。由於有些狀況下殘差項的方差甚至比季節變換後的序列還大，咱們令 $F_{T}$

$F_{T}$

季節性的強度定義以下

當季節強度 $F_{S}$

0x3：STL算法總體框架性架構

STL的總體框架性結構以下

筆者會在這個章節先進行框架性和原理性的討論，不涉及具體的算法細節，某個部分也會進行簡化處理，目的是總總體上對STL創建一個感性的認識。

1. 數據拆分規則

數據拆分規則是說咱們採用什麼模型來表徵時序數據的疊加方式。

1）加法方式

原始數據 = 平均季節數據 + 趨勢數據 + 殘差

這種方式，隨着時間的推移季節數據不會有太大的變化，在週期的業務數據更適合這樣的拆分方式，咱們在工業場景中用的最多的也是這種分解方式。

加法方式也是傅里葉變換時域信號分解的信號疊加方式。

2）乘法方式

原始數據 = 季節數據 * 趨勢數據 * 殘差

這種分解方式，隨着時間的推移季節數據波動會很是明顯。

2. 肯定數據週期

要從時序信號中提取出週期性信號，首要的事情是須要肯定時序週期，即多久循環一次的時間窗口長度。時序週期的肯定大概有兩個方法：

根據業務經驗人工設定：例如像博客園訪客趨勢這種業務場景，時序週期顯然應該設定爲7天。大部分項目中，咱們都是採用人工設定的方式
基於傅里葉變換技術從樣本數據中自動提取

3. 提取趨勢數據

從時序數據中抽取趨勢信息，最合理的方法就是尋找均值迴歸線，均值迴歸線反映了時序數據總體的趨勢走向。

上圖是分解後造成了：趨勢數據 + 其餘數據。
下圖是去到其餘數據最終的趨勢數據。

在STL中，趨勢信息的提取是經過LOWESS算法實現的。

4. 提取季節週期性數據

獲得了趨勢數據以後，接下來才能夠計算季節週期性數據。計算公式以下：

原始數據 - 趨勢數據 = 季節數據 + 殘差數據 = S

接下來經過週期均值化從S中提取季節週期性數據。

數據集S中對不一樣週期中，相同時間點數據的進行相加獲得數據集SP，而後除以原始數據的週期數c，獲得一個平均季節數據集H，而後複製c次數據集H，將獲得平均季節數據，以下圖：

5. 提取殘差數據

殘差數據的獲取公式以下：

原始數據 - 趨勢數據 - 平均季節數據 = 殘差數據

最終分解爲不一樣的信號部分的結果以下：

0x4：STL分解算法流程細節

首先，STL是一個迭代循環，整體思想有點相似EM算法的分佈優化過程，經過多輪迭代，逐步逼近「趨勢信號」和「週期信號」的機率分佈真值。

STL分爲內循環（inner loop）與外循環（outer loop），其中內循環主要作了趨勢擬合與週期份量的計算。外層循環主要用於調節robustness weight，即去除outlier離羣點。

1. STL內循環

下圖爲STL內循環的流程圖：

1）初始化

內循環從趨勢份量開始，初始化。

初始化完成後內循環。

2）Step-1：去趨勢份量

將本輪的總信號量，減去上一輪結果的趨勢份量，

3）Step-2：去除季節週期性信號

去除完趨勢份量以後，餘下的信號主要由兩部分組成，即」殘差信號+季節週期性信號「，接下來要作的事是提取出季節週期性信號。咱們分幾小步討論。

3.1）週期子序列平滑（Cycle-subseries smoothing）

將每一個週期間相同位置的樣本點（週期平移後相同位置）組成一個子序列集合（subseries），容易知道這樣的子序列個數等於Np個（一個週期內的樣本數爲Np） $n_{(} p)$

$n_{(} p)$

3.2）週期子序列的低通量過濾（Low-Pass Filtering）- 提取週期子序列中的非週期性信號

對上一個步驟的」temporary seasonal series」 $C_{v}^{(k + 1)}$

$C_{v}^{(k + 1)}$ $C_{v}^{(k + 1)}$

3.3）去除平滑週期子序列趨勢（Detrending of Smoothed Cycle-subseries）- 週期子序列週期信號提純

提純後獲得季節週期性信號。

3.4）去除週期性信號

4）Step-3：趨勢平滑（Trend Smoothing）- 提取下一輪迭代中須要用到的趨勢信號

對於去除週期以後的序列作 $q = n_{n (t)}$

以後須要判斷是否收斂，若是未收斂則繼續下一輪迭代，若是收斂則輸出本輪迭代的最終「趨勢份量結果」與「週期份量結果」。

$q = n_{n (t)}$

外層循環主要用於調節LOESS算法的robustness weight。若是數據序列中有outlier，則餘項會較大。

具體細節能夠參閱文章前面對LOWESS的robustness weight的相關討論。

爲了使得算法具備足夠的robustness，因此設計了內循環與外循環。特別地，當內存循環數足夠大時，內循環結束時趨勢份量與週期份量已收斂；

同時，若時序數據中沒有明顯的outlier，能夠將外層循環數設爲0。

0x5：STL信號分解代碼示例

# -*- coding: utf-8 -*-

import statsmodels.api as sm
import matplotlib.pyplot as plt
import pandas as pd

from datetime import datetime
from heapq import nlargest
from re import match
import pytz
import numpy as np

def datetimes_from_ts(column):
    return column.map(
        lambda datestring: datetime.fromtimestamp(int(datestring), tz=pytz.utc))

def date_format(column, format):
    return column.map(lambda datestring: datetime.strptime(datestring, format))

def format_timestamp(indf, index=0):
    if indf.dtypes[0].type is np.datetime64:
        return indf

    column = indf.iloc[:,index]

    if match("^\\d{4}-\\d{2}-\\d{2} \\d{2}:\\d{2}:\\d{2} \\+\\d{4}$",
             column[0]):
        column = date_format(column, "%Y-%m-%d %H:%M:%S")
    elif match("^\\d{4}-\\d{2}-\\d{2} \\d{2}:\\d{2}:\\d{2}$", column[0]):
        column = date_format(column, "%Y-%m-%d %H:%M:%S")
    elif match("^\\d{4}-\\d{2}-\\d{2} \\d{2}:\\d{2}$", column[0]):
        column = date_format(column, "%Y-%m-%d %H:%M")
    elif match("^\\d{2}/\\d{2}/\\d{2}$", column[0]):
        column = date_format(column, "%m/%d/%y")
    elif match("^\\d{2}/\\d{2}/\\d{4}$", column[0]):
        column = date_format(column, "%Y%m%d")
    elif match("^\\d{4}\\d{2}\\d{2}$", column[0]):
        column = date_format(column, "%Y/%m/%d/%H")
    elif match("^\\d{10}$", column[0]):
        column = datetimes_from_ts(column)

    indf.iloc[:,index] = column

    return indf

def get_gran(tsdf, index=0):
    col = tsdf.iloc[:,index]
    n = len(col)

    largest, second_largest = nlargest(2, col)
    gran = int(round(np.timedelta64(largest - second_largest) / np.timedelta64(1, 's')))

    if gran >= 86400:
        return "day"
    elif gran >= 3600:
        return "hr"
    elif gran >= 60:
        return "min"
    elif gran >= 1:
        return "sec"
    else:
        return "ms"



if __name__ == '__main__':
    # https://github.com/numenta/NAB/blob/master/data/artificialWithAnomaly/art_daily_flatmiddle.csv
    dta = pd.read_csv('./art_daily_flatmiddle.csv', usecols=['timestamp', 'value'])
    dta = format_timestamp(dta)
    dta = dta.set_index('timestamp')
    dta['value'] = dta['value'].apply(pd.to_numeric, errors='ignore')
    dta.value.interpolate(inplace=True)

    print "dta.value: ", dta.value
    res = sm.tsa.seasonal_decompose(dta.value, freq=288)
    res.plot()
    plt.show()

0x6：基於STL信號分解算法的應用

1. STL時序預測

STL時間序列分解算法不只本質上在學習時間序列中和在探索歷史隨之間變化中十分有用，它也能夠應用在預測中。

假設一個加法分解，分解後的時間序列能夠寫爲：

其中，St表明季節週期性信號，At表明趨勢性和殘差信號的混合。

對於時間序列的預測，咱們須要分別預測季節項 ${\hat{S}}_{t}$

${\hat{S}}_{t}$

另外一方面，可使用任意非季節性預測方法來預測份量At。例如，可使用帶漂移項的隨機遊走法，也可使用三次指數平滑法（Holt-Winters法），或者還能夠用非季節性的ARIMA模型。

咱們這裏以雙色球爲例，基於歷史雙色球開獎的記錄做爲時序數據，預測將來可能的雙色球開獎結果。

# -*- coding: utf-8 -*-

import statsmodels.api as sm
import matplotlib.pyplot as plt
import pandas as pd
import os
import numpy as np
import pickle
import pandas
from statsmodels.tsa.api import Holt


def load_historydata():
    # http://www.17500.cn/getData/ssq.TXT
    if not os.path.isfile("./ssq.pkl"):
        ori_data = np.loadtxt('./ssq.TXT', delimiter=' ', usecols=(0, 2, 3, 4, 5, 6, 7, 8), unpack=False)
        pickle.dump(ori_data, open("./ssq.pkl", "w"))
        df = pandas.DataFrame(data=ori_data)
        return df
    else:
        ori_data = pickle.load(open("./ssq.pkl", "r"))
        df = pandas.DataFrame(data=ori_data)
        return df


if __name__ == '__main__':
    dta = load_historydata()
    #print dta

    dta = dta.set_index(0)
    # Interpolate values
    for i in [1, 2, 3, 4, 5, 6, 7]:
        dta[i].interpolate(inplace=True)

    print dta

    # show the stl decomposition result
    for i in [1, 2, 3, 4, 5, 6, 7]:
        res = sm.tsa.seasonal_decompose(dta[i], freq=7)
        res.plot()
    plt.show()


    # predict the future time series data
    y_pred = []
    for j in range(1, 8):
        y_hat = []
        for i in [1, 2, 3, 4, 5, 6, 7]:
            fit = Holt(np.asarray(dta[i])).fit(smoothing_level=0.3, smoothing_slope=0.1)
            y_hat.append(int(fit.forecast(j)[-1]))
        y_pred.append(y_hat)

    print "y_pred: ", y_pred

從STL分解結果殘差信號能夠看出，雙色球的開獎結果的具有很大的隨機性，並無呈現出明顯的季節週期性和趨勢性。

這裏預測了8天的將來開獎結果以下：

[
    [5, 14, 18, 22, 26, 30, 8], 
    [5, 14, 18, 22, 26, 30, 8], 
    [5, 15, 18, 22, 26, 30, 8], 
    [5, 15, 18, 23, 26, 30, 8], 
    [5, 15, 18, 23, 26, 30, 7], 
    [5, 16, 19, 23, 26, 30, 7], 
    [5, 16, 19, 23, 26, 30, 7]
]

Relevant Link:

https://otexts.com/fppcn/forecasting-decomposition.html
http://www.sunsoda.fun/TS%E6%97%B6%E9%97%B4%E5%88%86%E6%9E%90(3).html

2. 時序信號異常檢測

從時序信號中提取出殘差數據後，咱們能夠基於此作時序異常檢測，時序數據若是有異常，都會體如今殘差數據集中。下面列舉幾種對殘差信號進行異常檢測的方法：

1）殘差信號假設檢驗

假設殘差數據知足正太分佈，或者近似正太分佈，能夠計算出殘差數據集的標準差。

若是數據點與均值的差值絕對值在3倍標準差外，則認爲是異常點，也就是3sigma異常檢驗，例以下面的殘差信號圖：

2）Generalized Extreme Studentized Deviate (GESD)

GESD，又名Grubb's Test檢驗，是一種多輪迭代發現多個異常點的異常發現算法。關於grubb's test的討論，能夠參閱另外一篇文章。

Relevant Link:

https://github.com/numenta/NAB
https://github.com/andreas-h/pyloess 
http://siye1982.github.io/2018/03/02/anomaly-detection1/
https://www.cnblogs.com/runner-ljt/p/5245080.html 
https://www.cnblogs.com/en-heng/p/7390310.html 
https://zhuanlan.zhihu.com/p/34342025

3. 時序異常檢測的應用思考

0x1：什麼問題場景下適合用時序異常檢測

這裏拉取了筆者在博客園放置的cnzz來訪者統計js插件的流量趨勢統計結果

能夠看到，個人博客園站點訪問呈現明顯的週期波動性，可是仔細觀察也會發現，在這個週期性之上，趨勢線中依然存在幾個「看起來不是那麼正常」的「離羣點」，這些離羣點就是由於筆者發表了一些新文章或者什麼其餘的緣由致使訪問流量發生非週期性的趨勢波動。

拋開具體的實際場景案例，抽象地來講，時間序列的異常檢測問題一般表示爲：在時變信息系統下，相對於某些標準信號或週期規律性信號的離羣點，

時序信號中有不少的異常類型，可是咱們只關注業務角度中最重要的類型，好比意外的峯值、降低、趨勢變化以及等級轉換。

咱們下面經過博客園的訪問趨勢來舉例說明典型的時序異常。

1. 附加異常

舉個例子來講，筆者發現每逢週末博客園訪問量都會降低，而到週一又會回到峯值，可是若是咱們發現某個週末訪問量忽然出現一個異常增長，看起來就像一個峯值。這些類型的異常一般稱爲附加異常。這種狀況極可能是由於筆者在週末發表了一篇新博客。

2. 時間變化異常

關於博客園的另外一個例子是，當某天博客園服務器宕機，咱們會發如今某個時間段內有零個或者很是少的用戶訪問，明顯和平時的趨勢不一樣，這些類型的異常一般被分類爲時間變化異常。

3. 水平位移或季節性水平位移異常

仍是以博客園爲例，咱們假設國家修改了法定假日的時間，從週六/日改到週一/二，那麼能夠想象，總體的訪問趨勢會平移2天，每週的低谷平移到週一/二。

這些類型的變化一般被稱爲水平位移或季節性水平位移異常。

整體來講，異常檢測算法應該將每一個時間點標記爲異常/非異常，或者預測某個點的信號，並衡量這個點的真實值與預測值的差值是否足夠大，從而將其視爲異常。

0x2：基於STL時序分解算法作時序異常檢測的總體思路

這個章節筆者嘗試從具體應用場景的角度出發，討論下怎麼基於STL分解算法作時序異常檢測。

1. 分析數據源是否具有時序週期性

包括STL在內，不少時序算法都有一個默認的假設前提，即：數據源自己是由週期性的衝擊信號源產生的，即數據自己具有季節週期性。

典型的正例子包括網站訪問pv量，通常來講是存在【天單位週期性】，由於通常來講，網站的pv訪問量天天是呈現相對固定的起伏波動的。

典型的反例包括每期的雙色球開獎結果，其數據自己並不存在明顯的週期性，因此STL沒法從其中提取出任何有機率顯著意義的週期性信號，天然也沒法基於殘差信號作異常檢測了。

從筆者本身的經驗來看，這個假設前提就極大限制了時序算法在安全領域裏的應用，從這個角度來分析，時序異常算法可能只適合在網站cc，服務器ddos，異常運維操做，kpi異常檢測等場景應用。

2. 肯定時序週期

確認了問題場景的數據適合於時序算法的假設前提後，接下來的問題就是：數據的週期性窗口多大。

例如網站pv訪問量的週期能夠設置爲1天，外賣銷量時序時序數據的週期性能夠設置爲7天，某地區降雨量時序數據的週期性能夠設置爲1年，等等等等。

一個好的建議是基於你本身的領域經驗來設置這個週期窗口。

3. 基於STL分解獲得的殘差信號進行異常離羣點檢測

對於時序數據來講，趨勢信號和週期信號是須要過濾去除的部分，餘下的殘差信號纔是作異常檢測的數據源。

對殘差信號的異常檢測是以週期窗口爲單位的，在每個定長的週期窗口內，對時序信號的指標進行異常檢測（假設檢驗或無監督異常檢測）

4. 基於STL時序預測的異常檢測

能夠基於對歷史數據的擬合，對將來的時序指標進行預測，將預測的結果做爲一個基線。當將來的時序指標發生時，將其和基線進行對比，經過度量偏離程度，以此判斷是不是異常離羣點。

0x3：基於STL時序分解算法的時序異常檢測在SSH異常登陸行爲檢測裏的應用

筆者這裏拉取了某臺網站服務器一段時間內的ssh登陸時序流水記錄，經過STL進行了分解，x軸爲分鐘時序，縱軸爲每分鐘內登陸次數。

週期窗口配置爲60min

若是咱們將週期窗口拉長，時序數據中的週期性強度會下降，這本質上和STL的算法迭代過程有關，越長的週期窗口，週期信號被滑動平均的稀疏就越嚴重。

週期窗口配置爲120mi

在STL的基礎上，筆者對殘差信號進行grubb's test假設檢驗異常檢測。可是實驗的結果並不理想。

這也反過來講明瞭ssh/rdp登陸日誌並不具有明顯的時序週期性。

更通常地說，時序異常檢測要求待檢測的數據源自己就具有「時間密集性」，即不發生異常時，數據源自己就是一個包含週期性信號的密集時序信號，例如網站平常pv訪問，4層數據包訪問。

而時序中的異常定義，是說：某個時間點突發的忽然性時序指標突變，或者總體時序週期性的水平平移異常（這種比較多），在業務場景中，佔大多數的是前者，即某時刻時序指標的突變異常。

拿服務器的ssh/rdp登陸日誌來講，首先它是一個非時間密集型的時序數據，由於平常除了管理員以外，不多會有登陸行爲，若是一旦有大量登陸行爲，極可能就是有攻擊者在對服務器進行爆破攻擊。若是直接強行對ssh/rdp登陸日誌進行stl分解，獲得的殘差信號極可能就是攻擊者的流量數據，基於這份殘差數據對假設檢驗獲得的異常時間點，自己也沒有可解釋性，極可能只是當時攻擊者增長了攻擊機器或者增長了網絡帶寬而已。

因此，ssh/rdp登陸時序檢測這個場景自己就不符合時序信號的基本假設，。

Relevant Link:

https://juejin.im/post/5c19f4cb518825678a7bad4c
https://github.com/twitter/AnomalyDetection
https://mp.weixin.qq.com/s/w7SbAHxZsmHqFtTG8ZAXNg

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。

關於局部加權迴歸（Locally Weighted Scatterplot Smoothing，LOWESS）、STL（Seasonal-Trend decomposition procedure b

1. LOWESS（Locally Weighted Scatterplot Smoothing，局部加權迴歸）

0x1：lowess算法主要解決什麼問題

1. 非線性迴歸擬合問題

2. 局部平滑問題

0x2：從KNN線性迴歸逐漸演進爲LOWESS

1. 多元非線性迴歸模型

2. KNN線性迴歸

3. 引入kernel加權平滑的KNN線性迴歸

4. 基於最大似然估計的kernel加權平滑KNN迴歸

5. Robust LOWESS

0x3：LOWESS對sin正弦週期數據的擬合效果測試

0x4：LOWESS和深度神經網絡在思想上的同源與聯繫

1. 如何解決非線性問題

2. 如何解決離羣異常點問題

0x5：LOWESS用於預測模型

2. 時間序列分解算法（STL，Seasonal-Trend decomposition procedure based on Loess）

0x1：時域數據中週期性信號的分解本質是信號從時域到頻域的轉換

0x2：STL主要解決什麼問題

1. 信號分解與提取

2. 時序信號趨勢性和季節性統計量的定義

1）趨勢性強度

2）季節性強度

0x3：STL算法總體框架性架構

1. 數據拆分規則

1）加法方式

2）乘法方式

2. 肯定數據週期

3. 提取趨勢數據

4. 提取季節週期性數據

5. 提取殘差數據

0x4：STL分解算法流程細節

1. STL內循環

1）初始化

2）Step-1：去趨勢份量

3）Step-2：去除季節週期性信號

3.1）週期子序列平滑（Cycle-subseries smoothing）

3.2）週期子序列的低通量過濾（Low-Pass Filtering）- 提取週期子序列中的非週期性信號

3.3）去除平滑週期子序列趨勢（Detrending of Smoothed Cycle-subseries）- 週期子序列週期信號提純

3.4）去除週期性信號

4）Step-3：趨勢平滑（Trend Smoothing）- 提取下一輪迭代中須要用到的趨勢信號

2. STL外循環

0x5：STL信號分解代碼示例

0x6：基於STL信號分解算法的應用

1. STL時序預測

2. 時序信號異常檢測

1）殘差信號假設檢驗

2）Generalized Extreme Studentized Deviate (GESD)

3. 時序異常檢測的應用思考

0x1：什麼問題場景下適合用時序異常檢測

1. 附加異常

2. 時間變化異常

3. 水平位移或季節性水平位移異常

0x2：基於STL時序分解算法作時序異常檢測的總體思路

1. 分析數據源是否具有時序週期性

2. 肯定時序週期

3. 基於STL分解獲得的殘差信號進行異常離羣點檢測

4. 基於STL時序預測的異常檢測

0x3：基於STL時序分解算法的時序異常檢測在SSH異常登陸行爲檢測裏的應用

$F_{T}$

$q = n_{n (t)}$