"""
ガウス基底関数を用いたリッジ回帰を実行し、結果を可視化するスクリプトです。
このスクリプトは、指定されたデータセットに対してガウス基底関数を使用したリッジ回帰分析を実行します。
コマンドライン引数を通じて、入力ファイル、ガウス基底関数の数、幅、リッジ正則化パラメータを
調整することが可能です。計算された係数とフィッティング結果はグラフにプロットされ、視覚的に
分析結果を確認できます。
関連リンク: :doc:`ridge_gaussian_basis_usage`
"""
import sys
import numpy as np
from numpy import exp
import openpyxl
import pandas as pd
import matplotlib.pyplot as plt
infile = 'random-poly-Gauss.xlsx'
nG = 51
x0 = []
wG = 0.3
lmda = 0.0
xcal0 = None
xcal1 = None
ncal = 101
fontsize = 16
if __name__ == "__main__":
argv = sys.argv
narg = len(argv)
if narg >= 2:
infile = argv[1]
if narg >= 3:
nG = int(argv[2])
if narg >= 4:
wG = float(argv[3])
if narg >= 5:
lmda = float(argv[4])
[ドキュメント]
def lsqfunc(i, x):
"""
指定されたインデックスのガウス基底関数の値を計算します。
グローバル変数 `x0` (ガウス基底関数の中心) と `wG` (ガウス基底関数の幅) を使用して、
指定された `x` におけるガウス関数 `exp(-((x - x0[i]) / wG)^2)` の値を計算します。
:param i: int: ガウス基底関数の中心 `x0` のインデックス。
:param x: float: 関数値を計算する独立変数の値。
:returns: float: 計算されたガウス基底関数の値。
"""
global x0, wG
a = (x - x0[i]) / wG
return exp(-a**2)
[ドキュメント]
def Ridge(x, y, m, lmda = 0.0):
"""
ガウス基底関数を用いたリッジ回帰により、係数ベクトルを計算します。
入力データ `(x, y)` と `m` 個のガウス基底関数に対し、以下の正規方程式を解きます。
`Sij` は基底関数の内積行列、`Si` はデータと基底関数の内積ベクトルです。
正則化パラメータ `lmda` を対角要素に追加することで、過学習を抑制します。
最終的に `c = (Sij + lmda * I)^-1 * Si` を計算します。
:param x: numpy.ndarray or list: 独立変数のデータポイントの配列。
:param y: numpy.ndarray or list: 従属変数のデータポイントの配列。
:param m: int: 使用するガウス基底関数の数。
:param lmda: float, optional: リッジ正則化のパラメータ。デフォルトは0.0 (通常の最小二乗法)。
:returns: list: 計算されたガウス基底関数の係数のリスト。
"""
n = len(x)
Si = np.empty([m, 1])
Sij = np.empty([m, m])
for l in range(0, m):
v = 0.0
for i in range(0, n):
v += y[i] * lsqfunc(l, x[i])
Si[l, 0] = v
for j in range(0, m):
for l in range(j, m):
v = 0.0
for i in range(0, n):
v += lsqfunc(j, x[i]) * lsqfunc(l, x[i])
Sij[j, l] = Sij[l, j] = v
Sij[j, j] += lmda
print("Vector and Matrix:")
print("Si=")
print(Si)
print("Sij=")
print(Sij)
print("")
ci = np.linalg.inv(Sij) @ Si
ci = ci.transpose().tolist()
return ci[0]
[ドキュメント]
def main():
"""
プログラムのエントリポイントであり、リッジ回帰の実行と結果の可視化を行います。
この関数は以下のステップを実行します。
1. コマンドライン引数を解析し、入力ファイル名、ガウス基底関数の数、幅、リッジパラメータを設定します。
2. Excelファイルからデータを読み込みます。
3. ガウス基底関数の中心 `x0` をデータの範囲に基づいて生成します。また、グローバル変数 `x0` と `wG` を設定します。
4. `Ridge` 関数を呼び出してガウス基底関数の係数を計算します。
5. 計算された係数を用いて、滑らかな曲線 (フィッティング結果) を生成します。
6. 入力データとフィッティング結果、および計算された係数をグラフにプロットします。
7. ユーザーの入力があるまでプロットを表示し続けます。
:param: なし
:returns: なし
"""
global x0, wG
print("Ridge Gaussian regression")
print(f"infile={infile}")
print(f"nGauss={nG}")
print(f" width={wG}")
print(f"Ridge lambda={lmda}")
print("")
print(f"Read [{infile}]")
df = pd.read_excel(infile, engine = 'openpyxl')
labels = df.columns.to_list()
x = df[labels[0]]
y = df[labels[1]]
ndata = len(x)
xcal0 = min(x)
xcal1 = max(x)
xcalstep = (xcal1 - xcal0) / (ncal - 1)
xGmin = xcal0
xGmax = xcal1
xGstep = (xGmax - xGmin) / (nG - 1)
x0 = [xGmin + i * xGstep for i in range(nG)]
print("")
print(f"Execute linear least-squares method")
ci = Ridge(x, y, nG, lmda)
for i in range(nG):
print(f" x0[{i}]={x0[i]:6.3g}: c[{i}]={ci[i]:6.3g}")
xcal = [xcal0 + i * xcalstep for i in range(ncal)]
ycal = []
for i in range(ncal):
_x = xcal[i]
yl = 0.0
for k in range(nG):
yl += ci[k] * lsqfunc(k, _x)
ycal.append(yl)
#================================================================
# Plot
#================================================================
fig, axes = plt.subplots(1, 2, figsize = (8, 6))
axes[0].plot(x, y, label = 'input', linestyle = '', marker = 'o')
axes[0].plot(xcal, ycal, label = 'fit', linestyle = '-')
axes[0].tick_params(labelsize = fontsize)
axes[0].set_xlabel('$x$', fontsize = fontsize)
axes[0].set_ylabel('$y$', fontsize = fontsize)
axes[0].legend(fontsize = fontsize)
axes[1].plot(range(nG), ci, label = 'coeff', linestyle = '-', linewidth = 0.5, marker = 'o')
axes[1].tick_params(labelsize = fontsize)
axes[1].set_xlabel('$i$', fontsize = fontsize)
axes[1].set_ylabel('$c_i$', fontsize = fontsize)
# axes[1].legend(fontsize = fontsize)
plt.tight_layout()
plt.pause(0.1)
print("")
print("Press ENTER to terminate")
input()
if __name__ == "__main__":
main()