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

プログラムの動作

lsq-general.py は、乱数を含む合成データに対して、指定された複数の基底関数（定数項、正弦関数、余弦関数）の線形結合を用いて最小二乗法でフィッティングを行うPythonプログラムです。

主な機能:

データ生成: func(x) 関数を使用して、特定の周期関数にランダムなノイズを加えたサンプルデータを生成します。このデータは $x$ 軸の指定された範囲で等間隔に生成されます。
基底関数の定義: lsqfunc(i, x) 関数で、最小二乗法に使用する基底関数（$1, \sin(x), \cos(x), \sin(2x), \cos(2x), \sin(3x), \cos(3x)$）を定義します。使用する基底関数の数はコマンドライン引数で変更可能です。
最小二乗法によるフィッティング: mlsq(x, y, m) 関数は、生成されたデータに対して正規方程式を構築し、連立一次方程式を解くことで最適な係数を算出します。
結果の出力: フィッティングされた係数と、元のデータ、フィッティング結果のデータをCSVファイル (lsq-general.csv) に保存します。
グラフの表示: matplotlib を用いて、元のデータとフィッティング結果の曲線を同一のグラフ上にプロットし、視覚的に比較できるようにします。

解決する課題:

ノイズが含まれる観測データから、その背後にある周期的な変動パターン（フーリエ級数の一部）を近似的に抽出し、数学的なモデルで表現することを目的としています。

原理

このプログラムは、線形最小二乗法を用いてデータフィッティングを行います。与えられたデータ点 $(x_i, y_i)$ ($i=0, \dots, n-1$) に対して、モデル関数 $f(x)$ をフィットさせます。

モデル関数 $f(x)$ は、予め定義された $m$ 個の基底関数 $\phi_k(x)$ の線形結合として表現されます。

\[ f(x) = \sum_{k=0}^{m-1} c_k \phi_k(x) \]

ここで、$c_k$ は求めたい係数です。プログラムでは、lsqfunc(i, x) が基底関数 $\phi_i(x)$ に対応し、nfunc が $m$ に対応します。

基底関数: lsqfunc 関数では、以下の7つの基底関数が定義されています。

$\phi_0(x) = 1$
$\phi_1(x) = \sin(x)$
$\phi_2(x) = \cos(x)$
$\phi_3(x) = \sin(2x)$
$\phi_4(x) = \cos(2x)$
$\phi_5(x) = \sin(3x)$
$\phi_6(x) = \cos(3x)$

データ点とモデル関数との残差平方和 $S$ を最小化する係数 $c_k$ を求めます。

\[ S = \sum_{i=0}^{n-1} (y_i - f(x_i))^2 = \sum_{i=0}^{n-1} \left(y_i - \sum_{k=0}^{m-1} c_k \phi_k(x_i)\right)^2 \]

この $S$ を各係数 $c_j$ で偏微分し、その結果を0と置くことで、正規方程式と呼ばれる連立一次方程式が得られます。

\[ \frac{\partial S}{\partial c_j} = 0 \quad \Rightarrow \quad \sum_{k=0}^{m-1} \left( \sum_{i=0}^{n-1} \phi_j(x_i) \phi_k(x_i) \right) c_k = \sum_{i=0}^{n-1} y_i \phi_j(x_i) \]

これを行列形式で記述すると、$S_{jk} c_k = S_j$ となります。ここで、 $$ S_j = \sum_{i=0}^{n-1} y_i \phi_j(x_i) $$ S_{jk} = \sum_{i=0}^{n-1} \phi_j(x_i) \phi_k(x_i) $$

プログラム内の mlsq 関数では、Si がベクトル $S_j$ に、Sij が行列 $S_{jk}$ に対応します。 Sij 行列は定義上対称行列となるため、プログラムでは Sij[j, l] = Sij[l, j] = v のように計算の効率化を図っています。この連立一次方程式は、係数ベクトル $c = (c_0, c_1, \dots, c_{m-1})^T$ について以下のように解かれます。

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

プログラムでは、numpy.linalg.inv(Sij) @ Si を用いて、行列 $S_{jk}$ の逆行列を計算し、それをベクトル $S_j$ にかけることで係数 $c_k$ を求めています。

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

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

