lsq_func_error_argparse.py 技術ドキュメント
プログラムの動作
lsq_func_error_argparse.py は、一般的な非線形最小二乗フィッティングを実行するためのPythonスクリプトです。このプログラムは、ユーザーが定義したモデル関数を使用して、観測データに最もよく適合するパラメータを特定します。主な機能と解決する課題は以下の通りです。
目的: 観測データに対して、ユーザーがPythonの文字列として定義した任意の非線形モデル関数を最小二乗法でフィッティングし、モデルパラメータの最適値、その標準誤差、およびモデル予測の信頼区間を推定すること。
主な機能:
柔軟なモデル関数定義: モデル関数はコマンドライン引数としてPythonの文字列で直接指定できます。これにより、様々な数学的モデル（例：ガウス関数、指数関数など）に柔軟に対応できます。
多変量対応: 複数の独立変数 (
x[0],x[1], ...) をサポートし、多次元データに対するフィッティングが可能です。パラメータ最適化:
scipy.optimize.minimizeを利用して、指定された最適化手法（例：BFGS、L-BFGS-B）でモデルパラメータを最適化します。パラメータ誤差の推定: 数値的なヤコビアン行列から、最適化されたパラメータの標準誤差を推定します。
信頼区間/不確実性帯の計算: デルタ法を用いて、モデル予測に対するパラメータ起因の不確実性帯および予測不確実性帯を計算します。これは、フィッティング結果の信頼性を視覚的に評価する上で重要です。
Excel出力: フィッティング結果、パラメータ情報、収束履歴、および計算グリッド上のモデル予測と不確実性帯をExcelファイルとして保存します。
リアルタイムプロット: フィッティングの進行状況をリアルタイムでグラフ表示し、収束状況を視覚的に監視できます。
シミュレーションモード: フィッティングを実行せず、指定された初期パラメータでモデル関数を評価するモードも提供します。
解決する課題:
特定の解析ソフトウェアに依存せず、Python環境でカスタムの非線形モデルフィッティングを行いたい研究者やエンジニア向けに、柔軟かつ強力なツールを提供します。
フィッティング結果だけでなく、パラメータの信頼性（標準誤差、共分散、相関）やモデル予測の不確実性（信頼区間）も定量的に評価することで、結果の解釈を深めます。
コマンドラインインターフェースにより、バッチ処理や他のスクリプトとの連携が容易になります。
原理
lsq_func_error_argparse.py は、非線形最小二乗フィッティングにおいて以下の数学的原理とアルゴリズムを使用しています。
1. 目的関数 (Objective Function)
プログラムは、観測データとモデル予測との間の残差二乗和の平均 (Mean Squared Error, MSE) を最小化することを目指します。 データ点 \(i\) における残差 \(r_i\) は、従属変数 \(y_i\) とモデル関数 \(f(\mathbf{x}_i, \mathbf{p})\) の差として定義されます。 $\( r_i = y_i - f(\mathbf{x}_i, \mathbf{p}) \)\( ここで \)\mathbf{x}_i\( は独立変数ベクトル、\) \mathbf{p} $ はモデルパラメータベクトルです。
目的関数 \(L(\mathbf{p})\) は、全データ点における残差二乗和の平均として与えられます。 $\( L(\mathbf{p}) = \frac{1}{N} \sum_{i=1}^{N} r_i^2 = \frac{1}{N} \mathbf{r}^T \mathbf{r} \)\( ここで \)N\( はデータ点の総数です。フィッティングは、この \)L(\mathbf{p})\( を最小化するパラメータ \)\mathbf{p}$ を見つけることを目的とします。
2. 最適化アルゴリズム
パラメータの最適化には、scipy.optimize.minimize 関数が使用されます。これは、様々な最適化手法（例：BFGS、L-BFGS-B、Nelder-Meadなど）をサポートしています。プログラムは、ユーザーが指定したメソッドに応じて、目的関数 \(L(\mathbf{p})\) を最小化する \(\mathbf{p}\) を探索します。
3. 数値微分とヤコビアン
勾配 (Gradient): 最適化アルゴリズムが勾配ベースの手法を使用する場合、目的関数 \(L(\mathbf{p})\) の勾配 \( \nabla L(\mathbf{p}) \) が必要になります。プログラムは、各パラメータ \(q_k\) (アクティブなパラメータ) に関する勾配を中央差分近似を用いて数値的に計算します。 $\( \frac{\partial L}{\partial q_k} \approx \frac{L(\mathbf{q} + h_k \mathbf{e}_k) - L(\mathbf{q} - h_k \mathbf{e}_k)}{2 h_k} \)\( ここで \)h_k\( は \)q_k\( に対する小さなステップサイズ、\) \mathbf{e}_k \( は \)k$ 番目の要素が1でそれ以外が0の単位ベクトルです。
モデルヤコビアン (Model Jacobian): パラメータの標準誤差や不確実性帯を推定するためには、モデル関数 \(f(\mathbf{x}, \mathbf{p})\) の各パラメータに関する偏微分からなるヤコビアン行列 \(J\) が必要です。この行列 \(J\) は \(N \times M\) の形状を持ち、その要素 \(J_{ij}\) は以下のように定義されます。 $\( J_{ij} = \frac{\partial f(\mathbf{x}_i, \mathbf{p})}{\partial p_j} \)\( プログラムはこれも中央差分近似により数値的に計算します。 \)\( \frac{\partial f(\mathbf{x}, \mathbf{p})}{\partial p_j} \approx \frac{f(\mathbf{x}, \mathbf{p} + h_j \mathbf{e}_j) - f(\mathbf{x}, \mathbf{p} - h_j \mathbf{e}_j)}{2 h_j} \)$
ステップサイズ \(h_k\) は、絶対ステップと相対ステップの組み合わせで決定されます。 $\( h_k = \text{jac\_absstep} + \text{jac\_relstep} \cdot \max(|\text{value}|, 1.0) \)$
4. パラメータ共分散行列と標準誤差
フィッティング後のパラメータ \(\mathbf{p}\) の不確実性を評価するために、共分散行列 \( \mathbf{C} \) が推定されます。非線形最小二乗法では、線形近似に基づいて共分散行列が計算されます。 $\( \mathbf{C} = \sigma_r^2 (J^T J)^{-1} \)$ ここで：
\(J\) は最適化されたパラメータでのモデルヤコビアン行列です。
\(J^T\) は \(J\) の転置行列です。
\((J^T J)^{-1}\) は \(J^T J\) の逆行列（特異値分解を用いた擬似逆行列
np.linalg.pinvが使用されます）です。\(\sigma_r^2\) は残差の分散であり、以下のように計算されます。 $\( \sigma_r^2 = \frac{\sum_{i=1}^{N} r_i^2}{N - M} = \frac{\mathbf{r}^T \mathbf{r}}{N - M} \)\( ここで \)N\( はデータ点の数、\)M\( は最適化されたアクティブなパラメータの数です。\)N - M$ は自由度 (Degrees of Freedom, DOF) です。
パラメータ \(p_j\) の標準誤差 \(\text{std}(p_j)\) は、共分散行列の対角要素の平方根として与えられます。 $\( \text{std}(p_j) = \sqrt{\mathbf{C}_{jj}} \)$
5. 信頼区間と不確実性帯 (Delta Method)
モデル予測に対する不確実性帯は、デルタ法 (Delta Method) を使用して計算されます。これは、関数の分散をその関数の一次テイラー近似を用いて推定する手法です。
パラメータ起因の不確実性 (Parameter Uncertainty): モデル関数 \(f(\mathbf{x}, \mathbf{p})\) の予測値の分散 \( \text{Var}[f(\mathbf{x}, \mathbf{p})] \) は、ヤコビアン \(J_{cal}\) とパラメータ共分散行列 \( \mathbf{C} \) を用いて以下のように近似されます。 $\( \text{Var}[f(\mathbf{x}_{cal}, \mathbf{p})] \approx J_{cal} \mathbf{C} J_{cal}^T \)\( ここで \)J_{cal}\( は計算グリッド \)x_{cal}\( におけるモデルヤコビアンです。 パラメータ起因の不確実性 \)\sigma_{param}(\mathbf{x}_{cal})\( は、この分散の平方根です。 \)\( \sigma_{param}(\mathbf{x}_{cal}) = \sqrt{\text{diag}(J_{cal} \mathbf{C} J_{cal}^T)} \)$
予測不確実性 (Prediction Uncertainty): 予測不確実性 \(\sigma_{pred}(\mathbf{x}_{cal})\) は、パラメータ起因の不確実性に加えて、観測データの測定誤差（残差分散 \(\sigma_r^2\) で代表される）も考慮に入れます。 $\( \sigma_{pred}(\mathbf{x}_{cal}) = \sqrt{\sigma_{param}(\mathbf{x}_{cal})^2 + \sigma_r^2} \)$
信頼区間の幅: 信頼区間 (Confidence Interval, CI) の幅は、Student の t 分布を用いて計算されます。指定された信頼水準
confidenceと自由度 \(N - M\) に基づく \(t\) 値 (\(t_{\alpha/2, N-M}\)) を不確実性 \(\sigma\) に乗じることで得られます。 $\( \text{CI} = t_{\alpha/2, N-M} \cdot \sigma \)\( ここで \)\alpha = 1 - \text{confidence}\( です。\)t_{\alpha/2, N-M}\( は、Student の t 分布の累積分布関数 (CDF) が \)0.5 + \text{confidence}/2.0$ となる点です。
必要な非標準ライブラリとインストール方法
このプログラムを実行するために必要な非標準ライブラリは以下の通りです。
numpy: 数値計算を効率的に行うためのライブラリ。
pandas: データ構造とデータ解析ツールを提供するライブラリ。主にデータの読み込みとExcel/CSV形式での出力に使用されます。
matplotlib: データのプロットやグラフ描画を行うためのライブラリ。フィッティング結果の可視化に使用されます。
scipy: 科学技術計算のためのライブラリ。特に
scipy.optimize.minimizeによる最適化、scipy.stats.tによるt分布の計算に使用されます。openpyxl: Pandas が Excel ファイル (.xlsx, .xlsm, .xls) を読み書きする際に内部的に使用されるライブラリ。
インストール方法
これらのライブラリは、Pythonのパッケージ管理システム pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。
pip install numpy pandas matplotlib scipy openpyxl
tklib (オプション)
ソースコード中には tklib.tkapplication と tklib.tkvariousdata のインポートが試みられていますが、これはGUIアプリケーションを構築するためのライブラリであり、本スクリプトのコマンドライン実行には必須ではありません。インポートに失敗してもプログラムは動作します。
必要な入力ファイル
lsq_func_error_argparse.py は、フィッティング対象となるデータを含むファイルを入力として受け取ります。
コマンドライン引数:
--infileで指定します。サポートされるファイル形式:
Excelファイル (
.xlsx,.xlsm,.xls)CSVファイル (
.csv)空白またはタブ区切りのテキストファイル (自動判別を試みます)
データ構造: 入力ファイルは表形式のデータである必要があります。
独立変数 (Independent Variables): 1つ以上のカラムが独立変数として使用されます。
--xlabels引数で、これらのカラム名（またはインデックス）をカンマ区切りで指定します。従属変数 (Dependent Variable): 1つのカラムが従属変数として使用されます。
--ylabel引数で、このカラム名（またはインデックス）を指定します。データの欠損値 (
NaN) は自動的に除外されます。また、--fit_range引数で指定された範囲外のデータも除外されます。
入力ファイルの例 (peak.xlsx)
以下の内容を持つ peak.xlsx ファイルを想定します。
1列目が独立変数 (X)、2列目が従属変数 (Y) です。
X |
Y |
|---|---|
0.1 |
0.05 |
0.2 |
0.15 |
0.3 |
0.3 |
0.4 |
0.6 |
0.5 |
0.9 |
0.6 |
1.0 |
0.7 |
0.8 |
0.8 |
0.4 |
0.9 |
0.2 |
1.0 |
0.08 |
このファイルを指定する場合、--xlabels "X" と --ylabel "Y" のように引数を指定します。カラム名ではなくインデックス（0から始まる）を使用する場合は、--xlabels "0" と --ylabel "1" のように指定することも可能です。
生成される出力ファイル
プログラムは、フィッティング結果を複数のExcelファイルとして生成します。出力ファイル名には、入力ファイル名から決定される接頭辞（--out_prefix で変更可能）が使用されます。
1. [out_prefix]-fit.xlsx
フィッティングデータ、計算グリッド上のモデル予測、収束履歴、およびフィッティングの要約が含まれます。
fit_dataシート:フィッティングに使用された生のデータ点 (
xlabelsで指定されたカラム)。従属変数 (
ylabelで指定されたカラム)。フィッティング開始時の初期パラメータでのモデル予測値 (
[ylabel](initial)列)。最適化後のパラメータでのモデル予測値 (
[ylabel](fit)列)。観測値とフィット値の残差 (
[ylabel](residual)列)。
calculation_gridシート:モデル予測と不確実性帯を計算するために使用されたグリッド点 (
[xlabels](cal)列)。単一独立変数モデルの場合はncalで指定された数の等間隔点、多変量モデルの場合は元のデータ点。フィットモデルの平均予測値 (
[ylabel](mean)列)。パラメータの不確実性に起因する標準偏差 (
[ylabel](sigma_param)列)。予測の不確実性（パラメータと残差の両方を含む）の標準偏差 (
[ylabel](sigma_pred)列)。パラメータの信頼区間の半幅 (
[ylabel](ci_param)列)。予測の信頼区間の半幅 (
[ylabel](ci_pred)列)。信頼区間計算に使用されたt分布の乗数 (
t_multiplier列)。
convergenceシート:最適化の反復回数 (
iter列)。各反復での平均二乗誤差 (MSE) (
MSE列)。
summaryシート:フィッティングに使用されたモデル関数 (
func)。最適化メソッド (
method)。データ点数 (
ndata)。最適化されたアクティブなパラメータの数 (
n_active_parameters)。自由度 (
dof)。残差二乗和 (RSS)。
平均二乗誤差 (MSE)。
残差の分散 (
sigma2_resid)。残差の標準偏差 (
sigma_resid)。信頼水準 (
confidence)。共分散推定のステータスメッセージ (
covariance_status)。
2. [out_prefix]-parameters.xlsx
フィッティングされたパラメータの詳細、共分散行列、相関行列、およびフィッティングの要約が含まれます。
parametersシート:パラメータのインデックス (
index列)。初期パラメータ値 (
initial列)。最適化された最終パラメータ値 (
final列)。パラメータの標準偏差 (
std列)。パラメータが固定されたかどうかを示すフラグ (
fixed列、0=False, 1=True)。
covarianceシート:最適化されたアクティブなパラメータ間の共分散行列。行と列のインデックスは
p0,p1, ... の形式。固定されたパラメータの行/列はNaNまたは0。
correlationシート:最適化されたアクティブなパラメータ間の相関行列。形式は共分散行列と同様。
summaryシート:[out_prefix]-fit.xlsxのsummaryシートと同じ内容。
3. [out_prefix]-convergence.xlsx
最適化の収束履歴のみを含む簡易的なExcelファイルです。内容は [out_prefix]-fit.xlsx の convergence シートと同じです。
コマンドラインでの使用例 (Usage)
lsq_func_error_argparse.py は、以下のコマンドライン引数を受け付けます。
usage: lsq_func_error_argparse.py [-h] [--mode {fit,sim}] [--infile INFILE] [--func FUNC] [--p0 P0] [--xlabels XLABELS] [--ylabel YLABEL] [--fit_range FIT_RANGE] [--method METHOD] [--maxiter MAXITER] [--tol TOL] [--bounds BOUNDS] [--fix FIX] [--jac_absstep JAC_ABSSTEP] [--jac_relstep JAC_RELSTEP] [--use_jac {0,1}] [--xcalmin XCALMIN] [--xcalmax XCALMAX] [--ncal NCAL] [--out_prefix OUT_PREFIX] [--xlsm_template XLSM_TEMPLATE] [--confidence CONFIDENCE] [--plot {0,1}] [--nplot_interval NPLOT_INTERVAL] [--plot_pause PLOT_PAUSE] [--plot_ci {0,1}] [--plot_sigma_param {0,1}] [--plot_sigma_pred {0,1}] [--plot_sigma_combined {0,1}] [--figsize FIGSIZE FIGSIZE] [--fontsize FONTSIZE] [--fontsize_legend FONTSIZE_LEGEND] [--pause {0,1}]
Generic nonlinear least-squares fitting with parameter errors and confidence bands
options:
-h, --help show this help message and exit
Basic IO and model:
--mode {fit,sim} fit or sim (default: fit)
--infile INFILE Input Excel/CSV/text file (default: peak.xlsx)
--func FUNC Model function. Use x[0], x[1], ... and p[0], p[1], ... (default: p[0]*exp(-((x[0]-p[1])/p[2])**2))
--p0 P0 Initial parameters, comma separated (default: 1.3,0.6,0.1)
--xlabels XLABELS Independent variable labels or column indices, comma separated (default: 0)
--ylabel YLABEL Dependent variable label or column index (default: 1)
--fit_range FIT_RANGE
Fit range for each x as xmin:xmax, comma separated. Example: -1:1,0:10 (default: *)
Optimizer:
--method METHOD scipy.optimize.minimize method (default: BFGS)
--maxiter MAXITER Maximum optimizer iterations (default: 1000)
--tol TOL Optimizer tolerance (default: 1e-08)
--bounds BOUNDS Bounds for parameters as lower:upper, comma separated. Use with L-BFGS-B/SLSQP/etc. (default: *)
--fix FIX Fixed parameter indices, comma separated. Example: 0,2 (default: )
Numerical derivatives:
--jac_absstep JAC_ABSSTEP
Absolute step for numerical derivative (default: 1e-08)
--jac_relstep JAC_RELSTEP
Relative step for numerical derivative (default: 1e-05)
--use_jac {0,1} Use numerical Jacobian for optimizer, 0/1 (default: 1)
Calculation grid and output:
--xcalmin XCALMIN Calculation x min for 1D plot/grid (default: *)
--xcalmax XCALMAX Calculation x max for 1D plot/grid (default: *)
--ncal NCAL Number of calculation grid points for 1D plot (default: 201)
--out_prefix OUT_PREFIX
Output prefix. Default: input file body (default: )
--xlsm_template XLSM_TEMPLATE
Reserved for compatibility (default: )
Uncertainty and plot options, Launcher friendly 0/1 flags:
--confidence CONFIDENCE
Confidence level. 0.6827 is about 1 sigma (default: 0.682689492)
--plot {0,1} Show plot, 0/1 (default: 1)
--nplot_interval NPLOT_INTERVAL
Update live fitting plot every N iterations. 0 disables live update (default: 10)
--plot_pause PLOT_PAUSE
Pause time for live plot update (default: 0.01)
--plot_ci {0,1} Plot uncertainty bands, 0/1 (default: 1)
--plot_sigma_param {0,1}
Plot parameter uncertainty band, 0/1 (default: 1)
--plot_sigma_pred {0,1}
Plot prediction uncertainty band, 0/1 (default: 1)
--plot_sigma_combined {0,1}
Plot param + measured-y scatter band, 0/1 (default: 0)
--figsize FIGSIZE FIGSIZE
Figure size (default: [10.0, 5.0])
--fontsize FONTSIZE Font size (default: 14)
--fontsize_legend FONTSIZE_LEGEND
Legend font size (default: 10)
--pause {0,1} Pause before termination, 0/1 (default: 1)
コマンドラインでの具体的な使用例
ここでは、前述の「必要な入力ファイル」で示した peak.xlsx を使用して、ガウス関数にフィッティングする例を示します。
peak.xlsx の内容（再掲）：
X |
Y |
|---|---|
0.1 |
0.05 |
0.2 |
0.15 |
0.3 |
0.3 |
0.4 |
0.6 |
0.5 |
0.9 |
0.6 |
1.0 |
0.7 |
0.8 |
0.8 |
0.4 |
0.9 |
0.2 |
1.0 |
0.08 |
実行コマンド:
python lsq_func_error_argparse.py \
--infile peak.xlsx \
--func "p[0]*exp(-((x[0]-p[1])/p[2])**2)" \
--p0 "1.0,0.6,0.1" \
--xlabels "X" \
--ylabel "Y" \
--confidence 0.95 \
--plot 1 \
--out_prefix my_peak_fit
コマンド引数の説明:
--infile peak.xlsx: 入力ファイルとしてpeak.xlsxを指定します。--func "p[0]*exp(-((x[0]-p[1])/p[2])**2)": フィッティングするガウス関数モデルを文字列で指定します。p[0]: 振幅p[1]: 中心p[2]: ピーク幅に関連するパラメータx[0]: 独立変数 (ここではXカラム)
--p0 "1.0,0.6,0.1": 初期パラメータp[0]=1.0,p[1]=0.6,p[2]=0.1をカンマ区切りで指定します。--xlabels "X": 独立変数のカラム名として "X" を指定します。--ylabel "Y": 従属変数のカラム名として "Y" を指定します。--confidence 0.95: 信頼区間の水準を95%に設定します。--plot 1: フィッティング後に結果のプロットを表示します。--out_prefix my_peak_fit: 出力ファイル名の接頭辞をmy_peak_fitに設定します。
実行結果:
標準出力: プログラムは、実行中に以下のような情報を標準出力に表示します（実際の数値は異なる場合があります）。
Generic nonlinear least-squares fitting infile : peak.xlsx func : p[0]*exp(-((x[0]-p[1])/p[2])**2) xlabels : ['X'] ylabel : Y method : BFGS p0 : [1. 0.6 0.1] fixed mask : [0 0 0] ndata : 10 fit ranges : [(-inf, inf)] Minimize: callback 1: MSE=8.70776953e-02 p=[1. 0.6 0.1 ] callback 2: MSE=2.97332213e-02 p=[1.06606018 0.59976361 0.08945625] ... (途中経過のMSEとパラメータ値が多数表示される) callback 15: MSE=1.04230191e-03 p=[1.00244675 0.60000000 0.10000000] Optimization result: success : True message : Optimization terminated successfully. nit : 14 MSE : 0.00104230191427 p : [1.00244675 0.6 0.1 ] Uncertainty estimate: status : OK ndata : 10 n_active_param: 3 dof : 7 RSS : 0.0104230191427 sigma_resid : 0.038596645362 Fitted parameters: p[0] = 1.002446754024 +- 0.0245089 p[1] = 0.600000000000 +- 0.0023249 p[2] = 0.100000000000 +- 0.0025988 Save fitting results to [my_peak_fit-fit.xlsx] Save parameter results to [my_peak_fit-parameters.xlsx] Save convergence history to [my_peak_fit-convergence.xlsx] Press ENTER to terminate >>
生成されるファイル: 指定した
--out_prefixに基づいて、以下のExcelファイルが生成されます。my_peak_fit-fit.xlsx: フィッティングデータ、計算グリッド、収束履歴、サマリー。my_peak_fit-parameters.xlsx: パラメータの詳細、共分散行列、相関行列、サマリー。my_peak_fit-convergence.xlsx: 収束履歴。
プロット:
--plot 1が指定されているため、以下の2つのグラフを含むMatplotlibウィンドウが表示されます。左側のグラフ (フィット結果):
元のデータ点が丸いマーカーで表示されます。
初期パラメータでのモデル関数が点線で表示されます。
最適化されたパラメータでのフィットモデルが実線で表示されます。
--confidence 0.95に基づく95%信頼区間の不確実性帯（パラメータ起因と予測の不確実性）が色付きの領域で表示されます。X軸には独立変数 "X" のラベル、Y軸には従属変数 "Y" のラベルが表示されます。
右側のグラフ (パラメータの誤差):
各パラメータの最適値が点とエラーバー（\(\pm 1\sigma\)）で表示されます。
X軸はパラメータのインデックス、Y軸はパラメータ値となります。
これらの出力により、フィッティングの品質、パラメータの信頼性、モデル予測の不確実性を総合的に評価できます。