adaptive_gaussian_ridge.py 技術ドキュメント

プログラムの動作

adaptive_gaussian_ridge.py は、スペクトルデータ（\(x-y\) データ）に対して、バックグラウンド成分と複数のガウス関数によるピーク成分を適応的にフィッティングするためのPythonプログラムです。特に、可変幅のガウス基底とチェビシェフ多項式によるバックグラウンドモデリングを組み合わせ、重み付きRidge正則化を用いて過学習を抑制します。オプションで、寄与の小さいガウス基底を段階的に除去する剪定（pruning）機能も備えています。

主な機能は以下の通りです。

  • Excelファイルから\(x-y\)データを読み込む、または内蔵のデモデータを生成して使用する。

  • \(y = \text{バックグラウンド} + \sum_i w_i \cdot \text{Gaussian}_i\) の形式でデータをフィッティングする。

  • バックグラウンド係数とガウス係数に異なるRidgeペナルティを適用できる。

  • ガウス基底の幅（\(\sigma\)）に応じてペナルティを調整し、幅の広い基底に強いペナルティを課すことができる。

  • オプションの剪定機能により、最も寄与の小さいガウス基底を段階的に削除し、最適なモデルを選択できる。

  • 計算結果をExcelファイルにエクスポートし、プロット用のデータを提供する。

    • 入力データ、フィッティングされた全スペクトル、バックグラウンドスペクトル、ガウス関数の合計スペクトル。

    • 基底パラメータ（中心、幅、係数）。

    • 剪定履歴（剪定モードの場合）。

    • 設定と結果の概要。

  • 結果の診断プロットを生成・表示・保存する。

このプログラムは、複雑なスペクトルデータからノイズを除去し、バックグラウンドとピークを分離し、各ピークの特性（位置、幅、強度）を定量化することを目的としています。

原理

adaptive_gaussian_ridge.py は、以下の数理的・アルゴリズム的原理に基づいてスペクトルデータの解析を行います。

1. モデル関数

プログラムは、入力データ \(y(x)\) をバックグラウンド成分とガウス関数の和で近似します。 $\( y(x) \approx f(x) = B(x) + \sum_{i=1}^{N_g} w_i G(x; \mu_i, \sigma_i) \)\( ここで、\)B(x)\( はバックグラウンド成分、\)G(x; \mu_i, \sigma_i)\( は\)i\(番目のガウス基底関数、\)\mu_i\( はその中心、\)w_i\( は係数、\)N_g$ はガウス基底の数です。

  • バックグラウンド \(B(x)\): チェビシェフ多項式 \(T_j(x')\) を用いてモデル化されます。 $\( B(x) = \sum_{j=0}^{\text{bg_order}} c_j T_j(x') \)\( \)x'\( は、元の \)x\( 座標をチェビシェフ多項式の定義域である \)[-1, 1]\( にスケーリングしたものです。`scale_x_to_chebyshev_domain` 関数によって行われます。 \)\( x' = \frac{2(x - x_{\min})}{x_{\max} - x_{\min}} - 1 \)\( ここで、\)c_j$ はチェビシェフ多項式の係数です。

  • ガウス基底 \(G(x; \mu_i, \sigma_i)\): 標準的なガウス関数が用いられます。 $\( G(x; \mu_i, \sigma_i) = \exp\left( - \frac{1}{2} \left( \frac{x - \mu_i}{\sigma_i} \right)^2 \right) \)$

2. ガウス基底の適応的な配置（中心 \(\mu_i\)）

ガウス基底の中心 \(\mu_i\) は、データの特性に応じて適応的に決定されます。adaptive_centers_from_derivatives 関数は、以下の手順で基底の中心を配置します。

  1. データの平滑化: Savitzky-Golayフィルター (scipy.signal.savgol_filter) を用いて、入力データ \(y\) を平滑化し、\(y_{smooth}\) を得ます。

  2. 微分の計算: 平滑化されたデータ \(y_{smooth}\) の1次微分 \(\frac{dy}{dx}\) と2次微分 \(\frac{d^2y}{dx^2}\) を計算します。

  3. 密度関数の構築: データの変化が大きい（勾配や曲率が大きい）領域に基底を密に配置するため、以下の密度関数 \(\rho(x)\) を定義します。 $\( \rho(x) = \epsilon + w_{grad} \left| \frac{dy_{smooth}}{dx} \right| + w_{curv} \left| \frac{d^2y_{smooth}}{dx^2} \right| \)\( ここで、\)\epsilon, w_{grad}, w_{curv}$ はそれぞれ基底の最小密度、勾配、曲率に対する重みを制御するパラメータです。

  4. \(x\) 軸の再マッピング: 密度関数を積分することで、新しい「仮想座標」\(u(x)\) を作成します。 $\( u(x) = \int_{x_{\min}}^x \rho(\xi) d\xi \)\( これにより、\)u$ 軸上ではデータが等間隔に分布するように見えます。

  5. 中心の選択: \(u\) 軸上で等間隔に \(N_{centers}\) 個の点を選択し、それを元の \(x\) 軸に逆変換することで、ガウス基底の中心 \(\mu_i\) を決定します。

