技術ドキュメント: peakfit-scipy-minimize.py

プログラムの動作

peakfit-scipy-minimize.py は、与えられたデータ点に対してガウス関数モデルをフィッティングするPythonプログラムです。scipy.optimize.minimize ライブラリの強力な最適化機能を利用し、最小二乗法に基づいてモデルパラメータを最適化します。

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

ガウス関数フィッティング: 複数のガウスピークの線形結合で構成されるモデル関数を使用して、与えられたデータに最適なピークパラメータ（強度 $I_0$、中心 $x_0$、半値幅 $w$）を決定します。
非線形最適化: Scipy の minimize 関数（デフォルトでは共役勾配法 (CG) を使用）を用いて、観測データとモデル予測の二乗誤差の和を最小化します。
解析的勾配の利用: 目的関数の解析的な勾配（ヤコビアン）を提供することで、数値微分に比べて最適化の収束速度と精度を向上させます。
データ生成とシミュレーション: フィッティングの対象となる「観測データ」は、スクリプト内で定義された既知のガウス関数パラメータからシミュレートされます。これにより、フィッティングアルゴリズムの性能を評価できます。
結果の可視化: Matplotlib を使用して、元のデータ、初期推定値によるモデル、および最適化されたモデルをリアルタイムでグラフ表示します。また、最適化の収束履歴も表示されます。
CSV出力: フィッティングの初期データ、最終結果、および最適化の収束履歴をCSVファイルとして保存します。

本プログラムは、実験データやシミュレーション結果に含まれるピーク構造を定量的に解析し、その特性を抽出する際に役立ちます。

原理

モデル関数

本プログラムでは、以下のガウス関数（または複数のガウス関数の線形結合）をデータにフィッティングします。

\[ y_c(x) = \sum_{p} I_{0,p} \exp\left(-\left(\frac{C}{w_p}(x - x_{0,p})\right)^2\right) \]

ここで、各ピーク $p$ について、

$I_{0,p}$：ピークの高さ（強度）
$x_{0,p}$：ピークの中心位置
$w_p$：ピークの半値幅の半分（Half Width at Half Maximum, HWHM）
$C = 0.832554611 \approx \sqrt{\ln 2}$

このモデルにおいて、$x=x_{0,p} \pm w_p$ のときにピーク高さが半減することから、$w_p$ が半値幅の半分であることがわかります。

目的関数（最小二乗法）

フィッティングの目的は、モデル関数 $y_c(x)$ が観測データ $y_d(x)$ に最もよく一致するようなパラメータ $P = (I_0, x_0, w)$ を見つけることです。この一致度を評価するために、観測データとモデル予測の二乗誤差の和 $S^2$ を最小化します。

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

ここで、$N$ はデータ点の総数、$x_i$ と $y_d(x_i)$ はそれぞれ $i$ 番目のデータ点のX座標とY座標を表します。

最適化アルゴリズム

目的関数 $S^2(P)$ を最小化するために、scipy.optimize.minimize 関数が使用されます。この関数は、非線形最適化問題に対する様々なアルゴリズムを実装しており、本プログラムではデフォルトで method="cg" (共役勾配法) が選択されています。

最適化の効率と精度を向上させるため、目的関数 $S^2(P)$ の勾配（各パラメータに関する偏微分）が解析的に計算され、jac 引数を通じて minimize 関数に提供されます。勾配 $\frac{\partial S^2}{\partial P_j}$ は以下のように計算されます。

\[ \frac{\partial S^2}{\partial P_j} = -2 \sum_{i=1}^{N} (y_d(x_i) - y_c(x_i; P)) \frac{\partial y_c(x_i; P)}{\partial P_j} \]

ここで、$\frac{\partial y_c}{\partial P_j}$ は各パラメータに対するモデル関数の偏微分であり、具体的には以下のようになります（単一ピークの場合）。

$X = \frac{C}{w}(x - x_0)$ とおくと、$y_c = I_0 \exp(-X^2)$。

$I_0$ に関する偏微分: $$ \frac{\partial y_c}{\partial I_0} = \exp(-X^2) $$
$x_0$ に関する偏微分: $$ \frac{\partial y_c}{\partial x_0} = 2 X \frac{C}{w} I_0 \exp(-X^2) $$
$w$ に関する偏微分: $$ \frac{\partial y_c}{\partial w} = \frac{2 X^2}{w} I_0 \exp(-X^2) $$

これらの解析的な勾配を利用することで、数値微分に伴う誤差を回避し、より安定した高速な収束が期待できます。

リアルタイムプロット

最適化の進行状況は、callback 関数を通じてリアルタイムで可視化されます。minimize 関数が各イテレーションでパラメータを更新するたびに callback 関数が呼び出され、現在のパラメータセットによるフィッティング結果と目的関数の値が Matplotlib グラフに更新表示されます。

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

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

NumPy: 数値計算を効率的に行うためのライブラリ。
SciPy: 科学技術計算のためのライブラリ。特に最適化機能 (scipy.optimize) を利用します。
Matplotlib: グラフ描画のためのライブラリ。リアルタイムプロットに使用します。

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

pip install numpy scipy matplotlib

必要な入力ファイル

