はじめに

本稿では，ガウス動径基底関数（Gaussian Radial Basis Function, Gaussian RBF）を基底として用いた回帰分析の理論的背景と，Pythonによる実装について述べる．
本手法は，多項式回帰とは異なり，局所的な基底関数の重ね合わせによってデータのパターンを表現するものであり，非線形かつ複雑な関数形状への適合において高い柔軟性を示す．

なお，本稿は以下の記事を前提知識として参照されたい．

1. 理論的背景

1.1 ガウス動径基底関数とは

ガウス動径基底関数は，次式で定義される．

$$
\phi_p(x) = a_p \exp\left(-\frac{(x - x_p)^2}{2\sigma_p^2}\right)
$$

ここで，$x_p$ は各基底関数の中心（データ点），$a_p$ は振幅パラメータ，$\sigma_p$ は幅（標準偏差）を示す．
各基底関数はガウス分布の形状を持ち，中心 $x_p$ から離れるに従い値が指数的に減衰する局所的な性質を有する．

1.2 回帰モデルの定義

$N$ 個のデータ点 ${x_1, x_2, \ldots, x_N}$ に対して，各データ点を中心とした $N$ 個のガウス基底関数を配置する．
予測値 $\hat{y}(x)$ は，これらの線形結合として以下のように表される．

$$
\hat{y}(x) = \sum_{p=1}^{N} a_p \exp\left(-\frac{(x - x_p)^2}{2\sigma_p^2}\right)
$$

パラメータ ${a_p, \sigma_p}_{p=1}^{N}$ を学習データに対して最適化することが，本手法の目的である．

1.3 損失関数

モデルの予測誤差を評価する損失関数 $E$ を，平均二乗誤差（MSE）として定義する．

$$
E(\mathbf{a}, \boldsymbol{\sigma}) = \frac{1}{N} \sum_{i=1}^{N} \left(y_i - \hat{y}(x_i)\right)^2
= \frac{1}{N} \sum_{i=1}^{N} \left(y_i - \sum_{p=1}^{N} a_p \exp\left(-\frac{(x_i - x_p)^2}{2\sigma_p^2}\right)\right)^2
$$

2. 最急降下法によるパラメータ最適化

損失関数 $E$ を最小化するために，最急降下法（Gradient Descent）を適用する．
各パラメータの更新則は次式で与えられる．

$$
a_p \leftarrow a_p - \alpha \frac{\partial E}{\partial a_p}
$$

$$
\sigma_p \leftarrow \sigma_p - \alpha \frac{\partial E}{\partial \sigma_p}
$$

ここで $\alpha$ は学習率である．

偏微分の計算には数値微分を用いる．微小量 $h$ を導入し，

$$
\frac{\partial E}{\partial a_p} \approx \frac{E(\ldots, a_p + h, \ldots) - E(\ldots, a_p, \ldots)}{h}
$$

$$
\frac{\partial E}{\partial \sigma_p} \approx \frac{E(\ldots, \sigma_p + h, \ldots) - E(\ldots, \sigma_p, \ldots)}{h}
$$

として近似する．これにより，解析的な偏微分の導出を必要とせず，汎用的な勾配計算が実現される．

3. Pythonによる実装

3.1 コード全文

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

L = 2
num = 10
x = np.linspace(-L, L, num)
y = np.sin(x) + 0.1 * np.random.normal(0, 1, num)

# 数値微分に用いる微小変化分
h = 1.0 * 10**-5

A_ary = np.ones(num)
Sigma_ary = np.ones(num)


def gauss(a, sigma, x0, x):
    return a * np.exp(-((x - x0)**2) / (2 * sigma**2))


def error(A_ary, Sigma_ary):
    err = 0
    for i in range(num):
        X = 0
        for p in range(num):
            a = A_ary[p]
            sigma = Sigma_ary[p]
            X = X + gauss(a, sigma, x[p], x[i])
        err = err + (1 / num) * (y[i] - X)**2
    return err


# 積算回数
iterations = 1000
# 学習率
alpha = 0.01
# 損失関数の履歴
error_pre = []

for i in range(iterations):
    for p in range(0, num):
        A_ary2 = np.copy(A_ary)
        Sigma_ary2 = np.copy(Sigma_ary)
        A_ary2[p] = A_ary[p] + h
        Sigma_ary2[p] = Sigma_ary[p] + h
        d_err_A = (error(A_ary2, Sigma_ary) - error(A_ary, Sigma_ary)) / h
        d_err_Sigma = (error(A_ary, Sigma_ary2) - error(A_ary, Sigma_ary)) / h
        A_ary[p] = A_ary[p] - alpha * d_err_A
        Sigma_ary[p] = Sigma_ary[p] - alpha * d_err_Sigma

