線形基底関数モデル③　ガウス基底モデル

ガウス基底モデル

　前回までは多項式近似を扱いましたが、今回は
　
\[\phi(x)=\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}\]
で定義されるガウス関数を基底関数として用いた ガウス基底モデル を扱います。

ガウス基底モデル①

　$\sigma=10,\ \mu=20k\ (k=0,\ 1,\ 2,\ 3,\ 4)$ とした 5 個の ガウス関数
　
\[\phi_k(x)=\exp\left\{-\frac{(x-20k)^2}{2\cdot 10^2}\right\}\quad (k=0,\ 1,\ 2,\ 3,\ 4)\]
および、定数関数 $\phi_5(x)=1$ の線形結合
　
\[\sum_{k=0}^{5}a_k\phi_k(x)\]
を近似関数として与えてパラメータ $a_k$ を最適化してみます。前回記事 に掲載したリストM3-B-1 を実装したうえで、以下のコードを加えてください（当サイトにおいては、たとえばリスト X0-A-3 を実行するためには、X0-A-1, X0-A-2 の実装が前提になります）。

# https://python.atelierkobato.com/gauss-model/

# 線型基底関数モデル リストM3-B-3

# 入力データベクトル
x = np.array([35, 16, 22, 43, 5,
              66, 20, 13, 52, 1,
              39, 62, 45, 33, 8,
              28, 71, 24, 18, 3])

# 目標データベクトル
y = np.array([85.19, 58.93, 64.27, 68.91, 21.27,
              68.88, 60.07, 55.18, 88.08, 8.89,
              82.31, 81.18, 80.76, 78.98, 37.55,
              75.9, 72.39, 69.51, 62.04, 12.47])

# ガウス関数を定義
def gauss(x, mu = 0, sigma = 1):
    return np.exp(-(x - mu)**2 / (2 * sigma**2))

# 基底関数を定義
def phi_0(x):
    return gauss(x, 0, 10)

def phi_1(x):
    return gauss(x, 20, 10)

def phi_2(x):
    return gauss(x, 40, 10)

def phi_3(x):
    return gauss(x, 60, 10)

def phi_4(x):
    return gauss(x, 80, 10)

def phi_5(x):
    return 1

# 基底関数のリスト
bf = [phi_0, phi_1, phi_2, phi_3, phi_4, phi_5]

# Fit_funcオブジェクトを作成
z = Fit_func(x, y, bf)

# params()メソッドで最適化されたパラメータを取得
prm = np.round(z.params(), 3)
print("p = {}".format(prm))

# sd()メソッドで標準偏差を取得
sd = np.round(z.sd(), 3)

# Figureを作成
fig = plt.figure(figsize = (8, 6))

# FigureにAxesを１つ追加
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Gaussian base".format(sd), fontsize = 16)

# テキストボックスの書式辞書を作成
boxdic = {"facecolor" : "white", "edgecolor" : "gray",}

# テキストボックスを表示
ax.text(3, 90, "SD = {}".format(sd),
        size = 14, linespacing = 1.5, bbox = boxdic)

# 目盛線を表示
ax.grid()

# 軸範囲を設定
ax.set_xlim([0, 70])
ax.set_ylim([0, 100])

# 軸ラベルを設定
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)

# 身長と体重データの散布図
ax.scatter(x, y, color = "blue")

# line()メソッドを用いて近似曲線を描画
x2 = np.linspace(0, 70, 100)
y2 = z.line(x2)
ax.plot(x2, y2, color = "red")

plt.show()

p = [-101.892  -38.274  -34.094  -30.118  -54.606  122.477]

　標準偏差は $5.004$ です。これは以前に扱った 2次関数モデルより少し良い近似ですが、2次関数と対数関数の混合基底モデルには及びません。ただし、ガウス関数は $\mu$ と $\sigma$ というパラメータを含んでいるので、これを少しずつ調整することで改良できそうです。
　

独学プログラマー Python言語の基本から仕事のやり方まで 中古価格 ￥1,537から (2019/7/30 20:25時点)

ガウス基底モデル②

