"""
多項式最小二乗法を実行し、結果をプロットするスクリプト。
概要:
指定されたExcelファイルからデータを読み込み、多項式による最小二乗近似を行います。
近似結果の多項式係数と、入力データおよび近似曲線をグラフで表示します。
コマンドライン引数で入力ファイル名と多項式の次数を指定できます。
詳細説明:
本スクリプトは、データ `(x, y)` に基づいて `y = c_0 + c_1*x + ... + c_m*x^m` 形式の
多項式を最小二乗法でフィットさせます。
コマンドライン引数を使用することで、以下のパラメータを変更できます。
- 第1引数: 入力Excelファイルのパス (デフォルト: 'random-poly.xlsx')
- 第2引数: フィットする多項式の次数 (デフォルト: 3)
計算結果はプロットされ、入力データ、近似曲線、および計算された係数が表示されます。
関連リンク:
:doc:`lsq-polynomial_usage`
"""
import sys
import numpy as np
from numpy import sin, cos, tan, pi
from pprint import pprint
import csv
import openpyxl
import pandas as pd
import matplotlib.pyplot as plt
norder = 3
infile = 'random-poly.xlsx'
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:
norder = int(argv[2])
[ドキュメント]
def mlsq(x, y, m, iPrint = 0):
"""
与えられたデータに対して多項式最小二乗法を実行し、係数を計算する。
概要:
与えられたデータ `(x, y)` に基づき、指定された次数の多項式の係数を最小二乗法で計算します。
詳細説明:
この関数は、`y = c_0 + c_1*x + c_2*x^2 + ... + c_m*x^m` の形式で、
データ `(x, y)` に最もよくフィットする多項式の係数 `c_i` を計算します。
正規方程式 `Sij * ci = Si` を解いて係数を求めます。
オプションで中間計算のベクトルと行列を表示できます。
:param x: list or numpy.ndarray: 独立変数 (x軸) のデータポイントのリストまたはNumPy配列。
:param y: list or numpy.ndarray: 従属変数 (y軸) のデータポイントのリストまたはNumPy配列。
:param m: int: フィットする多項式の次数。
:param iPrint: int, optional: 1の場合、中間計算のベクトル (Si) と行列 (Sij) を表示する。デフォルトは0。
:returns: list: 計算された多項式の係数 `[c_0, c_1, ..., c_m]` のリスト。
"""
n = len(x)
Si = np.empty([m+1, 1])
Sij = np.empty([m+1, m+1])
for l in range(0, m+1):
v = 0.0
for i in range(0, n):
v += y[i] * pow(x[i], l)
Si[l, 0] = v
for j in range(0, m+1):
for l in range(j, m+1):
v = 0.0
for i in range(0, n):
v += pow(x[i], j+l)
Sij[j, l] = Sij[l, j] = v
if iPrint == 1:
print("Vector and Matrix:")
print("Si=")
pprint(Si)
print("Sij=")
pprint(Sij)
print("")
ci = np.linalg.inv(Sij) @ Si
ci = ci.transpose().tolist()
return ci[0]
[ドキュメント]
def main():
"""
スクリプトのメイン処理。データ読み込み、最小二乗法実行、結果表示、プロットを行う。
概要:
コマンドライン引数またはデフォルト設定に基づいて、入力データファイルと多項式の次数を決定し、
最小二乗法を実行して結果をプロットします。
詳細説明:
1. コマンドライン引数を解析し、入力ファイル名と多項式の次数を設定します。
2. pandasを使用して指定されたExcelファイルからx, yデータを読み込みます。
3. `mlsq` 関数を呼び出し、多項式の係数を計算します。
4. 計算された多項式関数と係数を標準出力に表示します。
5. 近似曲線を描画するための計算範囲とステップを決定し、データを生成します。
6. matplotlibを用いて2つのサブプロットを作成します。
- 1つ目のサブプロットには、入力データと近似曲線を表示します。
- 2つ目のサブプロットには、計算された多項式の係数を表示します。
7. プロットを表示し、ユーザーがEnterキーを押すまで待機します。
:returns: None: この関数は値を返しません。
"""
print("Least-squares method for polynomial order {}".format(norder))
print(f"norder={norder}")
print(f"infile={infile}")
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)
minx = min(x)
maxx = max(x)
if minx < 0.0:
xcal0 = minx * 1.2
else:
xcal0 = minx * 0.8
if maxx < 0.0:
xcal1 = maxx * 0.8
else:
xcal1 = maxx * 1.2
xcalstep = (xcal1 - xcal0) / (ncal - 1)
print(f"Cal range: {xcal0} - {xcal1} at {xcalstep} step, {ncal} points")
print("")
print(f"Execute linear least-squares method")
ci = mlsq(x, y, norder, iPrint = 1)
print("LSQ function")
print("f(x) = {}".format(ci[0]), end = '')
for i in range(1, norder+1):
print(" + {} * x^{}".format(ci[i], i), end = '')
print("")
xcal = [xcal0 + i * xcalstep for i in range(ncal)]
ycal = []
for i in range(ncal):
_x = xcal[i]
yl = ci[0]
for k in range(1, norder+1):
yl += ci[k] * pow(_x, k)
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(norder+1), ci, label = 'coeff', linestyle = '', 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()