LOWESS

LOWESS#

散布図に従う近似線を描くために開発された局所回帰。 LOESS (locally estimated scatterplot smoothing) や LOWESS (locally weighted scatterplot smoothing) と呼ばれる。

Hide code cell source 
import numpy as np
np.random.seed(0)
x = np.linspace(0, 7, 100)
y = np.sin(x) + np.random.normal(0, 0.2, 100)

import statsmodels.api as sm
smoothed = sm.nonparametric.lowess(exog=x, endog=y, frac=0.2, it=3)

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.scatter(x, y)
ax.plot(smoothed[:, 0], smoothed[:, 1], c="k")
ax.set(xlabel="x", ylabel="y")
fig.show()
../_images/ca54b0c10482b6625d1867633bc7f1cce9874b7abbd5873a7855dd37226f49e9.png

LOWESSのアルゴリズム#

Cleveland (1979). に記されたアルゴリズムは以下の通り。

i番目のサンプルの目的変数yiを特徴量xiとノンパラメトリックの平滑化関数g(xi)で近似することを考える。

yi=g(xi)+ϵi

ここでϵiは平均0で分散が一定の確率変数である。

1. 重みの計算とr個の最近傍サンプルの取得#

xiについて、j=1,…,nにわたって|xi−xj|で距離を測り、r番目に近いサンプルとの距離をhiとする。

重み関数W(⋅)を用いて、k=1,…,nについて

wk(xi)=W(hi−1(xk−xi))

を計算する

ここで重み関数W(⋅)は以下の性質を満たすものとする

  1. |x|<1について W(x)>0

  2. W(−x)=W(x)

  3. W(x) は x≥0について 非増加関数

  4. |x|≥1についてW(x)=0

hi−1(xk−xi)は分子のxk−xiの絶対値が分母のhiより大きければ|hi−1(xk−xi)|≥1になるので重みが0になる。つまり、サンプルとして回帰に使用されなくなる。 なので重み関数は近傍のr個のサンプルを取り出しつつ、r個のサンプルにも距離に応じた重みをかける操作となる。

Wの例として tricube functionが考えられる

W(x)={(1−|x|3)3 for |x|<10 for |x|⩾1
Hide code cell source 
def tricube(x: float) -> float:
    if np.abs(x) >= 1:
        return 0
    return (1 - np.abs(x)**3 )**3

x = np.linspace(-2, 2, 100)
w = [tricube(x_i) for x_i in x]

# 別のやり方
# w = (1 - np.abs(x)**3)**3
# w[np.abs(x) >= 1] = 0
fig, ax = plt.subplots(figsize=[4,3])
ax.plot(x, w)
ax.set(title="tricube weight function", xlabel="x", ylabel="W(x)")
fig.show()
../_images/75f9a72542007636aa648af7a7540d10a6086b34a2518287debe5525e8be7eaf.png

2. 多項式回帰のフィッティング#

非線形回帰としてd次の多項式回帰を行う

minβ0,…,βd ∑k=1nwk(xi)(yk−β0−β1xk−…−βdxkd)2
y^i=∑j=0dβ^j(xi)xij

3. ロバスト性重みδの計算#

続いて、外れ値の影響を除外するための重みを計算する。 bisquare weight function B(x)を以下のように定義する

B(x)={(1−x2)2 for |x|<10 for |x|⩾1

残差ei=yi−y^iの絶対値|ei|の中央値をsとする。ロバスト性重み（robustness weights）を

δk=B(ek/6s)

と定義する。B(x)もtricube functionと似た形状であり、残差の絶対値の中央値の6倍（6s）以上の絶対値の残差|ek/6s|≥1を持つ外れ値は重みδkがゼロになり、推定に含まれなくなるので、推定からハズレ値の影響を除外できる。

Hide code cell source 
def bisquare(x: float) -> float:
    if np.abs(x) >= 1:
        return 0
    return (1 - x**2 )**2

x = np.linspace(-2, 2, 100)
b = [bisquare(x_i) for x_i in x]
fig, ax = plt.subplots(figsize=[4,3])
ax.plot(x, b)
ax.set(title="bisquare weight function", xlabel="x", ylabel="B(x)")
fig.show()
../_images/c94a781fbbd34d0da92206e2680f29fc959b62ff02535a889e6d25d2867dc5cc.png

4. δで重み付け回帰を行う#

またd次多項式回帰を行い、新たな推定値y^iを得る。このとき、重みはδkwk(xi)を使う。

5. 繰り返す#

3.と4.のステップをt回繰り返す。

実装#

# サンプルデータ
import numpy as np
n = 100
np.random.seed(0)
x = np.linspace(0, 7, n)
y = np.sin(x) + np.random.normal(0, 0.2, n)

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.scatter(x, y)
ax.set(xlabel="x", ylabel="y")
fig.show()
../_images/480b8e93e5f8bd5a23e8e31ee3f0bf45cac0c246950c154933db067bd4cff6f8.png 
frac = 0.66 # 使用するサンプルの割合
r = int(frac * n)  # 使用する近傍のサンプル数
d = 3  # 多項式回帰の次数
t = 3  # iteration

def tricube(x: np.array) -> np.array:
    w = (1 - np.abs(x)**3)**3
    w[np.abs(x) >= 1] = 0
    return w

def bisquare(x: np.array) -> np.array:
    w = (1 - np.abs(x)**2)**2
    w[np.abs(x) >= 1] = 0
    return w

# d次多項式を作るための特徴量生成
X = np.vstack([x**j for j in range(d)]).T

n = X.shape[0]
delta = np.ones_like(x)
y_pred = np.zeros_like(x)
for _ in range(t):
    for i in range(n):
        # 重みの計算
        dist = x - x[i]
        idx = np.argsort(np.abs(dist))[:r]
        h_i = np.abs(dist[idx]).max() # r番目に近いdiff
        w = tricube(dist / h_i)
        W = np.diag(delta * w)
    
        # WLS
        beta = np.linalg.inv(X.T @ W @ X) @ X.T @ W @ y
        y_pred[i] = X[i,:] @ beta
    
    e = y - y_pred
    s = np.median(np.abs(e))
    delta = bisquare(e / (6 * s))

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.scatter(x, y)
ax.plot(x, y_pred, c="k")
ax.set(xlabel="x", ylabel="y")
fig.show()
../_images/697b76781ba409429044ed1a29b00698a41002881c538944c69fa5defdd06bb2.png

最近のLOWESSアルゴリズム#

Wikipediaには、重み関数がtricubeではなくGaussianを使うものが紹介されている

Gaussian weight functionとは、2つのデータ点の特徴量ベクトルx,x′∈Rm（mは特徴量の次元数）について、

w(x,x′,α)=exp⁡(−‖x−x′‖22α2)

といったもの。

参考#

Contents