季節性の変動要因をフーリエ成分で予測する

10096 ワード
Seasonality Python フーリエ級数 フーリエ変換 statsmodels Python テキストリンク

仰々しいタイトルを掲げていますが、がっつりとしたフーリエ級数の扱いはしません笑

今回はKaggleチュートリアル3章の
Seasonality
をもとに「季節性の変動のような周期性をモデル構築して予測する」ことをテーマとしています。

前回までは純粋な線形回帰のみでモデリングしてきましたが、周期性をうまく捉えて予測したい場合にどうすれば良いのかをみていきます。

今回Fourier成分という数学的知識が出てくるので、高度であり完全な理解というのはなかなか怪しいです。
（理解せずとも、そーなんだーくらいでいいと思います。）

Seasonal Plotsと季節性特徴量

（Seasonal Plotsの訳がうまく見つからないのでこのままにします。）

前回まではダミー変数や移動平均を使って区間ごとの特徴を把握しようとしてきました。
次は、季節的な周期を見つけようとするために曜日ごとにダミー変数を作成する方法があります

rng = np.random.default_rng()
x = np.linspace(100, 200, 100)
y = np.round(x + rng.integers(100, 470, size=100), 0)
dates = pd.date_range(start='2016-1-4', periods=100)
data = pd.DataFrame({'y': y}, index=dates)

data['day']= pd.to_datetime(data.index, format='%W')
data['day'] = [day.day_name() for day in data['day']]
data = pd.get_dummies(data, drop_first=True)
data.head()


(こんな感じ)

one-hotやラベルエンコーディングなどアプローチの方法はさまざまで、この変数を用いて回帰モデルを作成したりすることできます。

公式のチュートリアルから引用すると、以下のようなイメージです

（引用元）: https://www.kaggle.com/ryanholbrook/seasonality#Seasonal-Plots-and-Seasonal-Indicators

各曜日の変数はどれかが1であればその他全て0をとるため、切片+（曜日の変数）*1という値をとります。

フーリエ成分とPeriodogram

上記のアプローチは今までの方法を踏襲していましたが、ここからすこし踏み込みます。

先程のような曜日の指標よりも遥かに長い期間の周期性があるときに、フーリエ成分は役に立ちます。
しかもフーリエ成分は少ない特徴量で連続的な関数を表現できるため、機械学習と相性がよかったりします。

フーリエ級数

「フーリエ、フーリエっていってますけど、なんなんですか！」って方もいるはずなので、ほんとに簡単な公式とその説明だけして進めます。

公式は以下。
$$ f(x) = \frac{a_0}{2} + \sum_{n=1}^{∞}(a_n \cos{nx}+ b_n \sin{nx}) $$
任意の関数はこのf(x)、つまり三角関数の級数で表現できる
というのが、めちゃくちゃ噛み砕きに噛み砕いたフーリエ級数です。

大学数学の範囲になるので詳しい説明はできないですが、関数を周期的な関数（sin, cos）で実は表現できちゃうんだよね、っていう雰囲気でOKです。
（ものすごい余談ですが、これを機にフーリエ変換を数学的にもっと理解したくなったりしました。）

今回の季節性変動のような周期を持ちうる動向を見たいときに、フーリエ級数で表現しようとするアプローチは非常に効果的であるといえます。

ここでのフーリエ成分

区間は任意ですが、例えば1年を一つの区間とすると何回周期があるのか？ということをsin, cosine curveとかで観測したりします。

こちらも再度公式のチュートリアルから引用すると、

（引用元）：https://www.kaggle.com/ryanholbrook/seasonality#Seasonal-Plots-and-Seasonal-Indicators
上記は1年に1周期。下は1年に2周期ある、というイメージです。

これらの三角関数を使って線形回帰のモデルを用いるとそれぞれの関数に対する重みなどを計算してくれるため、将来的な予測もすることが可能となります。

反対に最初に見てきた曜日の指標などをつかったり季節的な指標を使うと、非常に多くの特徴量を作成する必要が出てきます。（過学習などにもつながるので危険）

フーリエ成分の関数

数学的理解を少しでも深めたいという方は公式チュートリアルにフーリエ成分の関数がありますので、少し使ってみましょう。

# https://www.kaggle.com/ryanholbrook/seasonality#Seasonal-Plots-and-Seasonal-Indicators
def fourier_features(index, freq, order):
    time = np.arange(len(index), dtype=np.float32)
    k = 2 * np.pi * (1 / freq) * time
    features = {}
    for i in range(1, order + 1):
        features.update({
            f"sin_{freq}_{i}": np.sin(i * k),
            f"cos_{freq}_{i}": np.cos(i * k),
        })
    return pd.DataFrame(features, index=index)


# 先程のデータ
fourier_features(data.index, freq=365.25, order=4)


関数だけ見ても何がなんやらだとは思いますが、導出としては非常にシンプルです。
（ベースはフーリエ変換だったりするので、なんでこんな式になっているのか？に関しては省略します。）

フーリエ成分を最適に選ぶためのPeriodogram

