"""tkfitdiag.py
フィット結果の診断ユーティリティ。
線形/非線形を問わず、パラメータ共分散行列が得られた後に使う。
主な用途:
- 共分散行列 -> 標準偏差・相関係数行列
- 共分散行列の固有値/固有ベクトルから、不確かな方向を調べる
- J^T J の固有値/条件数から、同定しにくいパラメータ結合を調べる
- 相対誤差・強相関・不確かな固有方向への寄与から固定候補を提案する
:doc:`tkfitdiag_usage`
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Sequence, Tuple
import math
import numpy as np
[ドキュメント]
@dataclass
class EigenSummary:
"""固有ベクトルの主要成分サマリ。
概要:
固有ベクトルの主要成分をまとめたデータ構造。
詳細説明:
固有値、その固有ベクトルに対応するランク、および主要なパラメータ成分を格納します。
:param rank: 固有値の順位 (1-indexed)。
:param eigenvalue: 固有値。
:param components: 主要なパラメータ成分とその寄与度（正規化されていない値）。
"""
rank: int
eigenvalue: float
components: List[Tuple[str, float]]
[ドキュメント]
@dataclass
class FixCandidate:
"""固定または外部拘束の候補。
概要:
パラメータの固定または外部拘束の候補を表すデータ構造。
詳細説明:
パラメータ名、スコア、値、標準誤差、相対誤差、および候補として提案された理由を保持します。
:param param: パラメータ名。
:param score: 固定候補としての総合スコア。
:param value: パラメータの値。
:param stderr: パラメータの標準誤差。
:param relerr: パラメータの相対標準誤差 (stderr / |value|)。値が0に近い場合はNoneになることがあります。
:param reasons: 候補として提案された理由のリスト。
"""
param: str
score: float
value: float
stderr: float
relerr: Optional[float]
reasons: List[str] = field(default_factory=list)
[ドキュメント]
@dataclass
class FitDiagnostics:
"""パラメータ共分散に基づく診断結果。
概要:
フィット診断の結果を格納するデータ構造。
詳細説明:
パラメータ名、値、誤差、共分散行列、相関係数行列、固有値解析の結果など、
フィットの不安定性や問題点を特定するための情報を含みます。
:param names: パラメータ名のリスト。
:param values: パラメータ値のNumpy配列。
:param stderr: パラメータの標準誤差のNumpy配列。
:param relerr: パラメータの相対標準誤差のNumpy配列。
:param cov: パラメータ共分散行列のNumpy配列。
:param corr: パラメータ相関係数行列のNumpy配列。
:param eig_cov_values_desc: 共分散行列の固有値（降順）のNumpy配列。
:param eig_cov_vectors_desc: 共分散行列の固有ベクトル（対応する固有値が降順）のNumpy配列。
:param eig_cov_summary: 共分散行列の主要な固有方向のサマリーリスト。
:param jtj: Jacobianの転置とJacobianの積 (J^T J) のNumpy配列。指定されない場合はNone。
:param eig_jtj_values_asc: J^T J の固有値（昇順）のNumpy配列。指定されない場合はNone。
:param eig_jtj_vectors_asc: J^T J の固有ベクトル（対応する固有値が昇順）のNumpy配列。指定されない場合はNone。
:param cond_jtj: J^T J の条件数。指定されない場合はNone。
:param warning: 診断中に発生した警告メッセージ。
"""
names: List[str]
values: np.ndarray
stderr: np.ndarray
relerr: np.ndarray
cov: np.ndarray
corr: np.ndarray
eig_cov_values_desc: np.ndarray
eig_cov_vectors_desc: np.ndarray
eig_cov_summary: List[EigenSummary]
jtj: Optional[np.ndarray] = None
eig_jtj_values_asc: Optional[np.ndarray] = None
eig_jtj_vectors_asc: Optional[np.ndarray] = None
cond_jtj: Optional[float] = None
warning: str = ""
[ドキュメント]
def covariance_to_correlation(cov: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
"""共分散行列から相関係数行列と標準偏差を計算する。
概要:
共分散行列から相関係数行列と標準偏差を計算します。
詳細説明:
共分散行列の対角要素の平方根から標準偏差を求め、それらを用いて相関係数行列を計算します。
標準偏差がゼロのパラメータが存在する場合、その相関係数はNaNとなります。
:param cov: 共分散行列。正方行列である必要があります。
:returns: 相関係数行列 (np.ndarray) と標準偏差 (np.ndarray) のタプル。
"""
cov = np.asarray(cov, dtype=float)
if cov.ndim != 2 or cov.shape[0] != cov.shape[1]:
raise ValueError("cov must be a square matrix")
stderr = np.sqrt(np.maximum(np.diag(cov), 0.0))
denom = np.outer(stderr, stderr)
with np.errstate(invalid="ignore", divide="ignore"):
corr = cov / denom
corr[denom == 0] = np.nan
np.fill_diagonal(corr, 1.0)
return corr, stderr
[ドキュメント]
def eigen_sorted_symmetric(A: np.ndarray, descending: bool = True) -> Tuple[np.ndarray, np.ndarray]:
"""対称行列の固有値・固有ベクトルをソートして返す。
概要:
対称行列の固有値と固有ベクトルを計算し、指定された順序でソートして返します。
詳細説明:
`numpy.linalg.eigh`を使用して固有値と固有ベクトルを計算した後、
`descending`引数に基づいて固有値をソートし、それに対応する固有ベクトルもソートします。
:param A: 対称行列。
:param descending: 固有値を降順にソートするかどうか。Trueの場合降順、Falseの場合昇順。デフォルトはTrue。
:returns: ソートされた固有値の配列 (np.ndarray) と、それに対応するソートされた固有ベクトルを行列として格納したタプル (np.ndarray)。
"""
A = np.asarray(A, dtype=float)
vals, vecs = np.linalg.eigh(A)
idx = np.argsort(vals)
if descending:
idx = idx[::-1]
return vals[idx], vecs[:, idx]
[ドキュメント]
def summarize_eigenvectors(
eigenvalues: Sequence[float],
eigenvectors: np.ndarray,
names: Sequence[str],
*,
topk: int = 3,
compk: int = 3,
) -> List[EigenSummary]:
"""固有ベクトルの主要成分を読みやすくまとめる。
概要:
固有値と固有ベクトルから、主要な固有方向とその構成要素を分かりやすい形式でサマライズします。
詳細説明:
与えられた固有値と固有ベクトルから、上位の`topk`個の固有方向と、
各固有方向において寄与の大きい`compk`個のパラメータを抽出し、
`EigenSummary`オブジェクトのリストとして生成します。
:param eigenvalues: 固有値のシーケンス。
:param eigenvectors: 固有ベクトルの行列。各列が1つの固有ベクトルに対応します。
:param names: パラメータ名のシーケンス。
:param topk: サマリーに含める上位の固有方向の数。デフォルトは3。
:param compk: 各固有方向で表示する主要なパラメータ成分の数。デフォルトは3。
:returns: 主要成分のサマリーを表す`EigenSummary`オブジェクトのリスト。
"""
names = list(names)
eigenvalues = np.asarray(eigenvalues, dtype=float)
eigenvectors = np.asarray(eigenvectors, dtype=float)
out: List[EigenSummary] = []
n_show = min(topk, len(eigenvalues))
for i in range(n_show):
v = eigenvectors[:, i]
order = np.argsort(np.abs(v))[::-1]
comps = [(names[j], float(v[j])) for j in order[: min(compk, len(order))]]
out.append(EigenSummary(rank=i + 1, eigenvalue=float(eigenvalues[i]), components=comps))
return out
[ドキュメント]
def diagnose_covariance(
names: Sequence[str],
values: Sequence[float],
cov: np.ndarray,
*,
jacobian: Optional[np.ndarray] = None,
topk_eigen: int = 3,
compk_eigen: int = 4,
) -> FitDiagnostics:
"""共分散行列からフィット診断を作る。
概要:
パラメータ共分散行列に基づいてフィットの診断結果を生成します。
詳細説明:
共分散行列から標準誤差、相対誤差、相関係数行列を計算し、さらに共分散行列の固有値解析を行います。
Jacobian行列が提供された場合は、J^T J行列の条件数と最小固有方向も評価し、
これらの診断情報を`FitDiagnostics`オブジェクトにまとめて返します。
:param names: パラメータ名のシーケンス。
:param values: パラメータ値のシーケンス。
:param cov: パラメータ共分散行列。
:param jacobian: 残差のJacobian行列。渡すと J^T J の条件数と最小固有方向も評価します。デフォルトはNone。
:param topk_eigen: 共分散行列の固有ベクトルサマリーで表示する上位の固有方向の数。デフォルトは3。
:param compk_eigen: 共分散行列の固有ベクトルサマリーで、各固有方向で表示する主要なパラメータ成分の数。デフォルトは4。
:returns: フィット診断結果を格納した`FitDiagnostics`オブジェクト。
"""
names = list(names)
values = np.asarray(values, dtype=float)
cov = np.asarray(cov, dtype=float)
if cov.shape != (len(names), len(names)):
raise ValueError("cov shape does not match names")
if values.size != len(names):
raise ValueError("values size does not match names")
corr, stderr = covariance_to_correlation(cov)
tiny = 1e-300
relerr = stderr / (np.abs(values) + tiny)
eig_vals, eig_vecs = eigen_sorted_symmetric(cov, descending=True)
eig_summary = summarize_eigenvectors(
eig_vals, eig_vecs, names, topk=topk_eigen, compk=compk_eigen
)
jtj = None
eig_jtj_vals = None
eig_jtj_vecs = None
cond_jtj = None
warning = ""
if jacobian is not None:
J = np.asarray(jacobian, dtype=float)
if J.ndim != 2 or J.shape[1] != len(names):
warning = "WARNING: jacobian shape does not match parameter count; J^T J skipped."
else:
jtj = J.T @ J
eig_jtj_vals, eig_jtj_vecs = eigen_sorted_symmetric(jtj, descending=False)
pos = eig_jtj_vals[eig_jtj_vals > 0]
if len(pos) == len(eig_jtj_vals):
cond_jtj = float(eig_jtj_vals[-1] / eig_jtj_vals[0])
else:
cond_jtj = math.inf
warning = "WARNING: J^T J has zero or negative eigenvalues; problem may be rank deficient."
return FitDiagnostics(
names=names,
values=values,
stderr=stderr,
relerr=relerr,
cov=cov,
corr=corr,
eig_cov_values_desc=eig_vals,
eig_cov_vectors_desc=eig_vecs,
eig_cov_summary=eig_summary,
jtj=jtj,
eig_jtj_values_asc=eig_jtj_vals,
eig_jtj_vectors_asc=eig_jtj_vecs,
cond_jtj=cond_jtj,
warning=warning,
)
[ドキュメント]
def propose_fix_candidates(
names: Sequence[str],
values: Sequence[float],
stderr: Sequence[float],
corr: np.ndarray,
*,
eig_cov_values: Optional[Sequence[float]] = None,
eig_cov_vectors: Optional[np.ndarray] = None,
corr_thr: float = 0.95,
relerr_thr: float = 0.5,
topn: int = 3,
) -> List[FixCandidate]:
"""固定/外部拘束候補をヒューリスティックに提案する。
概要:
フィットパラメータの固定または外部拘束の候補をヒューリスティックな基準で提案します。
詳細説明:
パラメータの相対標準誤差、他のパラメータとの強相関、共分散行列の最大固有値方向への寄与度に基づいて、
各パラメータの「固定しやすさ」スコアを計算します。
スコアが高く、相対誤差や強相関の閾値を超えるパラメータを候補としてリストアップし、
スコアの降順でソートして返します。
判断材料:
- 相対標準誤差 stderr / |value|
- 他パラメータとの強相関
- 共分散最大固有値方向への寄与
注意: これは数値的な「決まりにくさ」の診断であり、最終判断は物理的妥当性、
独立測定、文献値の有無と合わせて行う必要があります。
:param names: パラメータ名のシーケンス。
:param values: パラメータ値のシーケンス。
:param stderr: パラメータの標準誤差のシーケンス。
:param corr: 相関係数行列。
:param eig_cov_values: 共分散行列の固有値（降順）。Noneの場合、固有方向からの寄与は評価されません。デフォルトはNone。
:param eig_cov_vectors: 共分散行列の固有ベクトル（対応する固有値が降順）。Noneの場合、固有方向からの寄与は評価されません。デフォルトはNone。
:param corr_thr: 強相関と見なす相関係数の絶対値の閾値。デフォルトは0.95。
:param relerr_thr: 相対標準誤差が高いと見なす閾値。デフォルトは0.5。
:param topn: 提案する候補の最大数。デフォルトは3。
:returns: `FixCandidate`オブジェクトのリスト。スコアの高い順にソートされます。
"""
names = list(names)
values = np.asarray(values, dtype=float)
stderr = np.asarray(stderr, dtype=float)
corr = np.asarray(corr, dtype=float)
if len(names) != values.size or values.size != stderr.size:
raise ValueError("names, values, stderr sizes must match")
if corr.shape != (len(names), len(names)):
raise ValueError("corr shape does not match names")
tiny = 1e-300
relerr = stderr / (np.abs(values) + tiny)
score = np.nan_to_num(relerr, nan=0.0, posinf=1e300, neginf=0.0).astype(float)
reasons: Dict[str, List[str]] = {n: [] for n in names}
for i, n in enumerate(names):
if np.isfinite(relerr[i]):
if relerr[i] >= relerr_thr:
reasons[n].append(f"relative stderr = {relerr[i]:.3g} (>= {relerr_thr})")
else:
reasons[n].append(f"relative stderr = {relerr[i]:.3g}")
else:
reasons[n].append("relative stderr is not finite")
for i in range(len(names)):
for j in range(i + 1, len(names)):
cij = corr[i, j]
if not np.isfinite(cij):
continue
if abs(cij) >= corr_thr:
bonus = (abs(cij) - corr_thr) * 5.0 + 0.5
score[i] += bonus
score[j] += bonus
reasons[names[i]].append(f"strong corr with {names[j]}: {cij:+.3f} (>= {corr_thr})")
reasons[names[j]].append(f"strong corr with {names[i]}: {cij:+.3f} (>= {corr_thr})")
if eig_cov_values is not None and eig_cov_vectors is not None:
vecs = np.asarray(eig_cov_vectors, dtype=float)
if vecs.ndim == 2 and vecs.shape[0] == len(names) and vecs.shape[1] > 0:
v = vecs[:, 0]
w = np.abs(v) / (np.max(np.abs(v)) + tiny)
for i, n in enumerate(names):
if w[i] >= 0.5:
score[i] += 0.7 * w[i]
reasons[n].append(
f"dominant in most-uncertain eigen-direction: |v|/max={w[i]:.2f}"
)
items = [
FixCandidate(
param=names[i],
score=float(score[i]),
value=float(values[i]),
stderr=float(stderr[i]),
relerr=float(relerr[i]) if np.isfinite(relerr[i]) else None,
reasons=reasons[names[i]],
)
for i in range(len(names))
]
items.sort(key=lambda x: x.score, reverse=True)
filtered = [
item
for item in items
if (item.relerr is not None and item.relerr >= relerr_thr)
or any("strong corr" in r for r in item.reasons)
]
if not filtered:
filtered = items
return filtered[: min(topn, len(filtered))]
[ドキュメント]
def propose_fix_candidates_from_diagnostics(
diag: FitDiagnostics,
*,
corr_thr: float = 0.95,
relerr_thr: float = 0.5,
topn: int = 3,
) -> List[FixCandidate]:
"""FitDiagnostics から固定候補を提案するショートカット。
概要:
`FitDiagnostics`オブジェクトを利用して、パラメータの固定または外部拘束の候補を提案します。
詳細説明:
`diagnose_covariance`関数によって生成された`FitDiagnostics`オブジェクトから必要な情報を取り出し、
`propose_fix_candidates`関数を呼び出すための便利なショートカット関数です。
:param diag: フィット診断結果を格納した`FitDiagnostics`オブジェクト。
:param corr_thr: 強相関と見なす相関係数の絶対値の閾値。デフォルトは0.95。
:param relerr_thr: 相対標準誤差が高いと見なす閾値。デフォルトは0.5。
:param topn: 提案する候補の最大数。デフォルトは3。
:returns: `FixCandidate`オブジェクトのリスト。スコアの高い順にソートされます。
"""
return propose_fix_candidates(
diag.names,
diag.values,
diag.stderr,
diag.corr,
eig_cov_values=diag.eig_cov_values_desc,
eig_cov_vectors=diag.eig_cov_vectors_desc,
corr_thr=corr_thr,
relerr_thr=relerr_thr,
topn=topn,
)
[ドキュメント]
def diagnostics_to_dict(diag: FitDiagnostics) -> dict:
"""JSON保存しやすい辞書に変換する。
概要:
`FitDiagnostics`オブジェクトを、JSON形式で保存するのに適した辞書形式に変換します。
詳細説明:
Numpy配列はPythonのリストに変換され、`EigenSummary`オブジェクトも入れ子になった辞書に変換されます。
これにより、標準的なJSONシリアライザで`FitDiagnostics`オブジェクトの内容を簡単に保存・復元できます。
`jtj`や`eig_jtj`関連のフィールドがNoneの場合は、出力辞書でもNoneとして表現されます。
:param diag: 変換する`FitDiagnostics`オブジェクト。
:returns: 変換された辞書形式のデータ。
"""
return {
"names": diag.names,
"values": diag.values.tolist(),
"stderr": diag.stderr.tolist(),
"relerr": diag.relerr.tolist(),
"cov": diag.cov.tolist(),
"corr": diag.corr.tolist(),
"eig_cov": {
"eigenvalues_desc": diag.eig_cov_values_desc.tolist(),
"eigenvectors_desc": diag.eig_cov_vectors_desc.tolist(),
"summary": [
{"rank": s.rank, "eigenvalue": s.eigenvalue, "components": s.components}
for s in diag.eig_cov_summary
],
},
"eig_jtj": None
if diag.eig_jtj_values_asc is None
else {
"eigenvalues_asc": diag.eig_jtj_values_asc.tolist(),
"eigenvectors_asc": diag.eig_jtj_vectors_asc.tolist(),
"cond": diag.cond_jtj,
},
"warning": diag.warning,
}