# 予測データ描写用
Num = 100
xx_ary = np.linspace(-L, L, Num)
yy_ary = np.zeros(Num)
for k in range(Num):
    yyy = 0
    for p in range(num):
        a = A_ary[p]
        sigma = Sigma_ary[p]
        x0 = x[p]
        yyy += gauss(a, sigma, x0, xx_ary[k])
    yy_ary[k] = yyy

title = "ガウス動径基底関数による回帰分析"
plt.title(title)
plt.scatter(x, y, label="データ", color="blue")
plt.plot(xx_ary, yy_ary, label="単回帰予測", color="red")
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.savefig(title + ".png")
plt.show()

3.2 実装上のポイント

データ生成

x = np.linspace(-L, L, num)
y = np.sin(x) + 0.1 * np.random.normal(0, 1, num)

区間 $[-2, 2]$ に $N=10$ 点を等間隔に配置し，$\sin(x)$ に標準偏差 $0.1$ のガウスノイズを加算したものを観測データとして用いる．

ガウス基底関数の定義

def gauss(a, sigma, x0, x):
    return a * np.exp(-((x - x0)**2) / (2 * sigma**2))

振幅 a，幅 sigma，中心 x0 をパラメータとして受け取り，スカラーまたは配列 x に対してガウス関数値を返す．

損失関数の計算

def error(A_ary, Sigma_ary):
    err = 0
    for i in range(num):
        X = 0
        for p in range(num):
            X = X + gauss(A_ary[p], Sigma_ary[p], x[p], x[i])
        err = err + (1 / num) * (y[i] - X)**2
    return err

全データ点 $x_i$ における予測値と観測値の平均二乗誤差を計算する．

勾配の数値計算とパラメータ更新

for p in range(0, num):
    A_ary2 = np.copy(A_ary)
    Sigma_ary2 = np.copy(Sigma_ary)
    A_ary2[p] = A_ary[p] + h
    Sigma_ary2[p] = Sigma_ary[p] + h
    d_err_A     = (error(A_ary2,  Sigma_ary ) - error(A_ary, Sigma_ary)) / h
    d_err_Sigma = (error(A_ary,  Sigma_ary2) - error(A_ary, Sigma_ary)) / h
    A_ary[p]     = A_ary[p]     - alpha * d_err_A
    Sigma_ary[p] = Sigma_ary[p] - alpha * d_err_Sigma

$p$ 番目のパラメータのみ $h$ だけ変化させたコピーを作成し，前進差分によって偏微分を数値的に近似する．
np.copy() を用いることで，他のパラメータへの影響を排除していることに注意されたい．

4. 実験条件とハイパーパラメータ

パラメータ	値
データ点数 $N$	10
入力区間 $[-L, L]$	$[-2, 2]$
ノイズ標準偏差	0.1
反復回数	1000
学習率 $\alpha$	0.01
数値微分刻み幅 $h$	$10^{-5}$
初期パラメータ $(a_p, \sigma_p)$	すべて 1.0

5. 実験結果

以下のグラフは，1000回の反復最適化後に得られた回帰曲線を示す．
青点が観測データ（ノイズ付き $\sin$ 関数），赤線が学習済みモデルによる予測曲線である．

グラフより，$\sin(x)$ の形状を概ね捉えた曲線が得られていることが確認できる．
各データ点を中心とするガウス基底関数が局所的に寄与し合うことで，多項式回帰では表現困難な非線形パターンへの適合が実現されている．

6. 考察

6.1 本手法の利点

局所的な表現力：各基底関数が特定領域に限定された影響を持つため，局所的なパターンの変化に対応しやすい．
数値微分による汎用実装：解析的な偏微分を導出することなく，任意の損失関数に対して勾配計算が可能である．
直感的なパラメータ解釈：$a_p$ は各基底の寄与度，$\sigma_p$ は影響範囲を直接制御する．

6.2 本手法の課題

計算コスト：$N$ 個の基底関数に対してそれぞれ独立した数値微分を行うため，計算量は $O(N^2)$ のオーダーとなる．データ点数の増大に対して計算時間が急増する点に注意が必要である．
局所最適解への収束：最急降下法は凸でない損失関数曲面において局所最適解に陥る可能性があり，初期値の設定が収束結果に影響を与える．
$\sigma_p$ の符号制約：本実装では $\sigma_p$ の非負制約を明示的に設けていない．学習の過程で $\sigma_p$ が負値をとりうる点は，実用上の注意事項となる．

7. おわりに

本稿では，ガウス動径基底関数を用いた回帰分析の理論的背景を整理し，最急降下法と数値微分を組み合わせた実装を示した．
本手法は，基底関数の種類やパラメータ最適化のアルゴリズムを変更することで，様々な非線形回帰問題へ容易に拡張可能である．

今後の発展として，以下のような展開が考えられる．

正則化項の導入による過学習の抑制
Adam などの適応的最適化アルゴリズムへの置き換え
カーネル回帰やガウス過程回帰との比較

参考文献・関連記事