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

プログラムの動作

lsq-general.py は、与えられたデータセット \( (x, y) \) に対して、複数の基底関数を用いた線形最小二乗法により曲線フィッティングを行うPythonプログラムです。 主な機能と解決する課題は以下の通りです。

  • データフィッティング: Excelファイルから読み込んだ \(x, y\) 座標データに対し、指定された数の基底関数を組み合わせて最適な近似関数を決定します。

  • 基底関数の選択: プログラムに組み込まれた複数の基底関数(定数、三角関数、多項式関数、指数関数)の中から、指定した数だけフィッティングに使用できます。

  • パラメータ推定: 線形最小二乗法を用いて、各基底関数に掛かる係数を算出します。

  • 結果の可視化: 入力データ、フィッティング曲線、および計算された係数をグラフで表示し、結果を直感的に把握できるようにします。

  • 汎用性: 入力ファイル名と使用する基底関数の数をコマンドライン引数で指定できるため、様々なデータセットやフィッティング要件に対応できます。

このプログラムは、実験データや観測データに特定の数学的モデルを当てはめたい場合や、データの傾向を分析したい場合に有用です。

原理

このプログラムは、線形最小二乗法に基づいてデータフィッティングを行います。線形最小二乗法は、観測データ \( (x_i, y_i) \) (\(i=1, \dots, n\)) に対して、基底関数 \( \phi_j(x) \) の線形結合で表される近似関数 $\( f(x) = \sum_{j=0}^{m-1} c_j \phi_j(x) \)\( を構築し、観測値 \)y_i\( と近似値 \)f(x_i)\( の差の二乗和 \)S\( を最小化する係数 \)c_j\( を求める手法です。 \)\( S = \sum_{i=1}^n \left( y_i - \sum_{j=0}^{m-1} c_j \phi_j(x_i) \right)^2 \)\( ここで、\)m$ は使用する基底関数の数 (nfunc) です。

\(S\) を最小化するためには、各係数 \(c_k\) (\(k=0, \dots, m-1\)) について \( \frac{\partial S}{\partial c_k} = 0 \) となる連立方程式を解きます。これにより、以下の正規方程式が得られます。

\[ \sum_{j=0}^{m-1} \left( \sum_{i=1}^n \phi_k(x_i) \phi_j(x_i) \right) c_j = \sum_{i=1}^n y_i \phi_k(x_i) \]

この正規方程式は、行列形式で \(S_{kj} c_j = S_k\) と表すことができます。 プログラム内では、以下の行列とベクトルを構築しています。

  • ベクトル \(S_k\) (Si): $\( S_k = \sum_{i=1}^n y_i \phi_k(x_i) \)$

  • 行列 \(S_{kj}\) (Sij): $\( S_{kj} = \sum_{i=1}^n \phi_k(x_i) \phi_j(x_i) \)$

これらの SiSij を用いて、プログラムは np.linalg.inv(Sij) @ Si の計算により、係数ベクトル \(c_j\) を求めています。

プログラムで利用可能な基底関数 \( \phi_j(x) \)lsqfunc 関数で定義されており、その内容は以下の通りです。

インデックス (j)

基底関数 \( \phi_j(x) \)

flabelでの表記

0

\(1\)

1

1

\( \cos(2x) \)

cos(2x)

2

\( \sin(2x) \)

sin(2x)

3

\( \cos(x) \)

cos(x)

4

\( \sin(x) \)

sin(x)

5

\( \cos(3x) \)

cos(3x)

6

\( \sin(3x) \)

sin(3x)

7

\( x \)

x

8

\( x^2 \)

x^2

9

\( e^x \)

exp(x)

nfunc で指定された数だけ、インデックス \(0\) から順に基底関数が使用されます。例えば、nfunc=3 の場合、\(1, \cos(2x), \sin(2x)\) の3つの基底関数が使用されます。

注意点: lsqfunc 内で $e^x$ を計算するために exp(x) が呼び出されていますが、numpy から exp が明示的にインポートされていないため、Pythonの標準の math.exp が利用できない環境下(math モジュールもインポートされていないため)では NameError が発生する可能性があります。これを避けるためには、from numpy import sin, cos, tan, pi, exp とするか、np.exp(x) と記述する必要があります。

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

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

  • numpy: 数値計算(特に配列操作や線形代数演算)に使用されます。

  • pandas: データフレーム形式でのデータ操作、特にExcelファイルの読み込みに使用されます。

  • openpyxl: PandasがExcelファイルを読み書きする際に内部的に使用します。

  • matplotlib: データのプロットと可視化に使用されます。

これらのライブラリは pip コマンドを使用してインストールできます。

pip install numpy pandas openpyxl matplotlib

必要な入力ファイル

プログラムはExcel形式 (.xlsx) のファイルを読み込みます。

  • ファイル名: デフォルトでは random-sin.xlsx ですが、コマンドライン引数で任意のファイル名を指定できます。

  • ファイル形式: Microsoft Excel 2007以降の形式(.xlsx)である必要があります。

  • データ構造:

    • データはExcelファイルの最初のシートに含まれている必要があります。

    • 最初の列が独立変数 \(x\) のデータ、2番目の列が従属変数 \(y\) のデータとして扱われます。

    • 各列にはヘッダー(ラベル)が必要です。