NumPy: 高度な数値計算（配列操作、線形代数演算など）のために使用されます。
Matplotlib: グラフの描画のために使用されます。

これらのライブラリは、Pythonのパッケージマネージャーである pip を使ってインストールできます。

pip install numpy matplotlib

必要な入力ファイル

lsq-general.py は、外部からの入力ファイルを必要としません。データはプログラム内部で定義された func(x) 関数によって生成されます。

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

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

lsq-general.csv:
- 内容: 元のデータ点 (x, y) と、最小二乗法でフィッティングされた曲線上の点 (y(LSQ)) がカンマ区切り形式で記述されます。
- フォーマット: 最初の行はヘッダ (x,y,y(LSQ)) で、それ以降の行は各データ点に対応する数値データが格納されます。

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

グラフウィンドウ:
- 内容: 横軸を x、縦軸を y とし、生成された元のデータ点がプロット（マーカー o）され、それに加えて最小二乗法でフィットされた曲線がプロット（線）されます。凡例によってどちらのデータが何を示しているかが明記されます。

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

lsq-general.py は、コマンドライン引数でフィッティングに使用する基底関数の数を指定できます。

python lsq-general.py [nfunc]

[nfunc]: オプション引数です。最小二乗法に使用する基底関数の数を整数で指定します。
- 指定しない場合、デフォルト値の 7 が使用されます。
- 有効な値は 1 から 7 までの整数です。
  - nfunc=1: 定数項のみ
  - nfunc=3: 定数項、$\sin(x)$, $\cos(x)$
  - nfunc=7: 全ての基底関数（定数項, $\sin(x)$, $\cos(x)$, $\sin(2x)$, $\cos(2x)$, $\sin(3x)$, $\cos(3x)$）

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

1. デフォルトの基底関数数 (7) で実行

nfunc を指定せずに実行すると、デフォルトの 7 個の基底関数（定数項、$\sin(x)$、$\cos(x)$、$\sin(2x)$、$\cos(2x)$、$\sin(3x)$、$\cos(3x)$）が使用されます。

python lsq-general.py

実行結果:

端末出力:

Least-squares method for sum of 7 functions
Make data by func(x) with random scattering
x = (0.0, 5.0, 0.05)
ndata=101

Vector and Matrix:
Si=
[[...]] (7x1ベクトル)
Sij=
[[...]] (7x7行列)

LSQ function
f(x) = 係数c0 + 係数c1 * sin(x) + ... + 係数c6 * cos(3x)
Press ENTER to exit>>

（Si と Sij の内容は実際の数値によって異なります）

ファイル出力: lsq-general.csv が生成され、x, y, y(LSQ) のデータが記録されます。
グラフ表示: デフォルトの7つの基底関数でフィッティングされた曲線と元のデータが表示されます。元のデータは 0.5 + 0.1 * sin(2.0*x) + 0.3 * cos(2.0*x) にノイズが加わっているため、sin(2x) と cos(2x) の項の係数が大きく、データによくフィットした曲線が描画されます。

2. 基底関数数 3 で実行

nfunc に 3 を指定して実行すると、3個の基底関数（定数項、$\sin(x)$、$\cos(x)$）が使用されます。

python lsq-general.py 3

実行結果:

端末出力:

Least-squares method for sum of 3 functions
Make data by func(x) with random scattering
x = (0.0, 5.0, 0.05)
ndata=101

Vector and Matrix:
Si=
[[...]] (3x1ベクトル)
Sij=
[[...]] (3x3行列)

LSQ function
f(x) = 係数c0 + 係数c1 * sin(x) + 係数c2 * cos(x)
Press ENTER to exit>>

（Si と Sij の内容は実際の数値によって異なります）

ファイル出力: lsq-general.csv が生成され、x, y, y(LSQ) のデータが記録されます。
グラフ表示: 3つの基底関数でフィッティングされた曲線と元のデータが表示されます。元のデータには sin(2x) や cos(2x) の成分が含まれていますが、フィッティングモデルにはそれらの項がないため、デフォルト（7つの基底関数）の場合に比べてフィットの精度が低い曲線が描画されることが予想されます。