"""tknlsq.py
非線形最小二乗の小さな汎用ラッパー。
:doc:`tknlsq_usage`
設計方針
--------
- アプリ固有の物理モデルは外に置く
- ここでは residual_func(params) -> residual vector を最小化する
- 数値 Jacobian と cov ≈ s^2 (J^T J)^(-1) を共通化する
- 固定パラメータを扱いやすくする
SciPy がある場合は scipy.optimize.least_squares を使う。
SciPy がない環境では ImportError を出す。
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Callable, Dict, List, Mapping, Optional, Sequence, Tuple, Union
import math
import numpy as np
ArrayLike = Union[Sequence[float], np.ndarray]
ParamsLike = Union[Sequence[float], Mapping[str, float]]
[ドキュメント]
@dataclass
class NonlinearLSQResult:
"""非線形最小二乗の結果を格納するデータクラス。
非線形最小二乗フィッティングの実行結果、計算されたパラメータ、統計情報などを保持します。
:param params: 最適化された全パラメータ。入力が辞書形式だった場合は辞書、それ以外はNumPy配列。
:type params: Union[np.ndarray, Dict[str, float]]
:param params_free: 最適化された自由パラメータのみのNumPy配列。
:type params_free: np.ndarray
:param names: 全パラメータの名前のリスト。
:type names: List[str]
:param free_names: 自由パラメータの名前のリスト。
:type free_names: List[str]
:param fixed_names: 固定パラメータの名前のリスト。
:type fixed_names: List[str]
:param residuals: 最適化後の残差ベクトル。
:type residuals: np.ndarray
:param jacobian: 最適化された自由パラメータにおける残差関数のヤコビアン行列。None の場合もある。
:type jacobian: Optional[np.ndarray]
:param cov_free: 自由パラメータの共分散行列。None の場合もある。
:type cov_free: Optional[np.ndarray]
:param stderr_free: 自由パラメータの標準誤差のNumPy配列。None の場合もある。
:type stderr_free: Optional[np.ndarray]
:param stderr: 全パラメータの標準誤差。入力が辞書形式だった場合は辞書、それ以外はNumPy配列。固定パラメータは0.0。None の場合もある。
:type stderr: Union[Optional[np.ndarray], Dict[str, Optional[float]]]
:param sigma2_resid: 残差の分散 (s^2 = RSS / dof)。None の場合もある。
:type sigma2_resid: Optional[float]
:param dof: 残差の自由度 (N - p_free)。
:type dof: int
:param N: 残差ベクトルの要素数。
:type N: int
:param p_free: 自由パラメータの数。
:type p_free: int
:param RSS: 残差平方和 (Residual Sum of Squares)。
:type RSS: float
:param cost: scipy.optimize.least_squares から返されるコスト (0.5 * RSS)。
:type cost: float
:param success: 最適化が成功したかどうかを示す真偽値。
:type success: bool
:param message: 最適化の終了メッセージ。
:type message: str
:param nfev: 関数評価の回数。
:type nfev: Optional[int]
:param njev: ヤコビアン評価の回数。
:type njev: Optional[int]
:param status: scipy.optimize.least_squares からの終了ステータスコード。
:type status: Optional[int]
:param warning: 発生した可能性のある警告メッセージ。
:type warning: str
:param raw_result: scipy.optimize.least_squares から返される生の結果オブジェクト。
:type raw_result: object
"""
params: Union[np.ndarray, Dict[str, float]]
params_free: np.ndarray
names: List[str]
free_names: List[str]
fixed_names: List[str]
residuals: np.ndarray
jacobian: Optional[np.ndarray]
cov_free: Optional[np.ndarray]
stderr_free: Optional[np.ndarray]
stderr: Union[Optional[np.ndarray], Dict[str, Optional[float]]]
sigma2_resid: Optional[float]
dof: int
N: int
p_free: int
RSS: float
cost: float
success: bool
message: str
nfev: Optional[int] = None
njev: Optional[int] = None
status: Optional[int] = None
warning: str = ""
raw_result: object = field(default=None, repr=False)
@property
def sigma_resid(self) -> Optional[float]:
"""残差の標準偏差 (s)。
:returns: 残差の標準偏差。sigma2_residがNoneの場合はNone。
:rtype: Optional[float]
"""
if self.sigma2_resid is None:
return None
return math.sqrt(max(float(self.sigma2_resid), 0.0))
[ドキュメント]
def numerical_jacobian(
fun_vec: Callable[[np.ndarray], ArrayLike],
p: ArrayLike,
*,
rel_step: float = 1e-6,
abs_step: float = 1e-12,
) -> np.ndarray:
"""中心差分法を用いてベクトル関数のヤコビアンを数値的に計算する。
fun_vec(p) -> 残差ベクトル が与えられたとき、J[i, j] = d residual_i / d p_j となるヤコビアン行列を計算する。
微分のステップサイズは `rel_step * (abs(p[j]) + 1.0) + abs_step` で計算され、
これは相対ステップと絶対ステップを組み合わせた頑健な方法である。
:param fun_vec: パラメータベクトルを受け取り、残差ベクトルを返す関数。
:type fun_vec: Callable[[np.ndarray], ArrayLike]
:param p: ヤコビアンを評価するパラメータの点。
:type p: ArrayLike
:param rel_step: 相対ステップサイズ。パラメータ値の大きさに応じてステップサイズが決定される。
:type rel_step: float
:param abs_step: 最小絶対ステップサイズ。パラメータ値が0に近い場合でも最小限のステップを保証する。
:type abs_step: float
:returns: 計算されたヤコビアン行列。J[i, j] = d residual_i / d p_j。
:rtype: np.ndarray
:raises ValueError: fun_vecの出力サイズが数値微分中に変更された場合。
"""
p = np.asarray(p, dtype=float).reshape(-1)
f0 = np.asarray(fun_vec(p), dtype=float).reshape(-1)
J = np.zeros((f0.size, p.size), dtype=float)
for j in range(p.size):
dp = rel_step * (abs(p[j]) + 1.0) + abs_step
p_plus = p.copy(); p_plus[j] += dp
p_minus = p.copy(); p_minus[j] -= dp
f_plus = np.asarray(fun_vec(p_plus), dtype=float).reshape(-1)
f_minus = np.asarray(fun_vec(p_minus), dtype=float).reshape(-1)
if f_plus.size != f0.size or f_minus.size != f0.size:
raise ValueError("fun_vec output size changed during numerical differentiation")
J[:, j] = (f_plus - f_minus) / (2.0 * dp)
return J
[ドキュメント]
def covariance_from_jacobian(
residuals: ArrayLike,
J: np.ndarray,
*,
dof: Optional[int] = None,
rcond: float = 1e-12,
) -> Tuple[Optional[np.ndarray], Optional[np.ndarray], Optional[float], str]:
"""残差とヤコビアンからパラメータの共分散行列を近似的に計算する。
共分散は `cov ≈ s^2 * (J^T J)^(-1)` の形式で近似される。ここで `s^2 = RSS / dof`。
自由度 (dof) が0以下の場合、共分散と標準誤差は `None` を返す。
`J^T J` が特異行列の場合、擬似逆行列 (pinv) を使用する。
:param residuals: 最適化後の残差ベクトル。
:type residuals: ArrayLike
:param J: 残差関数のヤコビアン行列。
:type J: np.ndarray
:param dof: 残差の自由度 (N - パラメータ数)。None の場合、`residuals.size - J.shape[1]` で計算される。
:type dof: Optional[int]
:param rcond: 擬似逆行列計算時に使用される、特異値の相対的なしきい値。
:type rcond: float
:returns: 共分散行列、標準誤差ベクトル、残差分散 `s^2`、警告メッセージのタプル。
自由度が0以下の場合は (None, None, None, warning) が返される。
:rtype: Tuple[Optional[np.ndarray], Optional[np.ndarray], Optional[float], str]
:raises ValueError: Jが2Dでない場合、またはJの行数と残差のサイズが一致しない場合。
"""
r = np.asarray(residuals, dtype=float).reshape(-1)
J = np.asarray(J, dtype=float)
if J.ndim != 2:
raise ValueError("J must be 2D")
if J.shape[0] != r.size:
raise ValueError("J rows and residual size mismatch")
if dof is None:
dof = r.size - J.shape[1]
RSS = float(r.T @ r)
if dof <= 0:
return None, None, None, (
f"WARNING: residual degrees of freedom N-p = {r.size}-{J.shape[1]} = {dof}; "
"covariance and parameter errors are disabled."
)
sigma2 = RSS / dof
JTJ = J.T @ J
warning = ""
try:
cov0 = np.linalg.inv(JTJ)
except np.linalg.LinAlgError:
cov0 = np.linalg.pinv(JTJ, rcond=rcond)
warning = "WARNING: J^T J is singular; pseudo-inverse was used for covariance."
cov = sigma2 * cov0
stderr = np.sqrt(np.maximum(np.diag(cov), 0.0))
return cov, stderr, float(sigma2), warning
def _normalize_p0(p0: ParamsLike, names: Optional[Sequence[str]]) -> Tuple[np.ndarray, List[str], bool]:
"""初期パラメータ `p0` をNumPy配列形式に正規化する。
`p0` が辞書形式の場合、`names` が指定されていなければそのキーが名前となる。
`p0` がシーケンス形式の場合、`names` は必須で、パラメータ名と値の順序を決定する。
戻り値の `bool` は、入力 `p0` が辞書形式であったかどうかを示す。
:param p0: 初期パラメータ。シーケンス（リスト、NumPy配列など）または名前と値のマップ（辞書）。
:type p0: ParamsLike
:param names: パラメータ名のシーケンス。`p0` がシーケンスの場合に必須。
`p0` が辞書の場合、この引数はオプションで、キーの順序を決定するために使用できる。
:type names: Optional[Sequence[str]]
:returns: パラメータ値のNumPy配列、正規化されたパラメータ名のリスト、`p0` が辞書だったかどうかの真偽値。
:rtype: Tuple[np.ndarray, List[str], bool]
:raises ValueError: `names` の長さが `p0` の長さと一致しない場合。
"""
if isinstance(p0, Mapping):
if names is None:
names = list(p0.keys())
names = list(names)
vals = np.array([float(p0[n]) for n in names], dtype=float)
return vals, names, True
vals = np.asarray(p0, dtype=float).reshape(-1)
if names is None:
names = [f"p{i}" for i in range(vals.size)]
names = list(names)
if len(names) != vals.size:
raise ValueError("len(names) must match len(p0)")
return vals, names, False
def _normalize_fixed(
names: Sequence[str],
p0_full: np.ndarray,
fixed: Optional[Union[Sequence[str], Mapping[str, float]]],
) -> Tuple[np.ndarray, List[str], List[str], Dict[str, float]]:
"""固定パラメータの定義を正規化し、初期パラメータベクトルを調整する。
固定パラメータのリストまたは辞書を受け取り、正規化された固定パラメータ情報と、
固定値で更新された初期全パラメータベクトルを生成する。
:param names: 全パラメータの名前のシーケンス。
:type names: Sequence[str]
:param p0_full: 全パラメータの初期値を含むNumPy配列。
:type p0_full: np.ndarray
:param fixed: 固定するパラメータの定義。
パラメータ名のシーケンス（`p0_full` から値を取得）。
または、`{name: value}` の辞書（指定された値で固定）。
:type fixed: Optional[Union[Sequence[str], Mapping[str, float]]]
:returns: 固定値が適用された初期全パラメータ配列、自由パラメータ名のリスト、固定パラメータ名のリスト、
固定パラメータ名とその値の辞書。
:rtype: Tuple[np.ndarray, List[str], List[str], Dict[str, float]]
:raises ValueError: `fixed` に含まれるパラメータ名が `names` に存在しない場合。
"""
names = list(names)
fixed_values: Dict[str, float] = {}
if fixed is None:
fixed_names: List[str] = []
elif isinstance(fixed, Mapping):
fixed_names = list(fixed.keys())
fixed_values = {str(k): float(v) for k, v in fixed.items()}
else:
fixed_names = [str(x) for x in fixed]
fixed_values = {n: float(p0_full[names.index(n)]) for n in fixed_names}
for n in fixed_names:
if n not in names:
raise ValueError(f"fixed parameter {n!r} is not in names")
p0_mod = p0_full.copy()
for n, v in fixed_values.items():
p0_mod[names.index(n)] = v
free_names = [n for n in names if n not in fixed_names]
return p0_mod, free_names, fixed_names, fixed_values
def _build_full_vector(
p_free: np.ndarray,
p_base: np.ndarray,
names: Sequence[str],
free_names: Sequence[str],
) -> np.ndarray:
"""自由パラメータとベースベクトルから完全なパラメータベクトルを構築する。
`p_base` には固定パラメータの値が含まれ、`p_free` の値で自由パラメータの部分を上書きすることで、
完全なパラメータベクトルが再構築される。
:param p_free: 自由パラメータの現在の値を含むNumPy配列。
:type p_free: np.ndarray
:param p_base: 固定パラメータの値を含むベースの全パラメータ配列。自由パラメータの位置にはダミー値が含まれる。
:type p_base: np.ndarray
:param names: 全パラメータの名前のシーケンス。
:type names: Sequence[str]
:param free_names: 自由パラメータの名前のシーケンス。
:type free_names: Sequence[str]
:returns: 自由パラメータと固定パラメータの両方を含む完全なパラメータベクトル。
:rtype: np.ndarray
"""
p_full = np.asarray(p_base, dtype=float).copy()
for n, v in zip(free_names, p_free):
p_full[list(names).index(n)] = float(v)
return p_full
def _vector_to_output(p_full: np.ndarray, names: Sequence[str], as_mapping: bool):
"""完全なパラメータベクトルを指定された出力形式に変換する。
`as_mapping` が `True` の場合、パラメータ名と値の辞書に変換する。
`False` の場合、NumPy配列をコピーして返す。
:param p_full: 全パラメータのNumPy配列。
:type p_full: np.ndarray
:param names: 全パラメータの名前のシーケンス。
:type names: Sequence[str]
:param as_mapping: 出力形式を辞書にするかどうかの真偽値。
:type as_mapping: bool
:returns: 指定された形式に変換されたパラメータ値。
:rtype: Union[np.ndarray, Dict[str, float]]
"""
if as_mapping:
return {n: float(v) for n, v in zip(names, p_full)}
return p_full.copy()
def _stderr_to_output(
stderr_free: Optional[np.ndarray],
names: Sequence[str],
free_names: Sequence[str],
fixed_names: Sequence[str],
as_mapping: bool,
):
"""自由パラメータの標準誤差から、全パラメータの標準誤差を出力形式に変換する。
`as_mapping` が `True` の場合、パラメータ名と標準誤差の辞書に変換する。
`False` の場合、NumPy配列に変換する。
固定パラメータの標準誤差は `0.0` となる。自由パラメータの標準誤差が利用できない場合は `None` または `NaN` を使用する。
:param stderr_free: 自由パラメータの標準誤差を含むNumPy配列。標準誤差が計算されなかった場合はNone。
:type stderr_free: Optional[np.ndarray]
:param names: 全パラメータの名前のシーケンス。
:type names: Sequence[str]
:param free_names: 自由パラメータの名前のシーケンス。
:type free_names: Sequence[str]
:param fixed_names: 固定パラメータの名前のシーケンス。
:type fixed_names: Sequence[str]
:param as_mapping: 出力形式を辞書にするかどうかの真偽値。
:type as_mapping: bool
:returns: 指定された形式に変換された全パラメータの標準誤差。
:rtype: Union[Optional[np.ndarray], Dict[str, Optional[float]]]
"""
if as_mapping:
out: Dict[str, Optional[float]] = {n: None for n in names}
if stderr_free is not None:
for n, se in zip(free_names, stderr_free):
out[n] = float(se)
for n in fixed_names:
out[n] = 0.0
return out
if stderr_free is None:
return None
arr = np.full(len(names), np.nan, dtype=float)
for n, se in zip(free_names, stderr_free):
arr[list(names).index(n)] = float(se)
for n in fixed_names:
arr[list(names).index(n)] = 0.0
return arr
def _normalize_bounds(
bounds,
names: Sequence[str],
free_names: Sequence[str],
) -> Tuple[np.ndarray, np.ndarray]:
"""境界条件を正規化し、自由パラメータに対応する下限と上限の配列を返す。
`bounds` が `None` の場合、無限大の境界が設定される。
`bounds` が辞書の場合、`{name: (lower, upper)}` の形式で自由パラメータの境界が取得される。
`bounds` がシーケンスのタプル `(lower_bounds, upper_bounds)` の場合、全パラメータまたは自由パラメータの数に一致する必要がある。
:param bounds: 境界条件の定義。`None`、`{name: (lower, upper)}` の辞書、または `(lower_bounds, upper_bounds)` のタプル。
:type bounds: Optional[Union[Mapping[str, Tuple[float, float]], Tuple[ArrayLike, ArrayLike]]]
:param names: 全パラメータの名前のシーケンス。
:type names: Sequence[str]
:param free_names: 自由パラメータの名前のシーケンス。
:type free_names: Sequence[str]
:returns: 自由パラメータに対応する下限値のNumPy配列と上限値のNumPy配列のタプル。
:rtype: Tuple[np.ndarray, np.ndarray]
:raises ValueError: `bounds` の長さが全パラメータまたは自由パラメータのどちらとも一致しない場合。
"""
if bounds is None:
return np.full(len(free_names), -np.inf), np.full(len(free_names), np.inf)
if isinstance(bounds, Mapping):
lo = []
hi = []
for n in free_names:
b = bounds.get(n, (-np.inf, np.inf))
lo.append(float(b[0]))
hi.append(float(b[1]))
return np.asarray(lo, dtype=float), np.asarray(hi, dtype=float)
lo, hi = bounds
lo = np.asarray(lo, dtype=float).reshape(-1)
hi = np.asarray(hi, dtype=float).reshape(-1)
if lo.size == len(names) and hi.size == len(names):
idx = [list(names).index(n) for n in free_names]
return lo[idx], hi[idx]
if lo.size == len(free_names) and hi.size == len(free_names):
return lo, hi
raise ValueError("bounds must match either all parameters or free parameters")
[ドキュメント]
def nonlinear_lsq(
residual_func: Callable[[Union[np.ndarray, Dict[str, float]]], ArrayLike],
p0: ParamsLike,
*,
names: Optional[Sequence[str]] = None,
fixed: Optional[Union[Sequence[str], Mapping[str, float]]] = None,
bounds=None,
loss: str = "linear",
max_nfev: Optional[int] = None,
rel_step: float = 1e-6,
abs_step: float = 1e-12,
use_result_jac: bool = True,
scipy_kwargs: Optional[dict] = None,
) -> NonlinearLSQResult:
"""非線形最小二乗フィッティングを実行する。
`residual_func(params) -> 1D residual vector` を最小化する。
内部で `scipy.optimize.least_squares` を使用する。
`p0` が辞書の場合、`residual_func` には辞書が渡される。
`p0` がリスト/配列の場合、`residual_func` には `np.ndarray` が渡される。
:param residual_func: パラメータを受け取り、残差の1Dベクトルを返す関数。
受け取るパラメータの型は `p0` の形式によって異なる (NumPy配列または辞書)。
:type residual_func: Callable[[Union[np.ndarray, Dict[str, float]]], ArrayLike]
:param p0: 初期パラメータ。シーケンス（リスト、NumPy配列など）または名前と値のマップ（辞書）。
:type p0: ParamsLike
:param names: パラメータ名のシーケンス。`p0` がシーケンスの場合に必須。
`p0` が辞書の場合、この引数はオプションで、キーの順序を決定するために使用できる。
:type names: Optional[Sequence[str]]
:param fixed: 固定するパラメータの定義。
パラメータ名のシーケンス（`p0` から値を取得）。
または、`{name: value}` の辞書（指定された値で固定）。
:type fixed: Optional[Union[Sequence[str], Mapping[str, float]]]
:param bounds: パラメータの境界条件。
`None` (無制限)、`{name: (lower, upper)}` の辞書、
または `(lower_bounds, upper_bounds)` のタプル。
`lower_bounds` と `upper_bounds` は `p0` と同じ形式か、自由パラメータの数に合わせたシーケンス。
:type bounds: Optional[Union[Mapping[str, Tuple[float, float]], Tuple[ArrayLike, ArrayLike]]]
:param loss: `scipy.optimize.least_squares` に渡される損失関数 ('linear', 'soft_l1', 'huber', 'cauchy', 'arctan')。
:type loss: str
:param max_nfev: 関数評価の最大回数。`None` の場合、`least_squares` のデフォルトを使用。
:type max_nfev: Optional[int]
:param rel_step: 数値ヤコビアン計算の相対ステップサイズ。
:type rel_step: float
:param abs_step: 数値ヤコビアン計算の最小絶対ステップサイズ。
:type abs_step: float
:param use_result_jac: `scipy.optimize.least_squares` から返されるヤコビアンを使用するかどうか。
`True` の場合、`res.jac` があればそれを使用し、なければ `numerical_jacobian` をフォールバックする。
`False` の場合、常に `numerical_jacobian` を使用する。
:type use_result_jac: bool
:param scipy_kwargs: `scipy.optimize.least_squares` に渡す追加のキーワード引数辞書。
:type scipy_kwargs: Optional[dict]
:returns: 非線形最小二乗フィッティングの結果を格納した `NonlinearLSQResult` オブジェクト。
:rtype: NonlinearLSQResult
:raises ImportError: `scipy` がインストールされていない場合。
"""
try:
from scipy.optimize import least_squares
except Exception as exc: # pragma: no cover
raise ImportError("nonlinear_lsq requires scipy. Install with: pip install scipy") from exc
p0_full, names, as_mapping = _normalize_p0(p0, names)
p_base, free_names, fixed_names, _fixed_values = _normalize_fixed(names, p0_full, fixed)
p0_free = np.array([p_base[names.index(n)] for n in free_names], dtype=float)
lb, ub = _normalize_bounds(bounds, names, free_names)
def call_residual_full(p_full: np.ndarray) -> np.ndarray:
params = _vector_to_output(p_full, names, as_mapping)
r = np.asarray(residual_func(params), dtype=float).reshape(-1)
if not np.all(np.isfinite(r)):
# least_squares dislikes non-finite residuals. Replace by large residuals.
r = np.where(np.isfinite(r), r, 1e100)
return r
def call_residual_free(p_free: np.ndarray) -> np.ndarray:
p_full = _build_full_vector(p_free, p_base, names, free_names)
return call_residual_full(p_full)
if scipy_kwargs is None:
scipy_kwargs = {}
if len(free_names) == 0:
r = call_residual_full(p_base)
RSS = float(r.T @ r)
return NonlinearLSQResult(
params=_vector_to_output(p_base, names, as_mapping),
params_free=np.array([], dtype=float),
names=list(names),
free_names=[],
fixed_names=list(fixed_names),
residuals=r,
jacobian=None,
cov_free=None,
stderr_free=None,
stderr=_stderr_to_output(None, names, free_names, fixed_names, as_mapping),
sigma2_resid=None,
dof=r.size,
N=r.size,
p_free=0,
RSS=RSS,
cost=0.5 * RSS,
success=True,
message="all parameters are fixed",
warning="All parameters are fixed; optimization and covariance are skipped.",
raw_result=None,
)
res = least_squares(
call_residual_free,
p0_free,
bounds=(lb, ub),
loss=loss,
max_nfev=max_nfev,
**scipy_kwargs,
)
p_full = _build_full_vector(res.x, p_base, names, free_names)
residuals = call_residual_free(res.x)
N = residuals.size
p = len(free_names)
dof = N - p
RSS = float(residuals.T @ residuals)
if use_result_jac and getattr(res, "jac", None) is not None:
J = np.asarray(res.jac, dtype=float)
else:
J = numerical_jacobian(call_residual_free, res.x, rel_step=rel_step, abs_step=abs_step)
cov, stderr_free, sigma2, cov_warning = covariance_from_jacobian(residuals, J, dof=dof)
stderr = _stderr_to_output(stderr_free, names, free_names, fixed_names, as_mapping)
return NonlinearLSQResult(
params=_vector_to_output(p_full, names, as_mapping),
params_free=np.asarray(res.x, dtype=float),
names=list(names),
free_names=list(free_names),
fixed_names=list(fixed_names),
residuals=residuals,
jacobian=J,
cov_free=cov,
stderr_free=stderr_free,
stderr=stderr,
sigma2_resid=sigma2,
dof=dof,
N=N,
p_free=p,
RSS=RSS,
cost=float(res.cost),
success=bool(res.success),
message=str(res.message),
nfev=getattr(res, "nfev", None),
njev=getattr(res, "njev", None),
status=getattr(res, "status", None),
warning=cov_warning,
raw_result=res,
)
[ドキュメント]
def delta_method_variance(
output_func: Callable[[np.ndarray], ArrayLike],
p: ArrayLike,
cov: np.ndarray,
*,
rel_step: float = 1e-6,
abs_step: float = 1e-12,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""デルタ法を用いて、パラメータの共分散から出力関数の分散を計算する。
出力関数 `y = f(p)` の分散 `Var(y)` を、
`Var(y) ≈ J * Cov(p) * J^T` (ここで `J` は `f` のヤコビアン) として近似する。
:param output_func: パラメータベクトルを受け取り、出力量のベクトルを返す関数。
:type output_func: Callable[[np.ndarray], ArrayLike]
:param p: 出力を評価するパラメータの点。
:type p: ArrayLike
:param cov: パラメータ `p` の共分散行列。
:type cov: np.ndarray
:param rel_step: 数値ヤコビアン計算の相対ステップサイズ。
:type rel_step: float
:param abs_step: 数値ヤコビアン計算の最小絶対ステップサイズ。
:type abs_step: float
:returns: 出力値のベクトル `y0`、出力の分散のベクトル `var_y`、出力関数のヤコビアン `grads` のタプル。
:rtype: Tuple[np.ndarray, np.ndarray, np.ndarray]
:raises ValueError: `cov` の形状が `p` のサイズと一致しない場合。
"""
p = np.asarray(p, dtype=float).reshape(-1)
cov = np.asarray(cov, dtype=float)
if cov.shape != (p.size, p.size):
raise ValueError("cov shape must match p")
def f(v):
return np.asarray(output_func(v), dtype=float).reshape(-1)
y0 = f(p)
grads = numerical_jacobian(f, p, rel_step=rel_step, abs_step=abs_step)
var_y = np.einsum("ni,ij,nj->n", grads, cov, grads)
var_y = np.maximum(var_y, 0.0)
return y0, var_y, grads
[ドキュメント]
def delta_method_band(
output_func: Callable[[np.ndarray], ArrayLike],
p: ArrayLike,
cov: np.ndarray,
*,
nsigma: float = 1.0,
rel_step: float = 1e-6,
abs_step: float = 1e-12,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""デルタ法を用いて、出力関数 y=f(p) の値、信頼帯の下限、上限、標準偏差を計算する。
`y0 - nsigma * sigma_y` と `y0 + nsigma * sigma_y` で信頼帯を計算する。
信頼帯は `nsigma` に応じて調整される（例: `nsigma=1` で1シグマ帯）。
:param output_func: パラメータベクトルを受け取り、出力量のベクトルを返す関数。
:type output_func: Callable[[np.ndarray], ArrayLike]
:param p: 出力を評価するパラメータの点。
:type p: ArrayLike
:param cov: パラメータ `p` の共分散行列。
:type cov: np.ndarray
:param nsigma: 信頼帯の幅を決定するためのシグマ値の倍数。
:type nsigma: float
:param rel_step: 数値ヤコビアン計算の相対ステップサイズ。
:type rel_step: float
:param abs_step: 数値ヤコビアン計算の最小絶対ステップサイズ。
:type abs_step: float
:returns: 出力値 `y0`、信頼帯の下限 `y_low`、信頼帯の上限 `y_high`、出力の標準偏差 `sigma_y` のタプル。
:rtype: Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]
"""
y0, var_y, _grads = delta_method_variance(
output_func, p, cov, rel_step=rel_step, abs_step=abs_step
)
sigma = np.sqrt(var_y)
return y0, y0 - nsigma * sigma, y0 + nsigma * sigma, sigma
[ドキュメント]
def result_free_values(result: NonlinearLSQResult) -> np.ndarray:
"""`NonlinearLSQResult` オブジェクトから最適化された自由パラメータ値のNumPy配列を返す。
:param result: 非線形最小二乗フィッティングの結果オブジェクト。
:type result: NonlinearLSQResult
:returns: 最適化された自由パラメータ値のNumPy配列。
:rtype: np.ndarray
"""
return np.asarray(result.params_free, dtype=float)
[ドキュメント]
def result_free_cov(result: NonlinearLSQResult) -> Optional[np.ndarray]:
"""`NonlinearLSQResult` オブジェクトから自由パラメータの共分散行列を返す。
:param result: 非線形最小二乗フィッティングの結果オブジェクト。
:type result: NonlinearLSQResult
:returns: 自由パラメータの共分散行列。共分散が計算されなかった場合は `None`。
:rtype: Optional[np.ndarray]
"""
return result.cov_free