区分的多項式を使ったモデルを勉強しています。そこでの発見を書いていきたいと思います。

区分的多項式とはなんぞや？

イメージ

区分的多項式によるモデルとは、複数の一次関数を局所的に用いることで予測精度を上げたモデルです。図1のように、細かく区切っていけば、一次関数を使ってもデータへの当てはまりがよくなるだろうという考え方のモデルです。

数式的理解

区分的多項式によるモデルは、加法モデルの一種です。加法モデルとは、基底関数 $h_k$ を使って次のように表現されるモデルです。

$\displaystyle{ y_i = \sum_{k=1}^K \beta_k h_k(x_i)+\varepsilon_i }$

ただし、 $i, x_i, y_i, \varepsilon_i$ はそれぞれデータの番号、説明変数、被説明変数、誤差項です。なお、簡単のため $x_i$ は一次元とします。

この加法モデルのうち、基底関数を次のように設定したモデルを区分的多項式によるモデルと呼びます。

$h_1(x) = 1 \\ h_2(x) = x \\ h_3(x) = (x-\xi_1) _ {+} \\ ... \\ h_K(x) = (x-\xi_{K-2}) _ {+} \\$

ただし、各 $\xi$ は節点（図1の直線が折れている点のx座標）、 $()_{+}$ は()の中が負であれば0、正の値であればその値を返す関数です。
この区分多項式を使うことで、

$\displaystyle{ \sum_{k=1}^K \beta_k h_k(x_i) }$

は節点で連続な関数となります（図2）。

パラメータの推定方法

ここでは最小二乗法を使ったパラメータの推定方法を述べます。先に結論を述べると

$\displaystyle{ \boldsymbol{\hat{\beta}} = (\boldsymbol{H}(\boldsymbol{x})\boldsymbol{H}(\boldsymbol{x})^\top)^{-1}\boldsymbol{H}(\boldsymbol{x})\boldsymbol{y} \tag{1} }$

で推定できます。ただし、

$\displaystyle{ \boldsymbol{\beta} = \begin{pmatrix} \beta_1 \\ \vdots \\ \beta_K \end{pmatrix}, \boldsymbol{x} = \begin{pmatrix} x_1 \\ \vdots \\ x_N \end{pmatrix}, \boldsymbol{H}(\boldsymbol{x}) = \begin{pmatrix} h_1(x_1) & \cdots & h_1(x_N) \\ \vdots & \cdots & \vdots \\ h_K(x_1) & \cdots & h_K(x_N) \end{pmatrix} }$

です。 (1)式の導出の仕方は最後におまけ的に述べようと思います。

実装

要件とコード

ここでは、実装コードを記載しようと思います。やりたいこととしては、以下の通りです。

基底関数 $h_1, h_2, \cdots, h_K$ を使った区分的多項式のモデルのパラメータ $\boldsymbol{\beta}$ を推定する
節点はデータの始点と終点間に等間隔に配置する
学習用データでパラメータを推定し、予測値とテスト用データの最小二乗誤差(Mean Square Error：MSE)で評価を行う
節点数を変えて最小のMSEとなる節点数を得る

これを実現したPythonコードが下記です（汚いコードですいません、、、）。コード中のx_train, y_trainがそれぞれ学習用データの説明変数、被説明変数です。また、x_test, y_testはそれぞれテスト用データの説明変数、被説明変数です。

# パラメータ
## MSEのリスト
MSE_list = []

## 節点の候補のリスト
node_candidate_list = range(1, 20)

# 関数
## 節点より小さい値を取る場合は0, それ以上のときはx-節点の値を取る関数を返す
def make_h(node):
    '''
    node:節点
    '''
    def h(x):
        if x<node:
            return 0
        return x-node

    return h

## 区分多項式
def h1(x):
        return 1

def h2(x):
    return x

## 区分多項式を計算して行列化したもの
def H(data_x, methods_list):
    H_tmp = np.array(range(0, len(data_x)))
    for h in methods_list:
        h_raw = np.array([h(x) for x in data_x])
        H_tmp = np.vstack((H_tmp, h_raw))

    return H_tmp[1:, :]

# 予測値を計算する関数
def predicted_y(x, beta, func_list):
    hx = np.array([h(x) for h in func_list])
    return np.dot(beta.T, hx)

