ほぼPython

Not技術ブログBut勉強ブログ 内容には誤りがあることが多いです

ガウス関数を基底関数にした線形基底関数モデル

これまでは、単なる線形回帰をしてきましたが、今回は、ガウス関数を基底関数として線形回帰を行いました。
特にデータとかもなかったので、今回は、正弦関数の値をそのまま使ってみました。

import numpy as np
import matplotlib.pyplot as plt

#データセット
X = np.linspace(0,np.pi,10) #0からπまで10等分
T = np.sin(X) #そのsinの値
X_n =X.shape[0] #xのサイズ

まず、ガウス関数の定義をしました。

#ガウス関数
def gauss(x,m,s):
    X = x - m
    gauss_calc = np.exp((-1*(X*X)/(2*s*s)))
    return gauss_calc


例えば、

print(gauss(np.array([1,2,4,5]),3,5))
#=> [0.92311635 0.98019867 0.98019867 0.92311635]

このように、ndarray型を投げると、ndarray型で返ってきます。ちなみに、第二引数はガウス関数の中心位置で、第三引数はガウス関数の広がり具合です。これらは設計者が決めるパラメータらしいです。数学に疎いので最適なパラメータの決め方は見当もつきません。

次に、それぞれの基底関数の最適な係数を決定します。最小二乗法を使います。各偏微分 = 0 から導き出される  w = (X^{T}X)^{-1}X^{T}t を使います。便利ですね。今回は基底関数を用いているので、  w = (Φ^{T}Φ)^{-1}Φ^{T}t となりますね。

ただし、Φは以下のようにあらわされます。

\[
{ \scriptsize
Φ =
\left(
\begin{array}{lllll}
Φ_0(x_0) &Φ_1(x_0) &\cdots &Φ_{ \scriptsize M}(x_0)\\
Φ_0(x_1) & Φ_1(x_1)&\cdots &Φ_ { \scriptsize M}(x_1) \\
\vdots & \vdots &\vdots & \vdots \\
Φ_0{ \scriptsize (x_{N-1}) }& Φ_1 { \scriptsize (x_{N-1}) }&\cdots &Φ_ { \scriptsize M}{ \scriptsize (x_{N-1}) }\\
\end{array}
\right)}
\]

ちなみに、一般化すると、各x_kは、D次元となっています。つまり、Φ中の  x_kはそれぞれ、D次元のベクトルです。
ただし、今回はXは一次元なので一次元です。

コードはこのようになります。

#最小二乗法
def fit_gauss():
    M = 2 #ガウス関数の数
    m = np.linspace(0,2,M) #ガウス関数の中心
    s = m[1] - m[0] #ガウス関数の広がり
    gauss_arr = np.ones((X_n,M+1)) #行がデータの数、列がガウス関数の数+1の1埋め行列を作る

    for i in range(len(m)):
        gauss_arr[:,i]=gauss(X,m[i],s) #ポイント

    a = np.linalg.inv(np.dot(gauss_arr.T,gauss_arr))
    b = np.dot(gauss_arr.T,T)
    w = np.dot(a,b) #最小二乗法によるwの計算

    recursion(X,w,s,m) #後述
    return m,s,w

#ポイント の部分が大事なところです。gauss_arrという行列に、各ガウス関数の計算結果を入れていきます。

まず、i=0について実行されて、以下のようになります。

[[1.         1.         1.        ]
 [0.98488453 1.         1.        ]
 [0.94089523 1.         1.        ]
 [0.87190236 1.         1.        ]
 [0.78372747 1.         1.        ]
 [0.68333383 1.         1.        ]
 [0.5779249  1.         1.        ]
 [0.47411154 1.         1.        ]
 [0.37727695 1.         1.        ]
 [0.29121293 1.         1.        ]]

これはどういう計算がされているかというと、Φ_0に、各Xの値が入れられ、縦に並べられています。

次に、i=1について実行され、

[[1.         0.60653066 1.        ]
 [0.98488453 0.71127372 1.        ]
 [0.94089523 0.80907988 1.        ]
 [0.87190236 0.89272289 1.        ]
 [0.78372747 0.95546014 1.        ]
 [0.68333383 0.99192563 1.        ]
 [0.5779249  0.99888682 1.        ]
 [0.47411154 0.97571748 1.        ]
 [0.37727695 0.92449065 1.        ]
 [0.29121293 0.84967256 1.        ]]

と、Φ_1に、各Xの値が入れられ、縦に並べられています。

ちなみに最後の「1」は切片のために必要な項です。また、今回は、Xが一次元なので各要素は数値になっていますが、、D次元の場合は各要素がベクトルで表されるはずです。(実行してないので自信はないです。)

これらを使って、wがうまい具合に計算されます。

最後に、得られたデータをもとにグラフを描いてみます。グラフを描くために recursion関数を作ります。

#回帰式で予測値を計算
def recursion(x,m,s,w):
    w_len = len(w) - 1
    y = np.zeros_like(x)
    for i in range(w_len):
        y = y + w[i]*gauss(x,m[i],s)
    y = y + w[w_len]
    return y

ここの計算方法は若干技巧的な気がします。

まず、0埋めの行列を引数x(つまりX)と同じサイズで作ります。今回、Xは一次元なので、1行10列の行列(配列)ですね。
また、今回、wは3つ求められていますが、w0, w1 ,w2のうち、最後のw2は切片として使うので、ループの範囲に入らないように w_len は wのサイズから1引いたサイズとします。それで、 yに計算結果を加算していきます。

どういうことかというと、 y = w0 * Φ_0(x) + w1 * Φ_1(x) + w2 という回帰式が求まったはずですが、ループの中では、以下のように回していくことで、結果的に各Xについて  y = w0 * Φ_0(x) + w1 * Φ_1(x) を計算します。

まず、各Xに対して w0 * Φ_0(x)を計算します。

[1.24619314 1.22735635 1.17253718 1.08655874 0.9766758  0.85156593
 0.72020604 0.59083454 0.47015995 0.36290756]  #1回目

次のループではこの計算結果に w1 * Φ_1(x) を加算します。これにより、各Xについて y = w0 * Φ_0(x) + w1 * Φ_1(x) が計算できたことになりますね。

[3.40887335 3.76351436 4.05743839 4.26970222 4.38351878 4.38843233
 4.28189366 4.06990826 3.76657637 3.39254836]

そして、最後に y = y + w[w_len] とし、各yから切片分引くことで、

 [-0.00223908  0.35240193  0.64632596  0.85858979  0.97240636  0.97731991
  0.87078124  0.65879584  0.35546395 -0.01856407]

という結果を得られました。

あとは、

#回帰式描画
plt.figure(figsize=(5,5))
M,S,W = fit_gauss()
Y = recursion(X,M,S,W)
W_format = ['{0:.2f}'.format(s) for s in W]
plt.plot(X,Y)
plt.plot(X,T,marker='o',linestyle='None')
plt.show()

とすれば、グラフが書けます。
f:id:short_4010:20180307210614p:plain
結構カクカクしてて正弦関数かと言われると微妙ですが、まぁ、いい感じですね。

手計算だったら絶対に無理なような計算でもpythonを使えば、比較的簡単に実装できました。
回帰の勉強はこれくらいにして次回からは教師あり学習の「分類」の分野に入っていきたいと思います。