はじめに

論文レベルのスペクトル解析では、ベストフィット値そのものよりも「その誤差をどう評価しているか」 が常に問われます。
フィット結果の表には当たり前のように誤差が並びますが、その裏側で 実際にどのような手順で誤差が計算されているのかは、意外と明示的に説明されることが少なく、初学者にとっては大きな障壁になりがちです。

特に pyXSPEC を用いた解析では、

XSPEC 対話モードとの対応関係が掴みにくい
error コマンドをどのタイミングで、どのパラメータに対して回しているのかが見えにくい
free / frozen / linked なパラメータの整理が煩雑になりやすい

といった点から、誤差評価の流れがブラックボックス化しやすいと感じています。

本記事の対象読者は、XSPEC の対話モードではなく、Python から操作する pyXSPEC を使って解析している人です。
というのも、私は XSPEC の対話主体の解析に慣れているユーザーではなく、実際の解析作業の多くを pyXSPEC で行っているためです。

XSPEC の基本操作や対話的な解析フローについては、山田さんや榎戸さんをはじめ、すでに多くの優れた解説記事がありますので、本記事ではそれらを改めて説明することはしません。

本記事では、

論文で一般的に行われている誤差評価を、pyXSPEC ではどのように実装しているのか

という点に絞り、一研究者の立場から、ペルセウス座銀河団の実データを用いて、その流れを一緒に整理します。
物理的解釈を主目的とするものではなく、解析作業を安定して進めるための実装寄りのメモとして活用してもらえれば幸いです。

本記事で使用するデータ

デモとして、
XRISM / Resolve ペルセウス座銀河団 Early Release Data
を使用します。

「とりあえず誰でも試せる」ことを重視してこのデータを選んでいますが、次の点を強調しておきます。

私はペルセウス座銀河団の解析を専門としていません
使用するモデルに物理的根拠を与えるつもりはありません
このモデルを「正しい」「推奨される」ものだとは考えていません

モデルはあくまで、誤差計算の流れを確認するための例です。

データ解析

ここでは、パワーロー成分に複数のガウシアン成分を加え、
ガウシアン成分の幅（Sigma）を共通とする制約を入れたモデルを用います。

繰り返しになりますが、このモデルは 物理的に正しいことを主張するものではありません。
あくまで、誤差計算の流れを確認するためのデモです。

import xspec as xs
from xspec import Xset
import matplotlib.pyplot as plt
import numpy as np

# XSPECのデータをクリア
xs.AllData.clear()

# スペクトルデータの読み込み
spec = xs.Spectrum('xa_merged_p0px1000_Hp.pi', respFile='xa_merged_p0px1000_HpS.rmf', arfFile='rsl_standard_GVclosed.arf')

# 指定したエネルギー範囲を無視
spec.ignore('**-6.46 6.65-**')

# モデルの定義
model = xs.Model('powerlaw + gaussian + gaussian + gaussian + gaussian + gaussian')
print('--------------------------------------------------------------')
print('Model component names:', model.componentNames)
print('--------------------------------------------------------------')

# パワーローモデルの初期パラメータ設定
model.powerlaw.PhoIndex = 2.0
model.powerlaw.norm = 1

# ガウシアンの初期パラメータ設定
model.gaussian.LineE = 6.521
model.gaussian.Sigma = 0.001
model.gaussian.norm = 1e-2
model.gaussian.LineE.frozen = True

model.gaussian_3.LineE = 6.535
model.gaussian_3.Sigma = 0.001
model.gaussian_3.norm = 1e-2
model.gaussian_3.LineE.frozen = True

model.gaussian_4.LineE = 6.549
model.gaussian_4.Sigma = 0.001
model.gaussian_4.norm = 1e-2
model.gaussian_4.LineE.frozen = True

model.gaussian_5.LineE = 6.565
model.gaussian_5.Sigma = 0.001
model.gaussian_5.norm = 1e-2
model.gaussian_5.LineE.frozen = True

model.gaussian_6.LineE = 6.583
model.gaussian_6.Sigma = 0.001
model.gaussian_6.norm = 2e-2
model.gaussian_6.LineE.frozen = True

# シグマパラメータをリンク
model.gaussian_3.Sigma.link = model.gaussian.Sigma
model.gaussian_4.Sigma.link = model.gaussian.Sigma
model.gaussian_5.Sigma.link = model.gaussian.Sigma
model.gaussian_6.Sigma.link = model.gaussian.Sigma