# 実際にパラメータを推定しMSEで評価
for n in node_candidate_list:
    
    # 節点数
    n_node = n+2 
    
    # 節点のリスト
    node_list = np.linspace(start, end, n_node)[1:-1] # 最初と最後の要素は捨てる(最初と最後は範囲の端点だから)

    # 区分多項式のリスト
    polynomial_func_list = [h1, h2]+[make_h(node) for node in node_list]
    
    # 推定量の計算
    H_matrix = H(x_train, polynomial_func_list)
    beta_hat = np.dot(np.dot(np.linalg.inv(np.dot(H_matrix, H_matrix.T)), H_matrix), y_train)
    
    # 予測
    predicted_y_test = np.array([predicted_y(x, beta_hat, polynomial_func_list) for x in x_test])

    MSE = np.sum((predicted_y_test-y_test)**2)
    
    MSE_list.append(MSE)

結果

データの設定と実験結果を図3に示します。

考察

$y={\rm sin}(x) (0 \leq x \leq 4\pi)$ を真の関数とした場合、節点が $0, \dfrac{\pi}{2}, \pi, \dfrac{3\pi}{2}, 2\pi, \dfrac{5\pi}{2}, 3\pi, \dfrac{7\pi}{2}, 4\pi$ のときに、MSEが最小になりました。つまり、 $y={\rm sin}(x)$ の一次導関数=0になる箇所か、二次導関数=0になる箇所に節点がある場合にMSEが最小になっています。理論的になぜこのようなことになるのかはわかりませんが、感覚的には曲線の曲率が変わる箇所で節点を持つというのは、感覚的には自然なことだと思いました。

次に、 $y=x ^ 3+x ^ 2-x+1$ を真の関数とした場合についての考察です。これは、節点が13個の場合にMSEが最小になっています。しかし、図3中の節点数とMSEのグラフを見ると節点が3個以上であれば、MSEの値はあまり変わらないことがわかります。 $y=x ^ 3+x ^ 2-x+1$ も一次導関数=0もしくは二次導関数=0となる箇所が3個なので、節点数が2個から3個になったときに、MSEが大幅に小さくなったと考えています。

以上のように、理由はわかりませんが、真の関数の一次導関数と二次導関数が0になる箇所と節点には、何かしらの関係性があるのではないかと思っています。このことについて、何か理論をご存じの方は是非コメントしていただきたいです。

まとめ

本記事では、区分的多項式を基底関数にした加法モデルの説明と実装について説明しました。また、実験の結果から真の関数の一次導関数、二次導関数と節点の間に何らかの関係性があるのではないかと考察しました。

参考

(1)式の導出

ここでは、(1)式の導出を行います。行列の微分さえ理解していれば簡単です。
まず、最小二乗誤差を計算します。

$\begin{align} L(\boldsymbol{\beta}) &= ||\boldsymbol{H(\boldsymbol{x})}^\top \boldsymbol{\beta}-\boldsymbol{y}||^2 \\ &= (\boldsymbol{H(\boldsymbol{x})}^\top \boldsymbol{\beta}-\boldsymbol{y})^\top (\boldsymbol{H(\boldsymbol{x})}^\top \boldsymbol{\beta}-\boldsymbol{y}) \\ &= \boldsymbol{\beta}^\top \boldsymbol{H(\boldsymbol{x})} \boldsymbol{H(\boldsymbol{x})}^\top \boldsymbol{\beta} - \boldsymbol{\beta}^\top \boldsymbol{H(\boldsymbol{x})} \boldsymbol{y} - \boldsymbol{y}^\top \boldsymbol{H(\boldsymbol{x})}^\top \boldsymbol{\beta} + \boldsymbol{y}^\top \boldsymbol{y} \end{align}$

次に、最小二乗誤差 $L(\boldsymbol{\beta})$ を最小化する $\boldsymbol{\beta}$ を見つけるために、 $L(\boldsymbol{\beta})$ を $\boldsymbol{\beta}$ で微分します。

$\displaystyle{\begin{align} \frac{\partial L(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} &= 2\boldsymbol{H(\boldsymbol{x})} \boldsymbol{H(\boldsymbol{x})}^\top \boldsymbol{\beta} - 2\boldsymbol{H(\boldsymbol{x})} \boldsymbol{y} \end{align}}$