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

プログラムの動作

lsq-polynomial.py は、与えられたデータ点に対して、指定された次数の多項式を最小二乗法でフィッティングするPythonプログラムです。

このプログラムの主な機能は以下の通りです。

  1. データ生成: 0から5の範囲で101点のx座標を生成し、それに3次多項式ベースの関数とランダムなノイズを加えたy座標を生成します。

  2. 最小二乗法によるフィッティング: 生成されたデータ点に対し、コマンドライン引数で指定された次数（デフォルトは3次）の多項式を最小二乗法でフィッティングし、その係数を算出します。

  3. 結果の出力: フィッティングされた多項式の係数をコンソールに表示し、元のデータとフィッティング結果をCSVファイル lsq-polynomial.csv に保存します。

  4. グラフ表示: matplotlib を使用して、元のデータ点とフィッティング曲線が重ねて描画されたグラフを表示します。

本プログラムは、実験データや観測データから基本的な傾向を示す多項式モデルを導き出す際の基本的な手法をデモンストレーションすることを目的としています。

原理

本プログラムは最小二乗法 (Least Squares Method) を用いて多項式フィッティングを行います。 与えられた \(N\) 個のデータ点 \((x_i, y_i)\) （\(i=1, \dots, N\)）に対し、\(m\) 次の多項式 \(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\]

このとき、各データ点における真値 \(y_i\) と多項式による予測値 \(P(x_i)\) との残差の二乗和 \(S\) を最小化します。

\[S = \sum_{i=1}^N (y_i - P(x_i))^2\]

\(S\) を最小化するために、各係数 \(c_j\) で偏微分したものが0になるという条件を用います。

\[\frac{\partial S}{\partial c_k} = \frac{\partial}{\partial c_k} \sum_{i=1}^N \left(y_i - \sum_{j=0}^m c_j x_i^j\right)^2 = 0\]

この式を整理すると、以下の連立一次方程式（正規方程式）が得られます。

\[\sum_{j=0}^m \left( \sum_{i=1}^N x_i^{k+j} \right) c_j = \sum_{i=1}^N y_i x_i^k \quad \text{for } k=0, \dots, m\]

この連立一次方程式は行列形式で \(\mathbf{S_{jl}} \mathbf{c_j} = \mathbf{S_j}\) と書くことができます。 ここで、係数行列 \(\mathbf{S_{jl}}\) と右辺ベクトル \(\mathbf{S_j}\) の要素はそれぞれ以下の通りです。

\[S_{jl} = \sum_{i=1}^N x_i^{j+l}\]
\[S_j = \sum_{i=1}^N y_i x_i^j\]

本プログラムでは、mlsq 関数内でこれらの行列とベクトルを構築し、numpy.linalg.inv を用いて係数行列の逆行列を計算し、係数ベクトル \(\mathbf{c}\) を求めています。

\[\mathbf{c} = \mathbf{S_{jl}}^{-1} \mathbf{S_j}\]

算出された係数 \(c_0, c_1, \dots, c_m\) が、与えられたデータ点に最もよくフィットする \(m\) 次多項式の係数となります。

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

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

  • NumPy: 数値計算、特に多次元配列の操作や線形代数の機能を提供します。

  • Matplotlib: グラフ描画機能を提供します。

これらのライブラリは、Pythonのパッケージマネージャ pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

このプログラムは、外部からの入力ファイルを必要としません。 プログラム内部で \(x\) 座標と \(y\) 座標の擬似データを生成します。

データ生成は func(x) 関数によって行われ、以下の形式の3次多項式にランダムなオフセットが加わった値を使用します。 func(x) = r + 0.5 + 0.1 * x - 0.2 * x * x + 0.5 * x*x*x ここで r は \(0.0\) から \(50.0\) の範囲のランダムな値です。

生成される出力ファイル