# フィットの実行（1回目）
xs.Fit.query = 'yes'
xs.Fit.perform()

# ラインエネルギーをフリーにして再度フィット
model.gaussian.LineE.frozen = False
model.gaussian_3.LineE.frozen = False
model.gaussian_4.LineE.frozen = False
model.gaussian_5.LineE.frozen = False

# フィットの実行（2回目）
xs.Fit.query = 'yes'
xs.Fit.perform()

# モデルを保存
Xset.save('best_fit_para.xcm')

# Matplotlibでデータをプロット
xs.Plot.device = '/null'
xs.Plot.xAxis = 'keV'
xs.Plot.add = True
xs.Plot('data delc')

# XSPECのプロットデータを取得
x_data = np.array(xs.Plot.x(1,1))
y_data = np.array(xs.Plot.y(1,1))
x_err = np.array(xs.Plot.xErr(1,1))
y_err = np.array(xs.Plot.yErr(1,1))

x_delc = np.array(xs.Plot.x(1,2))
y_delc = np.array(xs.Plot.y(1,2))
x_delc_err = np.array(xs.Plot.xErr(1,2))
y_delc_err = np.array(xs.Plot.yErr(1,2))

# ラベル情報取得
data_label = xs.Plot.labels(1)
delc_label = xs.Plot.labels(2)

y_model = np.array(xs.Plot.model())

# Matplotlibでプロット
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True, gridspec_kw={'height_ratios': [3, 1]})
fig.subplots_adjust(hspace=0.05)

# 観測データとモデルをプロット
ax1.errorbar(x_data, y_data, xerr=x_err, yerr=y_err, fmt='none', ecolor='k', capsize=0, label='data')
ax1.plot(x_data, y_model, color='r', label='model')

# 各コンポーネントをプロット
num_components = xs.Plot.nAddComps()
component_names = model.componentNames
for i in range(1, num_components + 1):
    y_component = np.array(xs.Plot.addComp(i))
    component_name = component_names[i - 1]
    ax1.plot(x_data, y_component, ':', label=component_name)

ax1.set_ylabel(data_label[1])
ax1.legend()

# 残差プロット
ax2.errorbar(x_delc, y_delc, xerr=x_delc_err, yerr=y_delc_err, fmt='none', ecolor='k', capsize=0, label='Data')
ax2.axhline(0, color='lime')
ax2.set_ylim(-5, 5)
ax2.set_xlabel(delc_label[0])
ax2.set_ylabel(delc_label[1])

plt.savefig('best_fit_para.png', bbox_inches='tight')
plt.show()

========================================================================
Model powerlaw<1> + gaussian<2> + gaussian<3> + gaussian<4> + gaussian<5> + gaussian<6> Source No.: 1   Active/On
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   powerlaw   PhoIndex            7.79771      +/-  3.98282      
   2    1   powerlaw   norm                3113.74      +/-  2.34025E+04  
   3    2   gaussian   LineE      keV      6.52100      +/-  3.86499E-03  
   4    2   gaussian   Sigma      keV      4.61565E-03  +/-  6.28477E-05  
   5    2   gaussian   norm                1.14767E-04  +/-  3.16128E-06  
   6    3   gaussian   LineE      keV      6.53500      +/-  3.19922E-03  
   7    3   gaussian   Sigma      keV      4.61565E-03  = p4
   8    3   gaussian   norm                5.01745E-05  +/-  2.51390E-06  
   9    4   gaussian   LineE      keV      6.54900      +/-  1.36247E-03  
  10    4   gaussian   Sigma      keV      4.61565E-03  = p4
  11    4   gaussian   norm                1.08535E-04  +/-  3.07518E-06  
  12    5   gaussian   LineE      keV      6.56500      +/-  2.08983E-03  
  13    5   gaussian   Sigma      keV      4.61565E-03  = p4
  14    5   gaussian   norm                7.75569E-05  +/-  2.69482E-06  
  15    6   gaussian   LineE      keV      6.58300      frozen
  16    6   gaussian   Sigma      keV      4.61565E-03  = p4
  17    6   gaussian   norm                2.79545E-04  +/-  5.54831E-06  
________________________________________________________________________


Fit statistic  : Chi-Squared                  543.99     using 379 bins.

Test statistic : Chi-Squared                  543.99     using 379 bins.
 Null hypothesis probability of 5.05e-09 with 367 degrees of freedom

