"""tklsq.py
汎用線形回帰・誤差評価・モデル選択ユーティリティ。
.. _tklsq_usage:
関連リンク:
:doc:`tklsq_usage`
設計方針
--------
- 基本は線形最小二乗: y = X beta + eps
- パラメータ誤差は cov(beta) = sigma^2 (X^T W X)^(-1)
- 推定値誤差は Var[y_mean] = diag(X_new cov(beta) X_new^T)
- 予測誤差は Var[y_pred] = Var[y_mean] + sigma^2
- モデル選択は AIC/BIC/AICc と、必要に応じてベイズ evidence を使う
このファイルは tklib / matplotlib / openpyxl に依存しない「計算コア」です。
入出力、グラフ、Excel保存は呼び出し側スクリプトに残す想定です。
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Dict, Iterable, List, Mapping, Optional, Sequence, Tuple, Union
import math
import numpy as np
ArrayLike = Union[Sequence[float], np.ndarray]
CandidateType = Union[Sequence[np.ndarray], Mapping[str, np.ndarray]]
# =====================================================================
# Dataclasses
# =====================================================================
[ドキュメント]
@dataclass
class LeastSquaresResult:
"""線形最小二乗の結果を格納するデータクラス。
線形最小二乗回帰の計算結果と、その統計的な情報、誤差推定の結果を保持します。
:param beta: np.ndarray: 推定された回帰係数ベクトル。
:param beta_std: Optional[np.ndarray]: 回帰係数の標準誤差。誤差推定が無効な場合はNone。
:param cov_beta: Optional[np.ndarray]: 回帰係数の共分散行列。誤差推定が無効な場合はNone。
:param sigma2_resid: Optional[float]: 残差分散 (sigma^2)。誤差推定が無効な場合はNone。
:param residuals: np.ndarray: 残差ベクトル (y - y_fit)。
:param y_fit: np.ndarray: モデルによる推定値 (X @ beta)。
:param N: int: データ点数。
:param p: int: 回帰係数の数（モデルパラメータの数）。
:param dof: int: 残差自由度 (N - p)。
:param rank: int: 設計行列のランク。
:param RSS: float: 残差平方和 (Residual Sum of Squares)。
:param WRSS: float: 重み付き残差平方和 (Weighted Residual Sum of Squares)。
:param singular_values: np.ndarray: 設計行列の特異値。
:param condition_number: float: 設計行列の条件数。
:param error_estimation_enabled: bool: 誤差推定が有効かどうか。
:param warning: str: 発生した警告メッセージ。デフォルトは空文字列。
:param known_sigma: bool: ノイズの分散が既知として扱われたかどうか。デフォルトはFalse。
"""
beta: np.ndarray
beta_std: Optional[np.ndarray]
cov_beta: Optional[np.ndarray]
sigma2_resid: Optional[float]
residuals: np.ndarray
y_fit: np.ndarray
N: int
p: int
dof: int
rank: int
RSS: float
WRSS: float
singular_values: np.ndarray
condition_number: float
error_estimation_enabled: bool
warning: str = ""
known_sigma: bool = False
@property
def sigma_resid(self) -> Optional[float]:
"""残差の標準偏差 (sigma_resid) を返すプロパティ。"""
if self.sigma2_resid is None:
return None
return float(math.sqrt(max(self.sigma2_resid, 0.0)))
[ドキュメント]
@dataclass
class PredictionResult:
"""線形モデルの予測結果と誤差バンドを格納するデータクラス。
新しい入力データに対するモデルの推定値と、その予測誤差の範囲を保持します。
:param y_mean: np.ndarray: 予測される応答変数の平均値。
:param sigma_param: Optional[np.ndarray]: パラメータ推定の不確かさに起因する標準誤差。
:param sigma_pred: Optional[np.ndarray]: 予測の総標準誤差（パラメータ不確かさ＋ノイズ）。
:param var_param: Optional[np.ndarray]: パラメータ推定の不確かさに起因する分散。
:param var_pred: Optional[np.ndarray]: 予測の総分散（パラメータ不確かさ＋ノイズ）。
"""
y_mean: np.ndarray
sigma_param: Optional[np.ndarray]
sigma_pred: Optional[np.ndarray]
var_param: Optional[np.ndarray]
var_pred: Optional[np.ndarray]
[ドキュメント]
@dataclass
class ModelSelectionResult:
"""複数モデルの比較結果を格納するデータクラス。
異なるモデル候補を比較し、選択基準（AIC、BIC、Evidenceなど）に基づいて評価した結果を保持します。
:param labels: List[str]: 各モデルのラベル。
:param scores: np.ndarray: 各モデルのスコア（例: AIC, BIC, log_evidence）。
:param weights: np.ndarray: 各モデルの相対的な重み（モデル確率）。
:param best_index: int: 最良モデルのインデックス。
:param criterion: str: モデル選択に使用された基準の名称。
:param results: Optional[List[LeastSquaresResult]]: 各モデルの線形最小二乗結果のリスト。デフォルトはNone。
:param log_evidences: Optional[np.ndarray]: 各モデルの対数エビデンス。デフォルトはNone。
"""
labels: List[str]
scores: np.ndarray
weights: np.ndarray
best_index: int
criterion: str
results: Optional[List[LeastSquaresResult]] = None
log_evidences: Optional[np.ndarray] = None
@property
def best_label(self) -> str:
"""最良モデルのラベルを返すプロパティ。"""
return self.labels[self.best_index]
[ドキュメント]
@dataclass
class BayesianLinearResult:
"""ベイズ線形回帰の事後分布の結果を格納するデータクラス。
ベイズ線形回帰によって得られたパラメータの事後平均、共分散、および関連するハイパーパラメータを保持します。
:param mean: np.ndarray: 事後分布の平均（推定された回帰係数）。
:param cov: np.ndarray: 事後分布の共分散行列。
:param alpha: float: パラメータの事前分布の精度（逆分散）ハイパーパラメータ。
:param beta: float: ノイズの精度（逆分散）ハイパーパラメータ。
:param sigma_noise: float: 推定されたノイズの標準偏差 (1/sqrt(beta))。
:param log_evidence: Optional[float]: モデルの対数エビデンス。デフォルトはNone。
:param n_iter: int: ハイパーパラメータ推定の繰り返し回数。デフォルトは0。
:param converged: bool: ハイパーパラメータ推定が収束したかどうか。デフォルトはFalse。
"""
mean: np.ndarray
cov: np.ndarray
alpha: float
beta: float
sigma_noise: float
log_evidence: Optional[float] = None
n_iter: int = 0
converged: bool = False
# =====================================================================
# Design matrices
# =====================================================================
[ドキュメント]
def polynomial_design_matrix(
x: ArrayLike,
order: Optional[int] = None,
basis_indices: Optional[Sequence[int]] = None,
) -> np.ndarray:
"""多項式基底に基づく設計行列を作成します。
与えられた1次元の説明変数 `x` に対して、指定された次数または基底インデックスの
多項式項を持つ設計行列 `X` を生成します。
:param x: ArrayLike: 1次元の説明変数データ。
:param order: Optional[int]: 0から `order` までのすべてのべき指数を持つ多項式基底を使用する場合の最大次数。
`basis_indices` が指定された場合は無視されます。デフォルトはNone。
:param basis_indices: Optional[Sequence[int]]: 使用する多項式のべき指数をリストで指定します。
例: [0, 1, 3] は定数項、1次項、3次項を意味します。
デフォルトはNone。
:returns: np.ndarray: 多項式基底によって構築された設計行列 `X` (shape: N, p)。
:raises ValueError: `order` と `basis_indices` の両方がNoneの場合。
"""
x = np.asarray(x, dtype=float).reshape(-1)
if basis_indices is None:
if order is None:
raise ValueError("Specify either order or basis_indices")
basis_indices = list(range(order + 1))
return np.vstack([x ** int(i) for i in basis_indices]).T
[ドキュメント]
def design_matrix_from_functions(
x: ArrayLike,
basis_functions: Sequence,
) -> np.ndarray:
"""任意の基底関数リストから設計行列を作成します。
与えられた1次元の説明変数 `x` と基底関数のシーケンスに基づいて、
設計行列 `X` を構築します。各基底関数は `x_array` または `x_scalar` を引数に取り、
対応する基底ベクトルの要素を返すと想定されます。
:param x: ArrayLike: 1次元の説明変数データ。
:param basis_functions: Sequence: 基底関数のリスト。各関数は `f(x_array)` を返すか、
`np.vectorize` 可能な `f(x_scalar)` であると想定されます。
:returns: np.ndarray: 基底関数によって構築された設計行列 `X` (shape: N, p)。
"""
x = np.asarray(x, dtype=float).reshape(-1)
cols = []
for f in basis_functions:
v = np.asarray(f(x), dtype=float)
if v.shape == () or v.shape != x.shape:
v = np.asarray([f(xi) for xi in x], dtype=float)
cols.append(v.reshape(-1))
return np.vstack(cols).T
[ドキュメント]
def add_intercept(X: np.ndarray) -> np.ndarray:
"""既存の設計行列の先頭列に定数項（切片）を追加します。
入力の設計行列 `X` の各行の先頭に1の列を追加し、定数項を含むモデルに対応させます。
:param X: np.ndarray: 既存の設計行列。
:returns: np.ndarray: 定数項が追加された新しい設計行列。
"""
X = np.asarray(X, dtype=float)
return np.column_stack([np.ones(X.shape[0]), X])
# =====================================================================
# Least squares and uncertainty
# =====================================================================
[ドキュメント]
def linear_lsq(
X: np.ndarray,
y: ArrayLike,
*,
weights: Optional[ArrayLike] = None,
known_sigma: bool = False,
rcond: Optional[float] = None,
cond_warning: float = 1.0e12,
) -> LeastSquaresResult:
"""線形最小二乗回帰を実行し、回帰係数と統計的な誤差情報を計算します。
モデル `y = X @ beta + eps` に基づいて回帰係数 `beta` を推定します。
設計行列がフルランクであり、残差自由度が正の場合、係数の共分散行列、
標準誤差、および残差分散も計算されます。重み付き最小二乗もサポートします。
:param X: np.ndarray: 設計行列 (shape: N, p)。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param weights: Optional[ArrayLike]: 各データ点に適用する重みベクトル。通常は `1 / sigma_y^2` を指定します。
Noneの場合、等重み（重み1）と見なされます。デフォルトはNone。
:param known_sigma: bool: Trueの場合、`weights` が絶対的な `1/sigma_y^2` であるとみなし、
`cov(beta) = (X^T W X)^(-1)` を使用します。
Falseの場合、残差から `sigma^2` を推定し、それを `(X^T W X)^(-1)` に乗じます。
デフォルトはFalse。
:param rcond: Optional[float]: `np.linalg.lstsq` における特異値のカットオフ値。デフォルトはNone。
:param cond_warning: float: 設計行列の条件数がこの値を超えた場合に、結果の`warning`フィールドに
警告メッセージを追記します。デフォルトは1.0e12。
:returns: LeastSquaresResult: 線形最小二乗回帰の結果を格納するデータクラス。
:raises ValueError: `X` が2次元配列でない場合、または `X` の行数と `y` の要素数が一致しない場合、
あるいは `weights` のサイズと `y` の要素数が一致しない場合、または `weights` が負の値を含む場合。
"""
X = np.asarray(X, dtype=float)
y = np.asarray(y, dtype=float).reshape(-1)
if X.ndim != 2:
raise ValueError("X must be a 2D array")
if X.shape[0] != y.size:
raise ValueError(f"X rows ({X.shape[0]}) and y size ({y.size}) mismatch")
N, p = X.shape
dof = N - p
if weights is None:
w = np.ones(N, dtype=float)
else:
w = np.asarray(weights, dtype=float).reshape(-1)
if w.size != N:
raise ValueError(f"weights size ({w.size}) and y size ({N}) mismatch")
if np.any(w < 0):
raise ValueError("weights must be non-negative")
sw = np.sqrt(w)
Xw = X * sw[:, None]
yw = y * sw
beta, _residuals_lstsq, rank, svals = np.linalg.lstsq(Xw, yw, rcond=rcond)
y_fit = X @ beta
residuals = y - y_fit
RSS = float(residuals.T @ residuals)
WRSS = float((sw * residuals).T @ (sw * residuals))
if len(svals) == 0 or np.min(svals) <= 0:
condition_number = float("inf")
else:
condition_number = float(np.max(svals) / np.min(svals))
warning_parts: List[str] = []
error_enabled = True
if rank < p:
error_enabled = False
warning_parts.append(
f"WARNING: Design matrix is rank deficient (rank={rank}, p={p}). "
"Least-squares coefficients were calculated, but parameter "
"errors/covariance and uncertainty bands are disabled."
)
if dof <= 0:
error_enabled = False
warning_parts.append(
f"WARNING: Residual degrees of freedom N-p = {N}-{p} = {dof}. "
"The unbiased residual variance is undefined. Parameter errors, "
"prediction errors, confidence intervals, and uncertainty bands are disabled."
)
if condition_number > cond_warning:
warning_parts.append(
f"WARNING: Ill-conditioned design matrix: cond={condition_number:.3e}. "
"Fitted coefficients may be numerically unstable."
)
if not error_enabled:
return LeastSquaresResult(
beta=beta,
beta_std=None,
cov_beta=None,
sigma2_resid=None,
residuals=residuals,
y_fit=y_fit,
N=N,
p=p,
dof=dof,
rank=rank,
RSS=RSS,
WRSS=WRSS,
singular_values=svals,
condition_number=condition_number,
error_estimation_enabled=False,
warning="\n".join(warning_parts),
known_sigma=known_sigma,
)
XtWX = Xw.T @ Xw
try:
XtWX_inv = np.linalg.inv(XtWX)
except np.linalg.LinAlgError:
# full rank のはずでも極端に悪条件の場合の保険。
XtWX_inv = np.linalg.pinv(XtWX, rcond=rcond if rcond is not None else 1e-15)
warning_parts.append(
"WARNING: np.linalg.inv failed; covariance was calculated with pseudo-inverse."
)
sigma2_resid = WRSS / dof
cov_scale = 1.0 if known_sigma else sigma2_resid
cov_beta = cov_scale * XtWX_inv
beta_var = np.diag(cov_beta)
beta_std = np.sqrt(np.maximum(beta_var, 0.0))
return LeastSquaresResult(
beta=beta,
beta_std=beta_std,
cov_beta=cov_beta,
sigma2_resid=float(sigma2_resid),
residuals=residuals,
y_fit=y_fit,
N=N,
p=p,
dof=dof,
rank=rank,
RSS=RSS,
WRSS=WRSS,
singular_values=svals,
condition_number=condition_number,
error_estimation_enabled=True,
warning="\n".join(warning_parts),
known_sigma=known_sigma,
)
[ドキュメント]
def polynomial_lsq(
x: ArrayLike,
y: ArrayLike,
order: Optional[int] = None,
*,
basis_indices: Optional[Sequence[int]] = None,
weights: Optional[ArrayLike] = None,
known_sigma: bool = False,
rcond: Optional[float] = None,
) -> LeastSquaresResult:
"""多項式線形回帰を実行するための便利な関数です。
`polynomial_design_matrix` を使用して多項式基底の設計行列を作成し、
その後 `linear_lsq` を呼び出して線形最小二乗回帰を実行します。
:param x: ArrayLike: 1次元の説明変数データ。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param order: Optional[int]: 0から `order` までのすべてのべき指数を持つ多項式基底を使用する場合の最大次数。
`basis_indices` が指定された場合は無視されます。デフォルトはNone。
:param basis_indices: Optional[Sequence[int]]: 使用する多項式のべき指数をリストで指定します。
例: [0, 1, 3]。デフォルトはNone。
:param weights: Optional[ArrayLike]: 各データ点に適用する重みベクトル。デフォルトはNone。
:param known_sigma: bool: ノイズの分散が既知として扱われるかどうか。デフォルトはFalse。
:param rcond: Optional[float]: `np.linalg.lstsq` における特異値のカットオフ値。デフォルトはNone。
:returns: LeastSquaresResult: 多項式線形最小二乗回帰の結果。
"""
X = polynomial_design_matrix(x, order=order, basis_indices=basis_indices)
return linear_lsq(X, y, weights=weights, known_sigma=known_sigma, rcond=rcond)
[ドキュメント]
def parameter_prediction_variance(X_new: np.ndarray, cov_beta: np.ndarray) -> np.ndarray:
"""新しい設計行列 `X_new` におけるパラメータ推定の不確かさに起因する分散を計算します。
`Var[y_mean] = diag(X_new @ cov_beta @ X_new.T)` の計算を行います。
これは、モデルの平均応答の分散を示します。
:param X_new: np.ndarray: 予測を行いたい新しいデータ点に対応する設計行列 (shape: N_new, p)。
:param cov_beta: np.ndarray: 推定された回帰係数の共分散行列 (shape: p, p)。
:returns: np.ndarray: 各予測点におけるパラメータ推定起因の分散のベクトル (shape: N_new,)。
"""
X_new = np.asarray(X_new, dtype=float)
cov_beta = np.asarray(cov_beta, dtype=float)
return np.sum((X_new @ cov_beta) * X_new, axis=1)
[ドキュメント]
def predict_linear(
X_new: np.ndarray,
result: LeastSquaresResult,
*,
sigma2_noise: Optional[float] = None,
include_noise: bool = True,
) -> PredictionResult:
"""線形回帰の結果から、新しいデータ点に対する予測値と誤差バンドを計算します。
指定された新しい設計行列 `X_new` に基づいて、線形モデルの平均応答と
その不確かさ（パラメータ推定に起因する誤差と、観測ノイズに起因する誤差）を計算します。
:param X_new: np.ndarray: 予測を行いたい新しいデータ点に対応する設計行列 (shape: N_new, p)。
:param result: LeastSquaresResult: `linear_lsq` 関数によって得られた線形最小二乗の結果。
:param sigma2_noise: Optional[float]: 予測に含めるノイズの分散。
Noneの場合、`result.sigma2_resid` が使用されます。デフォルトはNone。
:param include_noise: bool: 予測誤差に観測ノイズの分散を含めるかどうか。Trueの場合、`sigma2_noise` が加算されます。
デフォルトはTrue。
:returns: PredictionResult: 予測平均値、パラメータ誤差、予測誤差を格納するデータクラス。
"""
X_new = np.asarray(X_new, dtype=float)
y_mean = X_new @ result.beta
if result.cov_beta is None:
return PredictionResult(
y_mean=y_mean,
sigma_param=None,
sigma_pred=None,
var_param=None,
var_pred=None,
)
var_param = parameter_prediction_variance(X_new, result.cov_beta)
sigma_param = np.sqrt(np.maximum(var_param, 0.0))
if include_noise:
if sigma2_noise is None:
sigma2_noise = 0.0 if result.sigma2_resid is None else result.sigma2_resid
var_pred = var_param + float(sigma2_noise)
sigma_pred = np.sqrt(np.maximum(var_pred, 0.0))
else:
var_pred = None
sigma_pred = None
return PredictionResult(
y_mean=y_mean,
sigma_param=sigma_param,
sigma_pred=sigma_pred,
var_param=var_param,
var_pred=var_pred,
)
[ドキュメント]
def predict_polynomial(
x_new: ArrayLike,
result: LeastSquaresResult,
order: Optional[int] = None,
*,
basis_indices: Optional[Sequence[int]] = None,
sigma2_noise: Optional[float] = None,
include_noise: bool = True,
) -> PredictionResult:
"""多項式回帰の結果から、新しいデータ点に対する予測値と誤差バンドを計算します。
`polynomial_design_matrix` を使用して多項式基底の設計行列を作成し、
その後 `predict_linear` を呼び出して予測値を計算します。
:param x_new: ArrayLike: 予測を行いたい新しい説明変数データ点。
:param result: LeastSquaresResult: `polynomial_lsq` 関数によって得られた多項式線形最小二乗の結果。
:param order: Optional[int]: `polynomial_design_matrix` に渡す多項式の最大次数。
`basis_indices` がNoneで、かつ `order` がNoneの場合、`result.p - 1` が使用されます。
デフォルトはNone。
:param basis_indices: Optional[Sequence[int]]: `polynomial_design_matrix` に渡す多項式のべき指数リスト。デフォルトはNone。
:param sigma2_noise: Optional[float]: 予測に含めるノイズの分散。デフォルトはNone。
:param include_noise: bool: 予測誤差に観測ノイズの分散を含めるかどうか。デフォルトはTrue。
:returns: PredictionResult: 予測平均値、パラメータ誤差、予測誤差を格納するデータクラス。
"""
if basis_indices is None and order is None:
order = result.p - 1
X_new = polynomial_design_matrix(x_new, order=order, basis_indices=basis_indices)
return predict_linear(X_new, result, sigma2_noise=sigma2_noise, include_noise=include_noise)
[ドキュメント]
def measurement_std(y: ArrayLike, ddof: int = 1) -> float:
"""観測値の配列全体の標準偏差を計算し、測定ばらつきの目安として返します。
`numpy.std` を使用して標準偏差を計算します。`ddof` (Delta Degrees of Freedom) は、
計算で使う自由度の調整値であり、標本標準偏差を計算する場合は通常1とします。
:param y: ArrayLike: 観測値のベクトル。
:param ddof: int: 自由度のデルタ。標準偏差の計算でN-ddofで除算します。デフォルトは1。
:returns: float: 観測値の標準偏差。`y.size <= ddof` の場合は0.0を返します。
"""
y = np.asarray(y, dtype=float).reshape(-1)
if y.size <= ddof:
return 0.0
return float(np.std(y, ddof=ddof))
[ドキュメント]
def delta_method_variance(jacobian: np.ndarray, cov_beta: np.ndarray) -> np.ndarray:
"""デルタ法を用いて、変換された量の分散 `Var[g(beta)]` の対角要素を計算します。
デルタ法は、非線形関数 `g` の伝播誤差を近似するために使用されます。
`Var[g(beta)] = J @ cov_beta @ J.T` の対角要素を計算します。
ここで `J` は `g` の `beta` に関するヤコビ行列です。
活性化エネルギーや移動度などの派生量の、対数変換後の誤差伝播に有用です。
:param jacobian: np.ndarray: 関数 `g` の `beta` に関するヤコビ行列。
shape (N_deriv, p) または (p,) を想定します。
(p,) の場合は (1, p) に整形されます。
:param cov_beta: np.ndarray: パラメータ `beta` の共分散行列 (shape: p, p)。
:returns: np.ndarray: 各派生量に対する分散のベクトル (shape: N_deriv,)。
"""
J = np.asarray(jacobian, dtype=float)
if J.ndim == 1:
J = J.reshape(1, -1)
return np.sum((J @ cov_beta) * J, axis=1)
# =====================================================================
# Scores and information criteria
# =====================================================================
[ドキュメント]
def regression_scores(y: ArrayLike, y_fit: ArrayLike, p: int = 0) -> Dict[str, float]:
"""基本的な回帰分析のスコア（指標）を計算して返します。
残差平方和 (RSS)、平均二乗誤差 (RMSE)、決定係数 (R^2)、
自由度調整済み決定係数 (Adjusted R^2) などを計算します。
:param y: ArrayLike: 実際の観測値のベクトル。
:param y_fit: ArrayLike: モデルによる予測値または適合値のベクトル。
:param p: int: モデルのパラメータ数（説明変数の数 + 定数項の有無）。Adjusted R^2 の計算に使用されます。デフォルトは0。
:returns: Dict[str, float]: 計算された回帰スコアを含む辞書。キーは "N", "p", "RSS", "RMSE", "R2", "adj_R2"。
"""
y = np.asarray(y, dtype=float).reshape(-1)
y_fit = np.asarray(y_fit, dtype=float).reshape(-1)
residuals = y - y_fit
N = y.size
RSS = float(residuals.T @ residuals)
TSS = float(((y - np.mean(y)).T @ (y - np.mean(y))))
rmse = math.sqrt(RSS / N) if N > 0 else float("nan")
r2 = 1.0 - RSS / TSS if TSS > 0 else float("nan")
adj_r2 = 1.0 - (1.0 - r2) * (N - 1) / (N - p) if N > p and np.isfinite(r2) else float("nan")
return {"N": float(N), "p": float(p), "RSS": RSS, "RMSE": rmse, "R2": r2, "adj_R2": adj_r2}
[ドキュメント]
def normalized_weights_from_scores(scores: ArrayLike, *, lower_is_better: bool = True) -> np.ndarray:
"""与えられたモデルスコア（例: AIC, BIC, -log evidence）から、正規化された相対重み（モデル確率）を計算します。
Akaike weight や Bayesian model probability の計算に使用される、モデル間の相対的な強さを示す重みです。
スコアの差分を指数関数に適用し、合計が1になるように正規化します。
:param scores: ArrayLike: 各モデルのスコアの配列。
:param lower_is_better: bool: スコアが低いほど良いモデル（例: AIC, BIC）である場合はTrue。
スコアが高いほど良いモデル（例: log evidence）である場合はFalse。デフォルトはTrue。
:returns: np.ndarray: 各モデルの正規化された重みの配列。合計は1になります。
"""
s = np.asarray(scores, dtype=float)
if lower_is_better:
z = -0.5 * (s - np.nanmin(s))
else:
z = s - np.nanmax(s)
z = np.where(np.isfinite(z), z, -np.inf)
w = np.exp(z)
total = np.sum(w)
if total <= 0:
return np.full_like(w, 1.0 / len(w), dtype=float)
return w / total
def _normalize_candidates(
candidates: CandidateType,
labels: Optional[Sequence[str]] = None,
) -> Tuple[List[str], List[np.ndarray]]:
"""モデル候補の入力形式を正規化し、ラベルと設計行列のリストに変換します。
内部ヘルパー関数であり、`select_models_by_ic` や `select_models_by_evidence`
などのモデル選択関数で共通の入力処理を提供します。
:param candidates: CandidateType: モデル候補。
辞書 (キーはラベル、値は設計行列) または設計行列のシーケンス。
:param labels: Optional[Sequence[str]]: `candidates` がシーケンスの場合に使用するラベルのシーケンス。
Noneの場合、"model_0", "model_1" などのデフォルトラベルが生成されます。
デフォルトはNone。
:returns: Tuple[List[str], List[np.ndarray]]: モデルラベルのリストと、対応する設計行列のリストのタプル。
:raises ValueError: `labels` が指定されており、その長さが `candidates` の長さと一致しない場合。
"""
if isinstance(candidates, Mapping):
out_labels = [str(k) for k in candidates.keys()]
matrices = [np.asarray(v, dtype=float) for v in candidates.values()]
else:
matrices = [np.asarray(v, dtype=float) for v in candidates]
if labels is None:
out_labels = [f"model_{i}" for i in range(len(matrices))]
else:
if len(labels) != len(matrices):
raise ValueError("labels length and candidates length mismatch")
out_labels = [str(v) for v in labels]
return out_labels, matrices
[ドキュメント]
def select_models_by_ic(
candidates: CandidateType,
y: ArrayLike,
*,
labels: Optional[Sequence[str]] = None,
criterion: str = "BIC",
weights: Optional[ArrayLike] = None,
rcond: Optional[float] = None,
) -> ModelSelectionResult:
"""複数の候補モデルをAIC、AICc、またはBICに基づいて比較し、最良モデルを選択します。
各モデル候補に対して線形最小二乗回帰を実行し、選択された情報量基準を計算します。
その後、`normalized_weights_from_scores` を使用してモデル重み（Akaike weightsなど）を計算し、
最良のモデルを特定します。
:param candidates: CandidateType: モデル候補。設計行列のシーケンス、またはラベルをキーとする設計行列のマップ。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param labels: Optional[Sequence[str]]: `candidates` がシーケンスの場合に使用するモデルラベルのシーケンス。デフォルトはNone。
:param criterion: str: モデル選択に使用する情報量基準 ("AIC", "AICc", "BIC")。デフォルトは"BIC"。
:param weights: Optional[ArrayLike]: `linear_lsq` に渡す重みベクトル。デフォルトはNone。
:param rcond: Optional[float]: `linear_lsq` に渡す `np.linalg.lstsq` の `rcond`。デフォルトはNone。
:returns: ModelSelectionResult: モデル選択の結果、各モデルのスコア、重み、最良モデルのインデックスなど。
:raises ValueError: `criterion` が "AIC", "AICc", "BIC" のいずれでもない場合。
"""
criterion_key = criterion.upper()
if criterion_key not in {"AIC", "AICC", "BIC"}:
raise ValueError("criterion must be 'AIC', 'AICc', or 'BIC'")
if criterion_key == "AICC":
criterion_key = "AICc"
out_labels, matrices = _normalize_candidates(candidates, labels)
results = [linear_lsq(X, y, weights=weights, rcond=rcond) for X in matrices]
scores = np.array([information_criteria(r)[criterion_key] for r in results], dtype=float)
model_weights = normalized_weights_from_scores(scores, lower_is_better=True)
best_index = int(np.nanargmin(scores))
return ModelSelectionResult(
labels=out_labels,
scores=scores,
weights=model_weights,
best_index=best_index,
criterion=criterion_key,
results=results,
)
[ドキュメント]
def polynomial_candidates(x: ArrayLike, max_order: int, *, min_order: int = 0) -> Dict[str, np.ndarray]:
"""指定された次数範囲で、多項式基底に基づくモデル候補の辞書を作成します。
`min_order` から `max_order` までの各次数に対応する設計行列を生成し、
モデル選択関数で使用できる形式の辞書として返します。
:param x: ArrayLike: 説明変数データ。
:param max_order: int: 生成する多項式基底の最大次数。
:param min_order: int: 生成する多項式基底の最小次数。デフォルトは0。
:returns: Dict[str, np.ndarray]: キーがモデルラベル（例: "poly0", "poly1"）、値が設計行列の辞書。
"""
return {f"poly{d}": polynomial_design_matrix(x, order=d) for d in range(min_order, max_order + 1)}
[ドキュメント]
def basis_index_candidates(x: ArrayLike, models: Sequence[Sequence[int]]) -> Dict[str, np.ndarray]:
"""任意の多項式のべき指数リストに基づいて、モデル候補の設計行列の辞書を作成します。
特定のべき指数の組み合わせで構成される多項式モデルを柔軟に定義し、
モデル選択関数で使用できる形式の辞書として返します。
:param x: ArrayLike: 説明変数データ。
:param models: Sequence[Sequence[int]]: 各要素が1つのモデルに対応するべき指数（例: [[0, 1], [0, 1, 2]]）のシーケンス。
:returns: Dict[str, np.ndarray]: キーがモデルラベル（べき指数リストの文字列表現）、値が設計行列の辞書。
"""
return {str(list(m)): polynomial_design_matrix(x, basis_indices=m) for m in models}
# =====================================================================
# Bayesian linear regression / evidence
# =====================================================================
[ドキュメント]
def bayesian_posterior(
X: np.ndarray,
y: ArrayLike,
*,
alpha: float = 1.0,
beta: float = 1.0,
mu0: Optional[ArrayLike] = None,
prior_cov: Optional[np.ndarray] = None,
) -> BayesianLinearResult:
"""線形ガウスモデルのベイズ的な事後分布の平均と共分散を計算します。
パラメータ `beta_vec` にガウス事前分布 `N(mu0, prior_cov / alpha)` を、
ノイズにガウス分布 `N(X beta_vec, I / beta)` を仮定します。
事後分布 `p(beta_vec | y, X, alpha, beta)` の平均 `m` と共分散 `S` を算出します。
:param X: np.ndarray: 設計行列 (shape: N, p)。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param alpha: float: パラメータの事前分布の精度（逆分散）ハイパーパラメータ。デフォルトは1.0。
:param beta: float: ノイズの精度（逆分散）ハイパーパラメータ。デフォルトは1.0。
:param mu0: Optional[ArrayLike]: パラメータの事前平均ベクトル。Noneの場合、ゼロベクトルが使用されます。デフォルトはNone。
:param prior_cov: Optional[np.ndarray]: パラメータの事前共分散行列。`prior_cov / alpha` が実際の共分散となります。
Noneの場合、単位行列が使用されます。デフォルトはNone。
:returns: BayesianLinearResult: ベイズ線形回帰の事後平均、共分散、およびハイパーパラメータを含む結果。
:raises ValueError: `mu0` のサイズが `p` と一致しない場合、または `prior_cov` の形状が `(p, p)` と一致しない場合。
"""
X = np.asarray(X, dtype=float)
y = np.asarray(y, dtype=float).reshape(-1)
N, p = X.shape
if mu0 is None:
mu0_arr = np.zeros(p)
else:
mu0_arr = np.asarray(mu0, dtype=float).reshape(-1)
if mu0_arr.size != p:
raise ValueError("mu0 size mismatch")
if prior_cov is None:
prior_cov_arr = np.eye(p)
else:
prior_cov_arr = np.asarray(prior_cov, dtype=float)
if prior_cov_arr.shape != (p, p):
raise ValueError("prior_cov shape mismatch")
prior_precision = alpha * np.linalg.inv(prior_cov_arr)
S_inv = prior_precision + beta * (X.T @ X)
S = np.linalg.inv(S_inv)
m = S @ (beta * X.T @ y + prior_precision @ mu0_arr)
log_ev = bayesian_log_evidence(X, y, alpha=alpha, beta=beta, mu0=mu0_arr, prior_cov=prior_cov_arr)
return BayesianLinearResult(
mean=m,
cov=S,
alpha=float(alpha),
beta=float(beta),
sigma_noise=float(math.sqrt(1.0 / beta)),
log_evidence=log_ev,
)
[ドキュメント]
def bayesian_log_evidence(
X: np.ndarray,
y: ArrayLike,
*,
alpha: float = 1.0,
beta: float = 1.0,
mu0: Optional[ArrayLike] = None,
prior_cov: Optional[np.ndarray] = None,
) -> float:
"""線形ガウスモデルの対数エビデンス（モデル周辺尤度）を計算します。
モデル周辺尤度 `p(y | X, alpha, beta)` は、モデルがデータをどれだけよく説明するかを示す指標であり、
異なるモデルの比較に使用できます。
prior: beta_vec ~ N(mu0, prior_cov / alpha)
noise: y|beta_vec ~ N(X beta_vec, I / beta)
N が大きい場合は O(N^3) の計算コストがかかるため、多数の候補モデルや大規模データセットに
適用する場合は、高速化された実装を検討する必要があります。
:param X: np.ndarray: 設計行列 (shape: N, p)。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param alpha: float: パラメータの事前分布の精度ハイパーパラメータ。デフォルトは1.0。
:param beta: float: ノイズの精度ハイパーパラメータ。デフォルトは1.0。
:param mu0: Optional[ArrayLike]: パラメータの事前平均ベクトル。Noneの場合、ゼロベクトルが使用されます。デフォルトはNone。
:param prior_cov: Optional[np.ndarray]: パラメータの事前共分散行列。Noneの場合、単位行列が使用されます。デフォルトはNone。
:returns: float: 計算された対数エビデンス。行列式が非正の場合、`-inf` を返します。
"""
X = np.asarray(X, dtype=float)
y = np.asarray(y, dtype=float).reshape(-1)
N, p = X.shape
if mu0 is None:
mu0_arr = np.zeros(p)
else:
mu0_arr = np.asarray(mu0, dtype=float).reshape(-1)
if prior_cov is None:
prior_cov_arr = np.eye(p)
else:
prior_cov_arr = np.asarray(prior_cov, dtype=float)
actual_prior_cov = prior_cov_arr / float(alpha)
C = (1.0 / float(beta)) * np.eye(N) + X @ actual_prior_cov @ X.T
r = y - X @ mu0_arr
sign, logdet = np.linalg.slogdet(C)
if sign <= 0:
return -float("inf")
q = float(r.T @ np.linalg.solve(C, r))
return float(-0.5 * (N * math.log(2.0 * math.pi) + logdet + q))
[ドキュメント]
def empirical_bayes_linear(
X: np.ndarray,
y: ArrayLike,
*,
alpha0: float = 1.0,
beta0: Optional[float] = None,
sigma_noise0: float = 1.0,
mu0: Optional[ArrayLike] = None,
max_iter: int = 100,
tol: float = 1.0e-4,
eps: float = 1.0e-300,
) -> BayesianLinearResult:
"""MacKay型の evidence maximization を用いて、ベイズ線形回帰のハイパーパラメータ (alpha, beta) を推定します。
この関数は現在、パラメータの事前分布が等方性 `beta_vec ~ N(mu0, I/alpha)` である場合に
対応しています。期待値最大化アルゴリズムにより `alpha` と `beta` を繰り返し更新し、
最尤推定されたハイパーパラメータを求めます。
:param X: np.ndarray: 設計行列 (shape: N, p)。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param alpha0: float: `alpha` の初期値。デフォルトは1.0。
:param beta0: Optional[float]: `beta` の初期値。Noneの場合、`1.0 / (sigma_noise0 ** 2)` が使用されます。デフォルトはNone。
:param sigma_noise0: float: `beta` の初期値計算に使用されるノイズの初期標準偏差。`beta0` が指定された場合は無視されます。デフォルトは1.0。
:param mu0: Optional[ArrayLike]: パラメータの事前平均ベクトル。Noneの場合、ゼロベクトルが使用されます。デフォルトはNone。
:param max_iter: int: ハイパーパラメータ更新の最大繰り返し回数。デフォルトは100。
:param tol: float: 更新の相対変化がこの閾値より小さくなった場合に収束と見なす許容誤差。デフォルトは1.0e-4。
:param eps: float: 分母がゼロになるのを防ぐための微小な正の数。デフォルトは1.0e-300。
:returns: BayesianLinearResult: 推定されたハイパーパラメータと、それらを用いて計算された事後分布の結果。
:raises ValueError: `mu0` のサイズが `p` と一致しない場合。
"""
X = np.asarray(X, dtype=float)
y = np.asarray(y, dtype=float).reshape(-1)
N, p = X.shape
if mu0 is None:
mu0_arr = np.zeros(p)
else:
mu0_arr = np.asarray(mu0, dtype=float).reshape(-1)
if mu0_arr.size != p:
raise ValueError("mu0 size mismatch")
alpha = float(alpha0)
beta = float(beta0 if beta0 is not None else 1.0 / (sigma_noise0 ** 2))
eigvals = np.linalg.eigvalsh(X.T @ X)
converged = False
for it in range(1, max_iter + 1):
S_inv = alpha * np.eye(p) + beta * (X.T @ X)
S = np.linalg.inv(S_inv)
m = S @ (beta * X.T @ y + alpha * mu0_arr)
gamma = float(np.sum(beta * eigvals / (alpha + beta * eigvals)))
dm2 = float((m - mu0_arr).T @ (m - mu0_arr))
res = y - X @ m
rss = float(res.T @ res)
alpha_new = gamma / max(dm2, eps)
beta_new = max(N - gamma, eps) / max(rss, eps)
da = abs(alpha_new - alpha) / max(abs(alpha), eps)
db = abs(beta_new - beta) / max(abs(beta), eps)
alpha, beta = float(alpha_new), float(beta_new)
if da < tol and db < tol:
converged = True
break
posterior = bayesian_posterior(X, y, alpha=alpha, beta=beta, mu0=mu0_arr)
posterior.n_iter = it
posterior.converged = converged
return posterior
[ドキュメント]
def predict_bayesian(
X_new: np.ndarray,
result: BayesianLinearResult,
*,
include_noise: bool = True,
) -> PredictionResult:
"""ベイズ線形回帰の結果から、新しいデータ点に対する予測平均と分散を計算します。
`BayesianLinearResult` オブジェクトに含まれる事後平均と共分散、ノイズ精度を用いて、
新しい設計行列 `X_new` における予測値を計算します。
:param X_new: np.ndarray: 予測を行いたい新しいデータ点に対応する設計行列 (shape: N_new, p)。
:param result: BayesianLinearResult: `bayesian_posterior` または `empirical_bayes_linear` から得られたベイズ線形回帰の結果。
:param include_noise: bool: 予測分散に観測ノイズの分散を含めるかどうか。Trueの場合、`1/result.beta` が加算されます。
デフォルトはTrue。
:returns: PredictionResult: 予測平均値、パラメータ誤差、予測誤差を格納するデータクラス。
"""
X_new = np.asarray(X_new, dtype=float)
y_mean = X_new @ result.mean
var_param = parameter_prediction_variance(X_new, result.cov)
sigma_param = np.sqrt(np.maximum(var_param, 0.0))
if include_noise:
var_pred = var_param + 1.0 / result.beta
sigma_pred = np.sqrt(np.maximum(var_pred, 0.0))
else:
var_pred = None
sigma_pred = None
return PredictionResult(y_mean, sigma_param, sigma_pred, var_param, var_pred)
[ドキュメント]
def posterior_model_probabilities(log_evidences: ArrayLike, log_priors: Optional[ArrayLike] = None) -> np.ndarray:
"""対数エビデンスと対数事前確率から、各モデルの事後モデル確率を計算します。
ベイズモデル選択において、各モデルがデータによってどれだけ支持されているかを示す確率を計算します。
`P(M_i | D) = P(D | M_i) * P(M_i) / sum_j(P(D | M_j) * P(M_j))` の対数形式に基づきます。
:param log_evidences: ArrayLike: 各モデルの対数エビデンスの配列。
:param log_priors: Optional[ArrayLike]: 各モデルの対数事前確率の配列。
Noneの場合、すべてのモデルが等しい事前確率を持つと仮定されます。デフォルトはNone。
:returns: np.ndarray: 各モデルの事後モデル確率の配列。合計は1になります。
:raises ValueError: `log_priors` が指定されており、その形状が `log_evidences` の形状と一致しない場合。
"""
le = np.asarray(log_evidences, dtype=float)
if log_priors is None:
lp = np.zeros_like(le)
else:
lp = np.asarray(log_priors, dtype=float)
if lp.shape != le.shape:
raise ValueError("log_priors shape mismatch")
z = le + lp
z = z - np.nanmax(z)
w = np.exp(z)
total = np.sum(w)
if total <= 0:
return np.full_like(w, 1.0 / len(w), dtype=float)
return w / total
[ドキュメント]
def select_models_by_evidence(
candidates: CandidateType,
y: ArrayLike,
*,
labels: Optional[Sequence[str]] = None,
alpha: float = 1.0,
beta: float = 1.0,
log_priors: Optional[ArrayLike] = None,
) -> ModelSelectionResult:
"""複数の候補モデルをベイズ対数エビデンスに基づいて比較し、最良モデルを選択します。
各モデル候補に対して `bayesian_log_evidence` を実行し、対数エビデンスを計算します。
その後、オプションで事前モデル確率を考慮し、事後モデル確率を計算して最良のモデルを特定します。
:param candidates: CandidateType: モデル候補。設計行列のシーケンス、またはラベルをキーとする設計行列のマップ。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param labels: Optional[Sequence[str]]: `candidates` がシーケンスの場合に使用するモデルラベルのシーケンス。デフォルトはNone。
:param alpha: float: `bayesian_log_evidence` に渡すパラメータの事前分布の精度ハイパーパラメータ。デフォルトは1.0。
:param beta: float: `bayesian_log_evidence` に渡すノイズの精度ハイパーパラメータ。デフォルトは1.0。
:param log_priors: Optional[ArrayLike]: `posterior_model_probabilities` に渡す各モデルの対数事前確率。デフォルトはNone。
:returns: ModelSelectionResult: モデル選択の結果、各モデルの対数エビデンス、事後モデル確率、最良モデルのインデックスなど。
"""
out_labels, matrices = _normalize_candidates(candidates, labels)
log_evs = np.array([bayesian_log_evidence(X, y, alpha=alpha, beta=beta) for X in matrices])
probs = posterior_model_probabilities(log_evs, log_priors=log_priors)
best_index = int(np.nanargmax(probs))
return ModelSelectionResult(
labels=out_labels,
scores=log_evs,
weights=probs,
best_index=best_index,
criterion="log_evidence",
results=None,
log_evidences=log_evs,
)
# =====================================================================
# Optional sparse model selection
# =====================================================================
[ドキュメント]
def ard_select(
X: np.ndarray,
y: ArrayLike,
*,
threshold: float = 1.0e-3,
fit_intercept: bool = False,
**kwargs,
) -> Dict[str, object]:
"""sklearn の ARDRegression (Automatic Relevance Determination Regression) を使用して、
疎な基底選択（特徴量選択）を行います。
ARDRegression は、各回帰係数に異なる事前分布（精度パラメータを持つガウス分布）を仮定し、
エビデンス最大化を通じて不要な特徴量に対応する係数をゼロに近づけることで、
特徴量を自動的に選択します。
:param X: np.ndarray: 設計行列 (shape: N, p)。
:param y: ArrayLike: 観測値のベクトル (shape: N,)。
:param threshold: float: 選択される係数の絶対値の閾値。この値より小さい係数は選択されません。デフォルトは1.0e-3。
:param fit_intercept: bool: 定数項（切片）を計算に含めるかどうか。デフォルトはFalse。
:param kwargs: ARDRegression クラスに渡される追加のキーワード引数。
:returns: Dict[str, object]: 選択された係数、選択された基底のインデックス、モデルオブジェクト、
およびARDRegressionの他の属性（scores, alpha, lambda）を含む辞書。
:raises ImportError: scikit-learn ライブラリがインストールされていない場合。
"""
try:
from sklearn.linear_model import ARDRegression
except ImportError as exc:
raise ImportError("ard_select requires scikit-learn") from exc
X = np.asarray(X, dtype=float)
y = np.asarray(y, dtype=float).reshape(-1)
ard = ARDRegression(fit_intercept=fit_intercept, compute_score=True, **kwargs)
ard.fit(X, y)
coef = np.asarray(ard.coef_, dtype=float)
selected = [int(i) for i, c in enumerate(coef) if abs(c) > threshold]
return {
"model": ard,
"coef": coef,
"intercept": float(getattr(ard, "intercept_", 0.0)),
"selected": selected,
"scores": getattr(ard, "scores_", None),
"alpha": getattr(ard, "alpha_", None),
"lambda": getattr(ard, "lambda_", None),
}