3. ガウス基底の適応的な幅（\(\sigma_i\)）

ガウス基底の幅 \(\sigma_i\) は、estimate_variable_sigma 関数によって推定されます。複数の推定方法がサポートされており、--sigma-method 引数で選択できます。

  • distance メソッド: 各中心の周辺にある\(m\)個の隣接中心からの平均間隔に比例させて \(\sigma\) を決定します。 $\( \sigma_i = \alpha \cdot \frac{\mu_{i+m} - \mu_{i-m}}{2m} \)$

  • nearest メソッド: 各中心に最も近い左右の隣接中心との距離の最小値に比例させて \(\sigma\) を決定します。 $\( \sigma_i = \alpha \cdot \min(\mu_i - \mu_{i-1}, \mu_{i+1} - \mu_i) \)$

  • slope メソッド: 局所的な振幅と傾きに基づきます。 $\( \sigma_i = \alpha_{slope} \cdot \frac{A_i}{\left| \frac{dy_{smooth}}{dx}(\mu_i) \right| + \text{slope\_eps}} \)\( ここで、\)A_i\( は`estimate_local_amplitude`関数によって計算される、中心 \)\mu_i\( 周辺の平滑化されたデータ \)y_{smooth}$ の局所的な振幅（パーセンタイル幅または局所中央値からのずれ）です。これは、バックグラウンドの絶対値に依存しない頑健な振幅推定を目的としています。

  • hybrid_slope メソッド: nearest メソッドと slope メソッドで計算された \(\sigma\) のうち、小さい方を採用します。 $\( \sigma_i = \min(\sigma_{i, \text{nearest}}, \sigma_{i, \text{slope}}) \)\( 推定された \)\sigma_i$ は、--sigma-min および --sigma-max で指定された範囲にクリッピングされます。

4. 重み付きRidge回帰

モデルの係数（バックグラウンド係数 \(c_j\) とガウス係数 \(w_i\)）は、重み付きRidge回帰によって決定されます。fit_weighted_ridge 関数は、以下の最適化問題を解きます。 $\( \min_{\mathbf{\beta}} \left\| \mathbf{y} - \mathbf{X}\mathbf{\beta} \right\|^2 + \left\| \mathbf{P}\mathbf{\beta} \right\|^2 \)\( ここで、\)\mathbf{y}\( は入力データベクトル、\)\mathbf{X}\( はデザイン行列、\)\mathbf{\beta}\( は係数ベクトルです。\)\mathbf{P}$ はペナルティ項を定義する対角行列で、各係数に異なる重みを適用します。

  • デザイン行列 \(\mathbf{X}\): チェビシェフ多項式の基底とガウス基底を列とする行列です。 $\( \mathbf{X} = \begin{pmatrix} | & & | & & | & & | \\ T_0(\mathbf{x}') & \dots & T_{\text{bg_order}}(\mathbf{x}') & G(\mathbf{x}; \mu_1, \sigma_1) & \dots & G(\mathbf{x}; \mu_{N_g}, \sigma_{N_g}) \\ | & & | & & | & & | \end{pmatrix} \)$

  • ペナルティ行列 \(\mathbf{P}\): 対角要素 penalty_diag_ によって定義されます。

    • バックグラウンド係数には --lambda-bg で指定された定数ペナルティを適用。

    • ガウス係数には --lambda-gauss で指定された基本ペナルティに加えて、基底の幅 \(\sigma_i\) に応じた重み付けを適用します。 $\( \lambda_{gauss,i} = \lambda_{gauss} \cdot \left( \frac{\sigma_i}{\sigma_{ref}} \right)^{\text{sigma\_penalty\_power}} \)\( \)\sigma_{ref}$ はガウス基底の幅の中央値です。--sigma-penalty-power が正の場合、幅の広いガウス基底にはより強いペナルティが課され、より狭い基底が選択されやすくなります。

この最適化問題は、正規方程式 \(( \mathbf{X}^T \mathbf{X} + \mathbf{P}^T \mathbf{P} ) \mathbf{\beta} = \mathbf{X}^T \mathbf{y}\) を解くことで得られます。