フィットが終了すると、このような画面が表示されます。
ただし、ここに表示されているエラーは参考値に近いものであり、必ずしも厳密な誤差評価とは限りません。その点について次節で簡単に説明します。

なぜ追加の誤差計算が必要なのか

XSPEC のフィット結果として表示される誤差は、多くの場合 対称誤差です。
しかし実際には、

パラメータ空間の歪み
境界への衝突
他パラメータとの強い相関

などにより、誤差は非対称になることがほとんどです。

そのため、論文レベルの解析では error コマンドによる信頼区間評価 を行うのが一般的です。

誤差を計算する

誤差計算には error コマンドを用います。
XSPEC の対話モードでは、

error 4

のように指定したパラメータ番号の誤差を計算します。

このとき XSPEC は、対象パラメータを動かしつつ他のパラメータを再最適化し、
指定した Δstat に達するまで探索を行います。

理論背景については、以下の記事が非常に参考になります。
https://qiita.com/yamadasuzaku/items/3eac6df7d4ddbb55b0c0

ここで得られる誤差はデフォルトで 90% 信頼区間 であることを押さえておきましょう。

pyXSPEC で誤差計算を自動化する

モデルが複雑になってくると、

どのパラメータが free なのか
link されているパラメータはどれか

を 目視で確認するのはかなり手間になります。
特にガウシアン成分が多い場合や、複数の制約を入れている場合は、パラメータ番号を追うだけでも大変です。

そこで本記事では、以下の処理を pyXSPEC 上で自動化します。

.xcm ファイルを読み込む
free かつ unlinked なパラメータのみを自動抽出する
Fit.error をまとめて実行する
誤差計算結果を CSV として保存する

以下のスクリプトでは、事前に保存しておいた .xcm ファイルを読み込み、
誤差計算 → 結果整理 → 保存までを一括で行います。

import xspec as xs
from xspec import Xset, Fit
import pandas as pd
import numpy as np
import argparse


# 凍結されておらず、リンクされていないパラメータのみ取得する関数
def get_free_parameters_with_errors(model, free_param_indices):
    free_params = {}
    errors = {}

    for component in model.componentNames:
        comp = getattr(model, component)  # コンポーネント取得
        for param_name in comp.parameterNames:
            param = getattr(comp, param_name)  # パラメータ取得
            if not param.frozen and param.link == "":  # frozenされておらず、リンクされていない場合のみ
                param_key = f"{component}_{param_name}"
                free_params[param_key] = param.values[0]  # パラメータ名と値を保存
                error_values = param.error
                errors[param_key] = error_values

    return free_params, errors

# フィッティングを実行したファイル名を指定
filename_xcm = f"best_fit_para.xcm"

# 誤差計算した際の出力ファイルを指定
update_filename_xcm = f"best_fit_para_update.xcm"
param_csv = f"best_fit_para_update_error.csv"

# データとモデルを読み込む
Xset.restore(filename_xcm)
model = xs.AllModels(1)

# モデルフィットを実行
print("Fitting model...")
Fit.perform()

# 凍結されていないパラメータのインデックスを収集
free_param_indices = []
for component in model.componentNames:
    comp = getattr(model, component)
    for param_name in comp.parameterNames:
        param = getattr(comp, param_name)
        if not param.frozen and param.link == "":
            free_param_indices.append(param.index)

# エラーを一度に推定 (例: "2 3 6" の形式で実行)
if free_param_indices:
    index_string = " ".join(map(str, free_param_indices))
    print(f"Calculating errors for parameters: {index_string}")
    xs.Fit.error(f"{index_string}")

# パラメータとエラーを取得
free_params, errors = get_free_parameters_with_errors(model, free_param_indices)

# 結果を保存用のリストに追加
results = []
for param_name, value in free_params.items():
    if param_name in errors:
        error_values = errors[param_name]
        results.append({
            "parameter": param_name,
            "value": value,
            "error -": value - error_values[0],
            "error +": error_values[1] - value,
            "error flag": error_values[2],
            "stat_method": Fit.statMethod,
            "stat": Fit.statistic,
            "dof": Fit.dof
        })

# パラメータが更新されることもあるので最新を保存
Xset.save(update_filename_xcm, info='a')

# DataFrame 化して CSV に保存
df = pd.DataFrame(results)

# DataFrame の確認用の表示
df_text = df.to_string(index=False)
print(df_text)

df.to_csv(param_csv, index=False, encoding="utf-8")
print(f"Error results saved to {param_csv}")

実行時の注意点

