lsq-polynomial.py 技術ドキュメント

プログラムの動作

lsq-polynomial.py は、与えられた数値データに対して多項式による最小二乗フィッティングを実行するPythonスクリプトです。このプログラムは、Excelファイルから2列のデータを読み込み、指定された次数の多項式でデータを近似します。計算された多項式の係数はコンソールに表示され、元のデータとフィッティング曲線を比較するグラフがリアルタイムで可視化されます。これにより、実験データや観測データから特定の傾向を抽出し、その背後にある関係を多項式モデルとして表現するデータ解析タスクを支援します。

原理

このプログラムは、与えられたデータ点 \((x_k, y_k)\) （ここで \(k=0, \dots, n-1\) はデータ点のインデックス）に対し、次数の多項式 \(P(x)\) を用いた最小二乗法を適用します。多項式 \(P(x)\) は以下のように定義されます。

\[P(x) = c_0 + c_1 x + c_2 x^2 + \dots + c_m x^m = \sum_{j=0}^{m} c_j x^j\]

最小二乗法では、データ点と多項式との間の二乗誤差の合計 \(E\) を最小化するような係数 \(c_j\) を求めます。

\[E = \sum_{k=0}^{n-1} (y_k - P(x_k))^2\]

この誤差 \(E\) を各係数 \(c_i\) で偏微分してゼロと置くことで、以下の正規方程式が導出されます。

\[\sum_{j=0}^{m} \left( \sum_{k=0}^{n-1} x_k^{i+j} \right) c_j = \sum_{k=0}^{n-1} y_k x_k^i \quad \text{for } i = 0, \dots, m\]

この正規方程式は、行列形式で \(S_{ij} c_j = S_i\) と書くことができます。ここで、\(S_i\) は右辺のベクトル、\(S_{ij}\) は左辺の行列の要素を表します。

\[S_i = \sum_{k=0}^{n-1} y_k x_k^i\]

\[S_{ij} = \sum_{k=0}^{n-1} x_k^{i+j}\]

プログラム内の mlsq 関数では、これらの \(S_i\) ベクトル（Si）と \(S_{ij}\) 行列（Sij）を構築し、Numpyの線形代数ライブラリ np.linalg.inv を用いて \(S_{ij}\) の逆行列を計算し、係数ベクトル \(c_j\) （ci）を算出します。

\[c_j = (S_{ij})^{-1} S_i\]

必要な非標準ライブラリとインストール方法

このプログラムの実行には、以下のPython非標準ライブラリが必要です。

numpy: 数値計算を効率的に行うための基盤ライブラリ。
pandas: データ解析ライブラリ。Excelファイルからのデータ読み込みに利用します。
openpyxl: pandas がExcelファイルを読み書きする際に内部的に使用するライブラリ。
matplotlib: データの可視化（グラフ描画）に利用します。

これらのライブラリは、以下の pip コマンドでまとめてインストールできます。

pip install numpy pandas openpyxl matplotlib

必要な入力ファイル

プログラムはExcelファイル (.xlsx 形式) を入力として期待します。

ファイル名: デフォルトでは random-poly.xlsx ですが、コマンドライン引数で変更可能です。
ファイル形式: Microsoft Excel Workbook (.xlsx)。
データ構造: 少なくとも2列の数値データを含んでいる必要があります。プログラムは、最初の列を独立変数 \(x\)、2番目の列を従属変数 \(y\) として扱います。データにはヘッダー行を含んでいても問題ありません。

例 (random-poly.xlsx の内容):

x_value	y_value
1.0	2.1
2.0	4.3
3.0	7.8
4.0	12.0
...	...

生成される出力ファイル

lsq-polynomial.py は、ファイルを生成して保存することはありません。プログラムの実行結果は、以下の形式で表示されます。

コンソール出力:
- フィッティングに利用された多項式の次数と入力ファイル名が表示されます。
- 読み込まれたデータ範囲 (\(x_{min}, x_{max}\)) と、フィッティング曲線の計算範囲、点数が表示されます。
- mlsq 関数内で iPrint=1 が指定されているため、正規方程式を構成するベクトル Si と行列 Sij の内容が表示されます。
- 計算された多項式の係数 \(c_i\) を用いた、フィッティング関数 \(f(x)\) の形式が表示されます。
グラフィカル出力:
- matplotlib によって生成されたウィンドウに、2つのサブプロットを含むグラフが表示されます。
  - 左側のプロット: 入力データ点 (散布図) と、計算された多項式によるフィッティング曲線 (線) が重ねて表示されます。
  - 右側のプロット: 各次数 \(i\) に対応する多項式係数 \(c_i\) がプロットされます。

プログラムはグラフ表示後、ユーザーがEnterキーを押すまで一時停止します。

コマンドラインでの使用例 (Usage)

基本的な実行コマンドと引数の説明は以下の通りです。

python lsq-polynomial.py [入力ファイル名] [多項式の次数]

[入力ファイル名] (オプション):
- フィッティングに使用するExcelファイル (.xlsx) のパスを指定します。
- この引数を省略した場合、プログラムはデフォルトで random-poly.xlsx を読み込もうとします。
[多項式の次数] (オプション):
- データをフィッティングする多項式の次数を整数で指定します。
- この引数を省略した場合、プログラムはデフォルトで次数 3 を使用します。

コマンドラインでの具体的な使用例

例1: デフォルトの設定で実行

入力ファイル名 random-poly.xlsx がカレントディレクトリに存在し、次数 3 でフィッティングする場合。

python lsq-polynomial.py

実行結果の説明:

プログラムは random-poly.xlsx からデータを読み込み、3次多項式で最小二乗フィッティングを実行します。コンソールには、処理の開始メッセージ、Si と Sij 行列、そしてf(x) = c0 + c1 * x^1 + c2 * x^2 + c3 * x^3 の形式で計算された係数が表示されます。同時に、新しいウィンドウが開き、左側のグラフには元のデータと3次フィッティング曲線が表示され、右側のグラフには4つの係数 (\(c_0, c_1, c_2, c_3\)) がプロットされます。グラフが表示された後、「Press ENTER to terminate」というメッセージがコンソールに表示され、ユーザーがEnterキーを押すとプログラムが終了します。

例2: 特定のファイルと次数を指定して実行

my_experiment_data.xlsx というファイルに対して、5次多項式でフィッティングする場合。

python lsq-polynomial.py my_experiment_data.xlsx 5

実行結果の説明:

プログラムは my_experiment_data.xlsx からデータを読み込み、5次多項式で最小二乗フィッティングを実行します。コンソール出力は例1と同様ですが、norder=5 と表示され、フィッティング関数の形式が f(x) = c0 + ... + c5 * x^5 となります。グラフウィンドウでは、my_experiment_data.xlsx のデータと5次フィッティング曲線が左側のグラフに、6つの係数 (\(c_0\) から \(c_5\)) が右側のグラフにそれぞれ表示されます。グラフ表示後、ユーザーのEnterキー入力待ちとなります。