5. 剪定（Pruning）

--mode prune が指定された場合、prune_basis_by_ridge_weight 関数は、最も寄与の小さいガウス基底を段階的に削除し、モデルの複雑さを軽減します。

  1. 初期フィッティング: まず、全ての初期ガウス基底とバックグラウンドを用いてRidge回帰を行います。

  2. 基底の除去: 各ガウス基底の係数の絶対値 (np.abs(model.coef_[n_bg:])) が小さい順に、--prune-frac で指定された割合の基底を削除します。このプロセスは、残りのガウス基底の数が --min-basis に達するまで繰り返されます。

  3. モデル評価: 各剪定ステップにおいて、フィットの品質をMAE、RMSE、AIC、BICなどの指標で評価します。有効自由度 (ridge_effective_df 関数) は、Ridge回帰におけるモデルの複雑さを表現するために計算され、AIC/BICの計算に利用されます。

  4. 最適モデルの選択: 最終的に、--metric（MAEまたはRMSE）を最小化し、かつ基底数が最小のモデルが最適モデルとして選択されます。--allowed-ratio は、最適な評価指標値に対して許容される劣化の割合を定義し、より少ない基底数を持つモデルを選ぶ柔軟性を提供します。

6. 評価指標

プログラムでは、モデルの性能を評価するために以下の指標が使用されます。

  • MAE (Mean Absolute Error): \( \frac{1}{N} \sum_{k=1}^N |y_k - \hat{y}_k| \)

  • RMSE (Root Mean Squared Error): \( \sqrt{\frac{1}{N} \sum_{k=1}^N (y_k - \hat{y}_k)^2} \)

  • RSS (Residual Sum of Squares): \( \sum_{k=1}^N (y_k - \hat{y}_k)^2 \)

  • AIC (Akaike Information Criterion): \( N \ln\left(\frac{RSS}{N}\right) + 2 \cdot \text{df} \)

  • BIC (Bayesian Information Criterion): \( N \ln\left(\frac{RSS}{N}\right) + \ln(N) \cdot \text{df} \) ここで、\(N\) はデータ点数、\(y_k\) は真の値、\(\hat{y}_k\) は予測値、df はモデルの有効自由度です。

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

このプログラムは、以下の非標準Pythonライブラリを使用します。

  • numpy

  • pandas

  • matplotlib

  • scipy

  • openpyxl (pandasがExcelファイルを読み書きするために内部的に使用します)

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

pip install numpy pandas matplotlib scipy openpyxl

必要な入力ファイル

プログラムは、以下のいずれかの形式で入力データを受け入れます。

  1. Excelファイル:

    • --infile <ファイル名.xlsx> でExcelファイルを指定します。

    • --sheet <シート名またはインデックス> で対象シートを指定します（デフォルト: 最初のシート、インデックス 0）。

    • --xcol <列名またはインデックス> で\(x\)軸データを含む列を指定します（デフォルト: 最初の列）。

    • --ycol <列名またはインデックス> で\(y\)軸データを含む列を指定します（デフォルト: 2番目の列）。

    • ファイル内には、数値データを含む列が必要です。欠損値（NaN）を含む行は自動的に削除され、\(x\)軸でソートされます。重複する\(x\)値がある場合は、平均\(y\)値が使用されます。

  2. デモデータ:

    • --infile demo を指定すると、プログラムに内蔵されたシミュレーションデータが生成され、入力として使用されます。

    • --n-data、--noise、--seed などの引数で、デモデータの特性を調整できます。

例: --infile data.xlsx --sheet Sheet1 --xcol Energy --ycol Intensity --infile data.xlsx --xcol 0 --ycol 1 (0番目の列がx、1番目の列がyの場合) --infile demo --n-data 200 --noise 0.05

生成される出力ファイル