プログラムを実行すると、以下のファイルが生成されます。

  • lsq-polynomial.csv:

    • CSV (Comma Separated Values) 形式のテキストファイルです。

    • ヘッダー行として x, y, y(LSQ) が記述されます。

    • 各データ行には、生成された x 座標、それに対応する元の y 座標（ランダムノイズを含む）、および最小二乗法によってフィッティングされた多項式から計算された y の値が格納されます。

さらに、プログラムの実行中に matplotlib によってグラフウィンドウが表示されます。

  • グラフウィンドウ:

    • x 軸に「x」、y 軸に「y」のラベルが付けられた散布図です。

    • 生成された元のデータ点が青色の丸 (raw data) でプロットされます。

    • 最小二乗法によってフィッティングされた多項式曲線がオレンジ色の線 (fitted) でプロットされます。

    • このウィンドウはEnterキーが押されるまで表示され続けます。

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

本プログラムは、コマンドライン引数としてフィッティングする多項式の次数を指定できます。

python lsq-polynomial.py [norder]

  • [norder]: オプション引数です。フィッティングする多項式の次数（整数）を指定します。

    • この引数を省略した場合、デフォルトで3次の多項式が使用されます。

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

例1: デフォルトの3次多項式でフィッティング

多項式の次数を指定せずにプログラムを実行します。デフォルトの3次多項式が使用されます。

python lsq-polynomial.py

実行結果の説明: コンソールには、プログラムの実行状況、生成されたxの範囲、データ点数、そして最小二乗法で計算された多項式の係数が表示されます。

Least-squares method for polynomial order 3
Make data by func(x) with random scattering
x = (0.0, 5.0, 0.05)
ndata=101
Vector and Matrix:
Si=
[[  ... ]] (計算されたS_iベクトル)
Sij=
[[  ... ]] (計算されたS_ij行列)

LSQ function
f(x) = [係数c0] + [係数c1] * x^1 + [係数c2] * x^2 + [係数c3] * x^3 
Press ENTER to exit>>

（[係数c0] などは具体的な数値が表示されます。）

lsq-polynomial.csv ファイルには以下のようなデータが書き出されます（一部抜粋）。

x,y,y(LSQ)
0.0,29.83984180860548,29.83984180860548
0.05,29.983050162598375,30.01340156711956
0.1,30.31682337754898,30.18664177119106
...

同時に、matplotlib によるグラフウィンドウが表示され、元のランダムデータ点と、それにフィットした3次曲線が描画されます。この例では、元のデータ生成関数が3次式に近い形状であるため、フィッティング曲線はデータ点の傾向をよく捉えているはずです。

例2: 1次多項式（直線）でフィッティング

多項式の次数を1としてプログラムを実行します。

python lsq-polynomial.py 1

実行結果の説明: コンソール出力の Least-squares method for polynomial order 1 の部分が変化し、表示される多項式の係数も2つ（\(c_0, c_1\)）だけになります。

Least-squares method for polynomial order 1
Make data by func(x) with random scattering
x = (0.0, 5.0, 0.05)
ndata=101
Vector and Matrix:
Si=
[[  ... ]]
Sij=
[[  ... ]]

LSQ function
f(x) = [係数c0] + [係数c1] * x^1 
Press ENTER to exit>>

lsq-polynomial.csv ファイルにも同様に、1次多項式で計算された y(LSQ) の値が書き出されます。

x,y,y(LSQ)
0.0,29.83984180860548,29.83984180860548
0.05,29.983050162598375,29.87702844199996
0.1,30.31682337754898,29.91421507539445
...

グラフウィンドウには、元のデータ点と、それを近似した直線（1次多項式）が表示されます。元のデータは3次関数に基づいているため、この直線はデータ全体の傾向を大まかに捉えるものの、例1の3次多項式ほどの精度ではフィットしないことが視覚的に確認できます。これは、次数がデータ生成モデルに対して低すぎるため、モデルの複雑さを表現しきれていないことを示しています。