　そこで、基底関数 $\phi_k(x)\, (k=0,\ 1,\ 2,\ 3,\ 4)$ の $\mu$ に $-10$ ～ $10$ の乱数を、$\sigma$ に $-5$ ～ $5$ の乱数を加えて基底関数とします。
　
\[\begin{align*}
&\phi_0(x)=\exp\left\{-\frac{x-r_{\mu}^2}{2(10+ r_{\sigma})^2}\right\}\\[6pt]
&\phi_1(x)=\exp\left\{-\frac{x-(20+ r_{\mu})^2}{2(10+ r_{\sigma})^2}\right\}\\[6pt]
&\phi_2(x)=\exp\left\{-\frac{x-(40+ r_{\mu})^2}{2(10+ r_{\sigma})^2}\right\}\\[6pt]
&\phi_3(x)=\exp\left\{-\frac{x-(60+ r_{\mu})^2}{2(10+ r_{\sigma})^2}\right\}\\[6pt]
&\phi_4(x)=\exp\left\{-\frac{x-(80+ r_{\mu})^2}{2(10+ r_{\sigma})^2}\right\}\\[6pt]
&\phi_5(x)=1\end{align*}\]
　$r_{\mu}$ と $r_{\sigma}$ は、それぞれ相異なる $-10$ ～ $10$, $-5$ ～ $5$ の一様乱数です。これらの基底関数の線形結合で回帰を実行して標準偏差 SD を計算させます。次に乱数を入れ替えて同じ計算を実行します。これを10 回繰り返して、その中から標準偏差 SD が最も低い関数を近似関数として採用することにします。つまり、個々のガウス関数の中心位置と横幅を無作為にあれこれ変えながら重ね合わせて、その中から最適な組み合わせを見つけようとするアルゴリズムです。

# https://python.atelierkobato.com/gauss-model/

# 線型基底関数モデル リストM3-B-4

# 入力データベクトル
x = np.array([35, 16, 22, 43, 5,
              66, 20, 13, 52, 1,
              39, 62, 45, 33, 8,
              28, 71, 24, 18, 3])

# 目標データベクトル
y = np.array([85.19, 58.93, 64.27, 68.91, 21.27,
              68.88, 60.07, 55.18, 88.08, 8.89,
              82.31, 81.18, 80.76, 78.98, 37.55,
              75.9, 72.39, 69.51, 62.04, 12.47])

# ガウス関数を定義
def gauss(x, mu = 0, sigma = 1):
    return np.exp(-(x - mu)**2 / (2 * sigma**2))

# 標準偏差sdの初期値
sd = 10.0

# 乱数シードを固定
np.random.seed(10)

# 試行回数(Number of trials)を設定
nt = 10

# 乱数を使った基底関数の調整
for k in range(nt):
    rdmu = 20 * np.random.rand(5) - 10
    rdsg = 10 * np.random.rand(5) - 5
    
    def phi_0(x):
        return gauss(x, rdmu[0], 10 + rdsg[0])
    
    def phi_1(x):
        return gauss(x, 20 + rdmu[1], 10 + rdsg[1])
    
    def phi_2(x):
        return gauss(x, 40 + rdmu[2], 10 + rdsg[2])

    def phi_3(x):
        return gauss(x, 60 + rdmu[3], 10 + rdsg[3])
    
    def phi_4(x):
        return gauss(x, 80 + rdmu[4], 10 + rdsg[4])
    
    def phi_5(x):
        return 1

    # 基底関数のリスト
    func_list = [phi_0, phi_1, phi_2, phi_3, phi_4, phi_5]
    
    # Fit_funcオブジェクトを作成
    w = Fit_func(x, y, func_list)

    # 標準偏差と基底関数リストの更新
    if w.sd() < sd:
        sd = w.sd()
        bf_list = func_list

# 標準偏差の値を丸める
sd = np.round(sd, 3)
    
# 採用された基底関数でFit_funcオブジェクトを作成
z = Fit_func(x, y, bf_list)

# Figureを作成
fig = plt.figure(figsize = (8, 6))

# FigureにAxesを１つ追加
ax = fig.add_subplot(111)

# Axesのタイトルを設定
ax.set_title("Random Gaussian base", fontsize = 16)

# テキストボックス書式辞書を作成
boxdic = {"facecolor" : "white", "edgecolor" : "gray",}

# テキストボックスを表示
ax.text(3, 85, "NT = {}\nSD = {}".format(nt,sd),
        size = 14, linespacing = 1.5, bbox = boxdic)

# 目盛線を表示
ax.grid()

# 軸範囲を設定
ax.set_xlim([0, 70])
ax.set_ylim([0, 100])

# 軸ラベルを設定
ax.set_xlabel("Age", fontsize = 14)
ax.set_ylabel("Weight [kg]", fontsize = 14)

# 身長と体重データの散布図
ax.scatter(x, y, color = "blue")

# line()メソッドを用いて近似曲線を描画
x2 = np.linspace(0, 70, 100)
y2 = z.line(x2)
ax.plot(x2, y2, color = "red")

plt.show()

　$\mu$ と $\sigma$ を調整した結果、標準偏差は 4.359 まで下がりました。グラフの左上にある NT は Number of trials (試行回数) を表しています。もちろん試行回数を増やせば、それだけたくさんの組合わせを検討するので、標準偏差はさらに下がっていきますが、それが必ずしも良いモデルになるとは限りません。それについては次回記事で扱います。