プログラムは、計算結果を以下の形式で出力します。

  1. Excelファイル:

    • --outfile <ファイル名.xlsx> で出力ファイル名を指定します（指定しない場合、入力ファイル名から自動生成されます）。

    • このExcelファイルには、以下の複数のシートが含まれます。

      • InputData: 読み込まれた（または生成された）元の入力\(x, y\)データと、解析に使用されたソース列情報。

      • DataAndFit: 入力\(x\)座標における、入力\(y\)データ、平滑化された\(y\)データ、微係数、基底密度、フィッティング結果、残差データ。

      • FitSpectrum: フィッティングによって計算された、より高密度の\(x\)点での全フィットスペクトル、バックグラウンド、ガウス関数の合計スペクトル。デモデータの場合、真の関数も含まれます。

      • GaussianComponents: フィッティングされた個々のガウス関数の成分（バックグラウンドは含まない）。

      • BackgroundCoef: チェビシェフ多項式の係数とそのペナルティ値。

      • InitialGaussianBasis: 初期配置されたガウス基底の中心、幅、および推定方法ごとのシグマ候補値、最終的に選択されたかどうか。

      • SelectedGaussianBasis: 最終的に選択されたガウス基底の中心、幅、係数、絶対係数、ペナルティ、面積に類似する値。

      • PruneHistory (剪定モード (--mode prune) の場合のみ): 剪定ステップごとのガウス基底数、MAE、RMSE、RSS、有効自由度、AIC、BIC、そして選択されたステップを示す情報。

      • Summary: プログラムの実行パラメータ、入力データの概要、最終的なフィットの品質指標（MAE、RMSE、RSS）を含む概要情報。

  2. プロット画像ファイル:

    • --save-plot 1 オプションを指定した場合、--plotfile で指定されたファイル名（デフォルト: adaptive_gaussian_ridge_plot.png）で以下の3種類のプロット画像が保存されます。

      • *_diagnostics.png: データと平滑化、基底密度、\(\sigma\)候補、剪定履歴（エラー vs 基底数）。

      • *_mainfit.png: 入力データ、バックグラウンド、ガウス合計、全体フィット。

      • *_decomposition.png: 入力データ、全体フィット、ガウス合計、個々のガウス成分。

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

adaptive_gaussian_ridge.py は、多数の引数で動作をカスタマイズできます。

usage: adaptive_gaussian_ridge.py [-h] [--infile INFILE] [--sheet SHEET] [--xcol XCOL] [--ycol YCOL] [--outfile OUTFILE] [--mode {fit,prune}] [--n-data N_DATA] [--n-fit N_FIT] [--n-centers N_CENTERS] [--noise NOISE] [--seed SEED] [--bg-order BG_ORDER] [--lambda-bg LAMBDA_BG] [--lambda-gauss LAMBDA_GAUSS] [--sigma-penalty-power SIGMA_PENALTY_POWER] [--eps EPS] [--w-grad W_GRAD] [--w-curv W_CURV] [--sigma-method {distance,nearest,slope,hybrid_slope}] [--sigma-alpha SIGMA_ALPHA] [--sigma-m SIGMA_M] [--sigma-alpha-slope SIGMA_ALPHA_SLOPE] [--sigma-amp-window SIGMA_AMP_WINDOW] [--sigma-slope-eps SIGMA_SLOPE_EPS] [--sigma-min SIGMA_MIN] [--sigma-max SIGMA_MAX] [--smooth-window SMOOTH_WINDOW] [--smooth-order SMOOTH_ORDER] [--prune-frac PRUNE_FRAC] [--min-basis MIN_BASIS] [--metric {mae,rmse}] [--allowed-ratio ALLOWED_RATIO] [--show {0,1}] [--save-plot {0,1}] [--plotfile PLOTFILE]

Adaptive Gaussian basis + background + weighted Ridge + pruning for spectra

options:
  -h, --help            show this help message and exit
  --infile INFILE       Input Excel file. Use 'demo' for generated data.
  --sheet SHEET         Excel sheet name or index. Default: 0
  --xcol XCOL           x column name or zero-based column index. Default: first column
  --ycol YCOL           y column name or zero-based column index. Default: second column
  --outfile OUTFILE     Output Excel file. Default: created from infile name.
  --mode {fit,prune}    
  --n-data N_DATA       Number of generated data points for --infile demo.
  --n-fit N_FIT         Number of points in exported fitted spectra.
  --n-centers N_CENTERS
  --noise NOISE         
  --seed SEED           
  --bg-order BG_ORDER   Chebyshev background order. Use -1 for no background.
  --lambda-bg LAMBDA_BG 
  --lambda-gauss LAMBDA_GAUSS
  --sigma-penalty-power SIGMA_PENALTY_POWER
  --eps EPS             
  --w-grad W_GRAD       
  --w-curv W_CURV       
  --sigma-method {distance,nearest,slope,hybrid_slope}
                        Sigma estimation method. hybrid_slope uses min(nearest-distance, slope estimate).
  --sigma-alpha SIGMA_ALPHA
                        Distance-based sigma multiplier.
  --sigma-m SIGMA_M     Neighbor count for conventional distance sigma.
  --sigma-alpha-slope SIGMA_ALPHA_SLOPE
                        Slope-based sigma multiplier.
  --sigma-amp-window SIGMA_AMP_WINDOW
                        Window points for local amplitude estimation used in slope-based sigma.
  --sigma-slope-eps SIGMA_SLOPE_EPS
                        Small denominator for |y|/(|dy/dx|+eps) sigma estimate.
  --sigma-min SIGMA_MIN
  --sigma-max SIGMA_MAX
  --smooth-window SMOOTH_WINDOW
  --smooth-order SMOOTH_ORDER
  --prune-frac PRUNE_FRAC
  --min-basis MIN_BASIS
  --metric {mae,rmse}   
  --allowed-ratio ALLOWED_RATIO
  --show {0,1}          
  --save-plot {0,1}     
  --plotfile PLOTFILE

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