フーリエ級数の威力が（やんわり）理解したところで、どのくらいの関数を準備すれば良いのか？ということを考えなければいけません。

そこでperiodogramを用います。

y軸はそれぞれの$a_n, b_n$の2乗和の平均$\frac{a_n^2 + b_n^2}{2}$
x軸は各期間の頻度（annual（1年間）, quarterly（4半期）, weekly（1週間ごと）など）をとります


（引用元）https://www.kaggle.com/ryanholbrook/seasonality#Seasonal-Plots-and-Seasonal-Indicators

上記の公式のグラフにy軸はVariance（分散）と書かれていることから、このVarianceが低くなる頻度を探し、関数の数を決めていきます。

データを使って実装の練習

わかりにくい説明だったかもしれませんが、最後に実装して終わりましょう。

関数の準備（公式のものを利用）

# annotations: https://stackoverflow.com/a/49238256/5769929
def seasonal_plot(X, y, period, freq, ax=None):
    if ax is None:
        _, ax = plt.subplots()
    palette = sns.color_palette("husl", n_colors=X[period].nunique(),)
    ax = sns.lineplot(
        x=freq,
        y=y,
        hue=period,
        data=X,
        ci=False,
        ax=ax,
        palette=palette,
        legend=False,
    )
    ax.set_title(f"Seasonal Plot ({period}/{freq})")
    for line, name in zip(ax.lines, X[period].unique()):
        y_ = line.get_ydata()[-1]
        ax.annotate(
            name,
            xy=(1, y_),
            xytext=(6, 0),
            color=line.get_color(),
            xycoords=ax.get_yaxis_transform(),
            textcoords="offset points",
            size=14,
            va="center",
        )
    return ax


def plot_periodogram(ts, detrend='linear', ax=None):
    from scipy.signal import periodogram
    fs = pd.Timedelta("1Y") / pd.Timedelta("1D")
    print('fs:', fs)
    freqencies, spectrum = periodogram(
        ts,
        fs=fs,
        detrend=detrend,
        window="boxcar",
        scaling='spectrum',
    )
    if ax is None:
        _, ax = plt.subplots()
    ax.step(freqencies, spectrum, color="purple")
    ax.set_xscale("log")
    ax.set_xticks([1, 2, 4, 6, 12, 26, 52, 104])
    ax.set_xticklabels(
        [
            "Annual (1)",
            "Semiannual (2)",
            "Quarterly (4)",
            "Bimonthly (6)",
            "Monthly (12)",
            "Biweekly (26)",
            "Weekly (52)",
            "Semiweekly (104)",
        ],
        rotation=30,
    )
    ax.ticklabel_format(axis="y", style="sci", scilimits=(0, 0))
    ax.set_ylabel("Variance")
    ax.set_title("Periodogram")
    return ax

なっっがいですが、seasonal_plotに関しては描画の指定が複雑ですが、あくまでプロットしているだけなので、省略。

したのperiodogramに関して少し深掘りします。
from scipy.signal import periodogramですが、関数としては時系列データを指定し、そこから頻度を指定したりします。
（もろもろのパラメータは公式を参照する方がおすすめ）
scipy.signal.periodogram

実際にデータを自作して当てはめてみるとごちゃつく

plot_periodogram(data.Sales);


理想は公式チュートリアルの以下。

かなり綺麗ですな。。

statsmodelsで三角関数を生成

前回使ったDeterministicProcessでは定数やらTrendのような連番を振ったりしましたが、フーリエ成分を使うことで季節性の要因を追加することができます。

from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess

fourier = CalendarFourier(freq="A", order=10)  # 10 sin/cos pairs for "A"nnual seasonality
# print(fourier)
dp = DeterministicProcess(
    index=data.index,
    constant=True,               # dummy feature for bias (y-intercept)
    order=1,                     # trend (order 1 means linear)
    seasonal=True,               # weekly seasonality (indicators)
    additional_terms=[fourier],  # annual seasonality (fourier)
    drop=True,                   # drop terms to avoid collinearity
)

X = dp.in_sample()  # create features for dates in tunnel.index

予測

y = data["Sales"]

model = LinearRegression(fit_intercept=False)
_ = model.fit(X, y)

y_pred = pd.Series(model.predict(X), index=y.index)
X_fore = dp.out_of_sample(steps=90)
y_fore = pd.Series(model.predict(X_fore), index=X_fore.index)

ax = y.plot(color='0.25', style='.', title="Tunnel Traffic - Seasonal Forecast")
ax = y_pred.plot(ax=ax, label="Seasonal")
ax = y_fore.plot(ax=ax, label="Seasonal Forecast", color='C3')
_ = ax.legend()

Kaggleあるあるの描画の準備がいかついというのはありますが、なんとなく予測していることが読み取れています。

実際にpendrogramを1から自作することはかなり骨が折れるため、ある程度コピペになるかなとは思っています。

どちらかというとどのようなことをしていて、どのようなアプローチを行おうとしているのか？を理解しておくといいと思います。