スクリプトを実行すると、途中で

Number of trials exceeded: continue fitting? y

というメッセージが表示されることがあります。

これは 誤差探索がやや難しい状態にあることを示していますが、
必ずしも計算が失敗しているわけではありません。

今回のように挙動を確認する目的であれば、y を入力して処理を続行して問題ありません。

出力結果について

処理が完了すると、各パラメータについて

最適値
非対称誤差（−側 / ＋側）
error flag
使用した統計量（statistic）
stat 値と自由度

が、次のように一覧で出力されます。

parameter           value    error -   error +  error flag stat_method       stat  dof
powerlaw_PhoIndex   ...      ...        ...      FFFFFFFFF         chi ...
powerlaw_norm       ...      ...        ...      FTFFFFFFF         chi ...
...

同じ内容がCSV ファイルとしても保存されるため、そのまま

論文用の表作成
解析ログの管理
後続解析との比較

に利用できます。
（なお、LaTeX で非対称誤差（±）をきれいに表として整形するには、また少し壁があります。私自身もまだ整理しきれていないため、本記事では割愛します。）

出力結果の例と見方

実行すると、以下のような結果が得られました。

出力例

        parameter       value    error -      error + error flag stat_method       stat  dof
powerlaw_PhoIndex    7.788876   0.021710     0.023755  FFFFFFFFF         chi 535.580579  367
    powerlaw_norm 3085.772545 104.407305 19167.732710  FTFFFFFFF         chi 535.580579  367
   gaussian_LineE    6.521194   0.000228     0.000281  FFFFFFFFF         chi 535.580579  367
   gaussian_Sigma    0.004616   0.000121     0.000079  FFFFFFFFF         chi 535.580579  367
    gaussian_norm    0.000118   0.000007     0.000004  FFFFFFFFF         chi 535.580579  367
 gaussian_3_LineE    6.535831   0.000614     0.000529  FFFFFFFFF         chi 535.580579  367
  gaussian_3_norm    0.000050   0.000003     0.000006  FFFFFFFFF         chi 535.580579  367
 gaussian_4_LineE    6.549462   0.000220     0.000374  FFFFFFFFF         chi 535.580579  367
  gaussian_4_norm    0.000108   0.000005     0.000005  FFFFFFFFF         chi 535.580579  367
 gaussian_5_LineE    6.565367   0.000297     0.000503  FFFFFFFFF         chi 535.580579  367
  gaussian_5_norm    0.000077   0.000005     0.000004  FFFFFFFFF         chi 535.580579  367
  gaussian_6_norm    0.000277   0.000006     0.000008  FFFFFFFFF         chi 535.580579  367

error flag の確認が最重要

まず最初に見るべきなのは error flag です。

すべて F → 誤差計算は比較的きれいに収束している
T が含まれる → 何らかの問題が発生している

エラーフラグの意味は公式ドキュメントにまとめられています。
https://heasarc.gsfc.nasa.gov/docs/software/xspec/manual/node60.html#tclouterror

9 文字のフラグは、以下の条件に対応しています。

新しい最小値が見つかった
非単調性が検出された
最小化が問題に遭遇した可能性
下限に達した
上限に達した
パラメータが凍結されている
−方向の探索に失敗
＋方向の探索に失敗
reduced χ² が大きすぎる

例えば FFFFFFFFT は
「reduced χ² が高すぎるため誤差評価が失敗した」ことを意味します。

今回の例では powerlaw_norm に T が含まれており、
このパラメータの誤差は あまり信用できない可能性があることが分かります。

説明記事として「誤差計算がきれいに収束していない例」を示すことに違和感を覚える方もいるかもしれません。しかし実際の解析では、モデルの選び方、エネルギー範囲の設定、ビンまとめの方法によって、誤差評価は驚くほどセンシティブに振る舞います。そのため、error flag を手がかりに、次にどのパラメータや解析条件を見直すべきかを、物理的な背景も含めて検討していくことが重要になります。

おわりに

本記事では、
「pyXSPEC で論文レベルの誤差評価をどう回すか」
という一点に絞って整理しました。

物理の議論をあえて避け、実装と作業フローだけに集中することで、
「この誤差は、こういう手順で回して付けている」
と、実務として説明できる状態を確認することを目指しました。

pyXSPEC を使って解析を進めている方が、
誤差計算に過度な時間を取られることなく、
自分の解析手順の一部として納得しながら扱えるようになるための、
ひとつの足場になれば幸いです。

関連記事