例1: デモデータを使用し、デフォルト設定で剪定を行う

これは、プログラムが提供する組み込みのデモデータを用いて、最も基本的な剪定モードを実行する例です。

python adaptive_gaussian_ridge.py --infile demo --mode prune

実行結果の説明:

  • demo_adaptive_gaussian_ridge.xlsx というExcelファイルが生成されます。これには、デモデータ、フィッティング結果、剪定履歴、選択されたガウス基底の情報などがシートごとに含まれます。

  • コンソールには剪定の履歴と最終的に選択されたモデルの概要が表示されます。

  • 3つの診断プロットウィンドウがポップアップ表示されます。これらのウィンドウを閉じるには、コンソールでEnterキーを押す必要があります。

例2: デモデータでガウス基底数やバックグラウンド次数を指定して剪定

デモデータに対して、初期ガウス基底数とチェビシェフバックグラウンドの次数を明示的に指定して剪定を行います。

python adaptive_gaussian_ridge.py --infile demo --mode prune --n-centers 60 --bg-order 3 --save-plot 1 --plotfile demo_output_plot

実行結果の説明:

  • 初期に60個のガウス基底を配置し、3次のチェビシェフ多項式でバックグラウンドをモデル化して剪定を実行します。

  • demo_adaptive_gaussian_ridge.xlsx が生成されます（既存の場合上書き）。

  • demo_output_plot_diagnostics.png、demo_output_plot_mainfit.png、demo_output_plot_decomposition.png の3つのPNG画像ファイルがカレントディレクトリに保存されます。

  • プロットは表示されず、コンソールに保存された旨のメッセージが表示されます。

例3: 実際のExcelファイルからデータを読み込み、バックグラウンドなしでフィッティング

data.xlsx というExcelファイルのSheet1からEnergy列を\(x\)、Intensity列を\(y\)として読み込み、バックグラウンドなし（bg-order -1）で、剪定を行わずに直接フィットします。

python adaptive_gaussian_ridge.py --infile data.xlsx --sheet Sheet1 --xcol Energy --ycol Intensity --bg-order -1 --mode fit --outfile result_no_bg.xlsx --show 0

実行結果の説明:

  • data.xlsxファイルから指定された列のデータを読み込みます。

  • バックグラウンドモデリングは行われず、ガウス基底のみでデータがフィッティングされます。

  • 剪定は行われず、初期に配置された全てのガウス基底でフィットが実行されます。

  • result_no_bg.xlsx というExcelファイルが生成されます。

  • プロットは表示されず（--show 0）、ファイルも保存されません。

例4: シグマの推定方法を変更し、幅の広いガウスに強いペナルティを適用して剪定

デモデータに対して、シグマの推定方法をhybrid_slopeに、ガウス基底のペナルティ重み付けをsigma-penalty-power 3.0に変更して剪定を実行します。

python adaptive_gaussian_ridge.py --infile demo --mode prune --sigma-method hybrid_slope --sigma-penalty-power 3.0 --lambda-gauss 1e-2 --prune-frac 0.05 --min-basis 3

実行結果の説明:

  • hybrid_slope メソッドによってガウス基底の幅 \(\sigma\) が推定されます。これは、近接距離に基づく\(\sigma\)と、平滑化されたデータの傾きに基づく\(\sigma\)の小さい方を選択するものです。

  • sigma-penalty-power 3.0 の指定により、\(\sigma\)が大きくなるにつれてペナルティが指数的に強くなります。これは、より狭いピークを優先的に選択する傾向に繋がります。

  • ガウスのRidgeペナルティ (lambda-gauss) がデフォルトの1.0e-3から1.0e-2に強化されます。

  • 剪定の割合が5%に、最小基底数が3に設定されます。

  • demo_adaptive_gaussian_ridge.xlsx が生成され、コンソールに剪定結果と概要が出力されます。

  • 3つの診断プロットウィンドウがポップアップ表示されます。