入力ファイルの例 (random-sin.xlsx):

x

y

0.0

1.02

0.1

1.28

0.2

1.56

0.3

1.75

0.4

1.83

...

...

6.0

1.05

生成される出力ファイル

lsq-general.py は、計算結果をファイルとして保存しません。代わりに、以下の情報を標準出力(コンソール)とグラフィカルなウィンドウに出力します。

  • コンソール出力:

    • 線形最小二乗法の概要(使用される関数の数、入力ファイル名)。

    • 入力データの解析範囲(\(x\) データの最小値、最大値、計算ステップ数)。

    • 正規方程式を構成するベクトル Si および行列 Sij の内容。

    • 算出されたフィッティング関数の式 \(f(x)\) と、各基底関数に対応する係数 \(c_i\) の値。

  • グラフィカル出力:

    • matplotlib によって生成されるウィンドウに2つのプロットが表示されます。

      1. データとフィッティング曲線: 入力データ点(丸マーカー)と、計算された近似関数によるフィッティング曲線(実線)が \(x\)\(y\) の形式で表示されます。

      2. フィッティング係数: 各基底関数に対応する係数 \(c_i\) が、インデックス \(i\) に対してプロットされます。

プログラムは、グラフィカルウィンドウが表示された後、ユーザーがEnterキーを押すまで終了しません。

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

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

python lsq-general.py [入力ファイル名] [使用する基底関数の数]
  • [入力ファイル名] (オプション):

    • フィッティングに使用するExcelファイル (.xlsx) のパスを指定します。

    • 指定しない場合、デフォルトで random-sin.xlsx が使用されます。

  • [使用する基底関数の数] (オプション):

    • フィッティングに使用する基底関数の数を整数で指定します。

    • 基底関数は定義された順序(インデックス0から)で選択されます。

    • 指定しない場合、デフォルトで 4 (つまり \(1, \cos(2x), \sin(2x), \cos(x)\))が使用されます。

    • 指定できる最大値は 10 です。

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

ここでは、架空の sample_data.xlsx ファイルを準備し、それを用いてプログラムを実行する例を挙げます。

sample_data.xlsx の作成例:

sample_data.xlsx は、例えば \(y = 1.5 + 2.0 \cos(2x) - 0.5 \sin(x)\) のような関数に小さなノイズを加えたデータを含むExcelファイルであると仮定します。

x

y

0.0

1.51

0.1

1.69

0.2

1.78

0.3

1.75

0.4

1.60

...

...

実行例:

上記 sample_data.xlsx を用意し、基底関数を5つ(\(1, \cos(2x), \sin(2x), \cos(x), \sin(x)\))使用してフィッティングを行う場合。

python lsq-general.py sample_data.xlsx 5

実行結果の説明:

プログラムが実行されると、まずコンソールに以下のような情報が出力されます(具体的な数値はデータに依存します)。

Least-squares method for sum of 5 functions
nfunc=5
infile=sample_data.xlsx

Read [sample_data.xlsx]
Cal range: -0.5 - 7.5 at 0.07000000000000001 step, 101 points

Execute linear least-squares method
Vector and Matrix:
Si=
[[ 7.72886479]
 [ 0.98565171]
 [-0.1706245 ]
 [ 0.43093282]
 [-0.34757512]]
Sij=
[[10.         0.          0.         0.          0.        ]
 [ 0.         4.97543949 -0.01525381  0.03820295 -0.0345037 ]
 [ 0.        -0.01525381  5.02456051 -0.04005881  0.03434685]
 [ 0.         0.03820295 -0.04005881  4.97543949 -0.01525381]
 [ 0.        -0.0345037   0.03434685 -0.01525381  5.02456051]]

LSQ function
f(x) =    1.5 +    2.0 * cos(2x) +  -0.01 * sin(2x) +    0.0 * cos(x) +   -0.5 * sin(x)

このコンソール出力では、正規方程式の構築に用いられた Si ベクトルと Sij 行列が示され、最終的に計算された係数に基づくフィッティング関数 f(x) の式が表示されます。 上記の例では、\(c_0 = 1.5, c_1 = 2.0, c_2 = -0.01, c_3 = 0.0, c_4 = -0.5\) のような係数が得られたことがわかります。

同時に、Matplotlibによってグラフウィンドウが表示されます。

  • 左側のグラフには、入力データ点と、上記で得られた \(f(x)\) による曲線が重ねて描画されます。これにより、フィッティングの適合度が視覚的に確認できます。

  • 右側のグラフには、基底関数のインデックス \(i\) に対する係数 \(c_i\) の値がプロットされます。

グラフウィンドウが表示された後、ユーザーがコンソールで Enter キーを押すまでプログラムは待機し、Enterキーが押されると終了します。