本プログラムは、デフォルトでは外部からの入力ファイルを必要としません。フィッティングの対象となるデータ（yd）は、スクリプト内で定義された「真の」ピークパラメータ (peaks) を用いてガウス関数からシミュレートされます。フィッティングの初期推定値（x0）もスクリプト内で直接定義されますが、コマンドライン引数で上書きすることが可能です。

したがって、事前に用意すべきファイルはありません。

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

プログラムの実行後、以下の3つのCSVファイルが現在のディレクトリに生成されます。

initial.csv:
- 内容: フィッティング開始前のデータ。
- 列:
  - x: X座標
  - y(data): シミュレートされた（真の）観測データ
  - y(ini): 初期推定パラメータによるモデル予測値
final.csv:
- 内容: フィッティング完了後のデータ。
- 列:
  - x: X座標
  - y(data): シミュレートされた（真の）観測データ
  - y(ini): 初期推定パラメータによるモデル予測値
  - y(opt): 最適化された（フィッティング後の）パラメータによるモデル予測値
convergence.csv:
- 内容: 最適化の収束履歴。
- 列:
  - iter: 最適化のイテレーション回数
  - error: 各イテレーションにおける目的関数 $S^2$ の値

これらのファイルは、フィッティング結果の詳細な分析や、異なる初期値や設定でのフィッティング性能の比較に利用できます。

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

本プログラムは、コマンドライン引数でフィッティングの初期パラメータを指定できます。

基本的な実行コマンド:

python peakfit-scipy-minimize.py [initial_I0] [initial_x0] [initial_w]

引数の説明:

[initial_I0]: 初期推定のピーク強度 $I_0$ を指定します（浮動小数点数）。
[initial_x0]: 初期推定のピーク中心 $x_0$ を指定します（浮動小数点数）。
[initial_w]: 初期推定のピークの半値幅の半分 $w$ を指定します（浮動小数点数）。

これらの引数はオプションです。指定しない場合、プログラム内のデフォルト値 (x0 = [1.3, 0.6, 0.1]) が使用されます。

注意点: コマンドライン引数には iprint_interval という第4引数もコード上存在しますが、この変数はプログラムのどこにも使用されていないため、指定しても動作に影響はありません。

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

例1: デフォルトの初期パラメータで実行

python peakfit-scipy-minimize.py

実行結果:

プログラムが実行され、以下のような出力がコンソールに表示されます。（表示は scipy.optimize.minimize の設定によって異なりますが、収束情報や最終結果が含まれます）

Peak fitting by python scipy.optimize
Initial data are saved to [initial.csv]
callback 0: xk=[1.3 0.6 0.1]
   fmin=0.3341884488975945
callback 1: xk=[1.29956693 0.59969476 0.10014022]
   fmin=0.3341065184288673
... (最適化のイテレーションが続く) ...
callback X: xk=[1.10000000e+00 4.00000000e-01 4.00000000e-01]
   fmin=7.087455850901509e-28
       fun: 7.087455850901509e-28
       jac: [-3.36637153e-13 -8.08386348e-13  3.94821810e-13]
   message: 'Optimization terminated successfully.'
      nfev: 21
       nit: 14
      njev: 21
    status: 0
   success: True
         x: [1.10000000e+00 4.00000000e-01 4.00000000e-01]
Final data are saved to [final.csv]
Convergence data are saved to [convergence.csv]
Press ENTER to terminate:

同時に、リアルタイムプロットウィンドウが表示され、データ、初期フィット、最適化中のフィット（赤線）、および収束履歴が視覚的に確認できます。プログラムは最終的に「Press ENTER to terminate:」と表示され、Enterキーを押すまでグラフウィンドウを維持します。

この例では、初期パラメータ [1.3, 0.6, 0.1] から最適化が始まり、最終的に真のパラメータ [1.1, 0.4, 0.4] に非常に近い値に収束していることがわかります。

生成されるファイル：

initial.csv
final.csv
convergence.csv

例2: 初期パラメータを変更して実行

初期のピーク強度を 1.5、中心を 0.5、幅を 0.2 に変更して実行します。

python peakfit-scipy-minimize.py 1.5 0.5 0.2

実行結果:

コンソール出力は例1と同様に最適化のログが表示されますが、callback の最初の xk の値が [1.5 0.5 0.2] となり、そこから最適化が開始されます。

Peak fitting by python scipy.optimize
Initial data are saved to [initial.csv]
callback 0: xk=[1.5 0.5 0.2]
   fmin=0.5562772590483815
callback 1: xk=[1.49942767 0.49962772 0.20014819]
   fmin=0.5561339396328325
... (最適化が続く) ...
       fun: 7.087455850901509e-28
       jac: [-3.36637153e-13 -8.08386348e-13  3.94821810e-13]
   message: 'Optimization terminated successfully.'
      nfev: 21
       nit: 14
      njev: 21
    status: 0
   success: True
         x: [1.10000000e+00 4.00000000e-01 4.00000000e-01]
Final data are saved to [final.csv]
Convergence data are saved to [convergence.csv]
Press ENTER to terminate:

この場合も、最終的には同じ真のパラメータに近い値に収束することが期待されますが、初期値の違いにより収束経路やイテレーション回数が異なる可能性があります。

生成されるファイルは例1と同様です。initial.csv の y(ini) 列は、指定した初期パラメータに基づくデータとなります。