tklsq.py 技術ドキュメント
このドキュメントは、Pythonライブラリ tklsq.py の技術的な構成と使用方法について説明します。
ライブラリの機能と目的
tklsq.py は、汎用的な線形回帰、誤差評価、およびモデル選択のためのユーティリティを提供するPythonライブラリです。主な目的は、線形最小二乗法のフレームワークに基づき、モデルの適合、パラメータの誤差評価、予測値の信頼区間計算、そしてAIC/BIC/AICcやベイズevidenceといった情報量基準を用いたモデル選択をサポートすることです。
このライブラリの設計方針は以下の通りです。
基本は線形最小二乗: \(y = X \beta + \epsilon\) の形式に基づきます。
パラメータ誤差: パラメータの共分散行列は \(\text{cov}(\beta) = \sigma^2 (X^T W X)^{-1}\) として計算されます。
推定値誤差: 平均応答の分散は \(\text{Var}[y_{\text{mean}}] = \text{diag}(X_{\text{new}} \text{cov}(\beta) X_{\text{new}}^T)\) です。
予測誤差: 予測値の分散は \(\text{Var}[y_{\text{pred}}] = \text{Var}[y_{\text{mean}}] + \sigma^2\) です。
モデル選択: AIC、BIC、AICc、および必要に応じてベイズevidenceを使用します。
tklsq.py は、tklib、matplotlib、openpyxl などの特定のライブラリには依存せず、計算の「コア」機能に特化しています。入出力、グラフ描画、Excelへの保存などは、このライブラリを呼び出す側のスクリプトで実装することを想定しています。
importする方法
このライブラリを他のPythonプログラムから利用するには、通常通り import ステートメントを使用します。
import tklsq
あるいは、特定の関数やクラスのみをインポートすることも可能です。
from tklsq import linear_lsq, LeastSquaresResult
必要な非標準ライブラリとインストール方法
tklsq.py のコア機能は、主に数値計算ライブラリ numpy に依存しています。numpy は、以下のコマンドでインストールできます。
pip install numpy
また、ard_select 関数は、機械学習ライブラリ scikit-learn に依存しています。この関数を使用する場合のみ、scikit-learn のインストールが必要です。
pip install scikit-learn
importできる変数と関数
tklsq.py は、データ型エイリアス、データクラス、および多数の関数を提供します。
型エイリアス
ArrayLike:Union[Sequence[float], np.ndarray]NumPy配列またはfloatのシーケンス(リスト、タプルなど)のいずれかを許容する型ヒント。
CandidateType:Union[Sequence[np.ndarray], Mapping[str, np.ndarray]]NumPy配列のシーケンス、または文字列をキーとするNumPy配列のマッピング(辞書など)のいずれかを許容する型ヒント。
データクラス
線形最小二乗法、予測、モデル選択、およびベイズ線形回帰の結果を構造化して保持するためのデータクラスです。
LeastSquaresResult
線形最小二乗の結果を格納します。
beta:np.ndarray推定された回帰係数ベクトル。
beta_std:Optional[np.ndarray]各回帰係数の標準誤差。誤差推定が無効な場合は
None。
cov_beta:Optional[np.ndarray]回帰係数の共分散行列。誤差推定が無効な場合は
None。
sigma2_resid:Optional[float]残差の分散 (\(\sigma^2\)) の推定値。誤差推定が無効な場合は
None。
residuals:np.ndarray残差 (\(y - \hat{y}\))。
y_fit:np.ndarray適合値 (\(\hat{y}\))。
N:intデータ点の数。
p:intパラメータの数 (設計行列の列数)。
dof:int残差の自由度 (\(N-p\))。
rank:int設計行列のランク。
RSS:float残差平方和 (Residual Sum of Squares)。
WRSS:float重み付き残差平方和 (Weighted Residual Sum of Squares)。
singular_values:np.ndarray設計行列の特異値。
condition_number:float設計行列の条件数。
error_estimation_enabled:bool誤差推定が成功したかどうか。
warning:str計算中に発生した警告メッセージ。デフォルトは空文字列。
known_sigma:boolノイズ分散が既知として扱われたかどうか。デフォルトは
False。
sigma_resid:Optional[float](プロパティ)残差の標準偏差 (\(\sigma\)) の推定値。
sigma2_residがNoneの場合はNone。
PredictionResult
線形モデルの推定値と誤差バンドを格納します。
y_mean:np.ndarray予測された平均応答値。
sigma_param:Optional[np.ndarray]パラメータ推定の不確実性に起因する標準偏差 (信頼区間)。
var_paramがNoneの場合はNone。
sigma_pred:Optional[np.ndarray]予測値の標準偏差 (予測区間)。
var_predがNoneの場合はNone。
var_param:Optional[np.ndarray]パラメータ推定の不確実性に起因する分散。
var_pred:Optional[np.ndarray]予測値の分散。
ModelSelectionResult
複数モデル比較の結果を格納します。
labels:List[str]各モデルのラベル。
scores:np.ndarray各モデルのスコア (例: AIC, BIC, log evidence)。
weights:np.ndarray各モデルの相対的な重み (例: AIC/BIC重み、事後モデル確率)。
best_index:int最もスコアの良いモデルのインデックス。
criterion:strモデル選択に使用された基準の名称。
results:Optional[List[LeastSquaresResult]]各モデルの
LeastSquaresResultオブジェクトのリスト。
log_evidences:Optional[np.ndarray]各モデルのベイズ対数evidence。
best_label:str(プロパティ)最もスコアの良いモデルのラベル。
BayesianLinearResult
ベイズ線形回帰の事後分布パラメータを格納します。
mean:np.ndarray事後平均回帰係数ベクトル。
cov:np.ndarray事後共分散行列。
alpha:float事前分布の精度パラメータ (回帰係数の分散の逆数に比例)。
beta:floatノイズの精度パラメータ (ノイズ分散の逆数)。
sigma_noise:float推定されたノイズの標準偏差 (\(\sqrt{1/\beta}\))。
log_evidence:Optional[float]モデルの対数evidence。
n_iter:int経験的ベイズ法でのイテレーション数。デフォルトは0。
converged:bool経験的ベイズ法が収束したかどうか。デフォルトは
False。
関数
Design matrices
設計行列を生成するためのユーティリティ関数。
polynomial_design_matrix(x, order=None, basis_indices=None)
多項式基底の設計行列を作成します。
動作:
xの各要素を \(x^0, x^1, \dots, x^{\text{order}}\) または指定されたbasis_indicesのべき指数で変換し、設計行列を生成します。引数:
x:ArrayLike1次元の説明変数。
order:Optional[int]\(0, 1, \dots, \text{order}\) の全てのべき指数を基底として使用する場合に指定します。
basis_indices:Optional[Sequence[int]]使用するべき指数を任意のシーケンスで指定します (例:
[0, 1, 3])。orderまたはbasis_indicesのいずれかを指定する必要があります。
戻り値:
np.ndarray, shape(N, p)生成された設計行列。
design_matrix_from_functions(x, basis_functions)
任意の基底関数リストから設計行列を作成します。
動作: 提供された基底関数を
xの各要素に適用し、設計行列の各列を生成します。引数:
x:ArrayLike1次元の説明変数。
basis_functions:Sequence各要素が
f(x_array)を返す関数、またはnp.vectorize可能なf(x_scalar)関数であるシーケンス。
戻り値:
np.ndarray, shape(N, p)生成された設計行列。
add_intercept(X)
設計行列の先頭列に定数項 (切片) を追加します。
動作: 既存の設計行列
Xの左端に、全ての要素が1である列を追加します。引数:
X:np.ndarray元の設計行列。
戻り値:
np.ndarray定数項が追加された設計行列。
Least squares and uncertainty
線形最小二乗計算と不確実性評価のための関数。
linear_lsq(X, y, *, weights=None, known_sigma=False, rcond=None, cond_warning=1.0e12)
線形最小二乗を実行し、可能であれば共分散行列と誤差を計算します。
動作: 設計行列
Xと観測値yを用いて、重み付き線形最小二乗法により回帰係数 \(\beta\) を推定します。ランク落ちや自由度が不足している場合でも \(\beta\) は返しますが、誤差推定は無効になります。引数:
X:np.ndarray, shape(N, p)設計行列。
y:ArrayLike, shape(N,)観測値。
weights:Optional[ArrayLike]重み。通常は \(1/\sigma_y^2\) を指定します。
Noneの場合は等重みとして扱われます。
known_sigma:boolTrueの場合、weightsが絶対的な \(1/\sigma_y^2\) であるとみなし、\(\text{cov}(\beta) = (X^T W X)^{-1}\) を使用します。Falseの場合は、残差から \(\sigma^2\) を推定して共分散行列に乗じます。
rcond:Optional[float]np.linalg.lstsqにおける特異値のカットオフ値を指定します。
cond_warning:float設計行列の条件数がこの値を超えると、結果の
warningフィールドに警告メッセージが追記されます。
戻り値:
LeastSquaresResult線形最小二乗の結果を格納したデータクラス。
polynomial_lsq(x, y, order=None, *, basis_indices=None, weights=None, known_sigma=False, rcond=None)
多項式線形回帰の簡易関数。
動作:
polynomial_design_matrixを使用して設計行列を生成し、linear_lsqを呼び出して多項式回帰を実行します。引数:
x:ArrayLike1次元の説明変数。
y:ArrayLike観測値。
order:Optional[int]多項式の最大次数。
basis_indices:Optional[Sequence[int]]使用するべき指数。
weights:Optional[ArrayLike]linear_lsqに渡す重み。
known_sigma:boollinear_lsqに渡すknown_sigma。
rcond:Optional[float]linear_lsqに渡すrcond。
戻り値:
LeastSquaresResult多項式線形回帰の結果を格納したデータクラス。
parameter_prediction_variance(X_new, cov_beta)
\(\text{Var}[y_{\text{mean}}] = \text{diag}(X_{\text{new}} \text{cov}(\beta) X_{\text{new}}^T)\) を計算します。
動作: 新しい説明変数
X_newと回帰係数の共分散行列cov_betaを用いて、予測される平均応答値の分散を計算します。これは信頼区間の計算に使用されます。引数:
X_new:np.ndarray新しいデータ点の設計行列。
cov_beta:np.ndarray回帰係数の共分散行列。
戻り値:
np.ndarray各データ点における予測平均の分散値の配列。
predict_linear(X_new, result, *, sigma2_noise=None, include_noise=True)
線形回帰結果から、推定値と誤差バンドを計算します。
動作:
linear_lsqで得られたLeastSquaresResultを使用して、新しいデータ点X_newにおける平均応答値、信頼区間、および予測区間を計算します。引数:
X_new:np.ndarray予測を行いたい新しいデータ点の設計行列。
result:LeastSquaresResultlinear_lsqの結果オブジェクト。
sigma2_noise:Optional[float]ノイズの分散を明示的に指定する場合に用います。
Noneの場合、result.sigma2_residが使用されます。
include_noise:boolTrueの場合、予測誤差にノイズの分散を含めて予測区間を計算します。Falseの場合、パラメータ誤差のみに基づいた信頼区間のみを計算します。
戻り値:
PredictionResult予測結果を格納したデータクラス。
predict_polynomial(x_new, result, order=None, *, basis_indices=None, sigma2_noise=None, include_noise=True)
多項式回帰結果から、推定値と誤差バンドを計算します。
動作:
polynomial_design_matrixを使用してx_newから設計行列を生成し、predict_linearを呼び出して多項式回帰の予測値を計算します。引数:
x_new:ArrayLike予測を行いたい新しい説明変数。
result:LeastSquaresResultpolynomial_lsqの結果オブジェクト。
order:Optional[int]多項式の最大次数。
basis_indicesもorderもNoneの場合、result.p - 1が使用されます。
basis_indices:Optional[Sequence[int]]使用するべき指数。
sigma2_noise:Optional[float]predict_linearに渡すノイズ分散。
include_noise:boolpredict_linearに渡すinclude_noise。
戻り値:
PredictionResult多項式回帰の予測結果を格納したデータクラス。
measurement_std(y, ddof=1)
観測値全体の標準偏差を測定ばらつきの目安として返します。
動作: 与えられたデータ
yの標本標準偏差を計算します。引数:
y:ArrayLike観測値。
ddof:int自由度 (degrees of freedom) の調整。デフォルトは1 (標本標準偏差)。
戻り値:
float計算された標準偏差。
delta_method_variance(jacobian, cov_beta)
デルタ法を用いて \(\text{Var}[g(\beta)] = J \text{cov}(\beta) J^T\) の対角要素を計算します。
動作: パラメータ \(\beta\) の関数 \(g(\beta)\) の分散を、ヤコビアン \(J\) と \(\beta\) の共分散行列 \(\text{cov}(\beta)\) を用いて近似します。これは、変換された量の誤差伝播に有用です。
引数:
jacobian:np.ndarray関数 \(g\) のヤコビアン行列。shape は
(N, p)または(p,)。
cov_beta:np.ndarray回帰係数の共分散行列。
戻り値:
np.ndarray近似された分散値の配列。
Scores and information criteria
回帰モデルのスコアと情報量基準を計算するための関数。
regression_scores(y, y_fit, p=0)
基本的な回帰スコアを返します。
動作: 観測値
yと適合値y_fitを比較し、残差平方和 (RSS)、平均二乗誤差 (RMSE)、決定係数 (\(R^2\))、調整済み決定係数 (Adjusted \(R^2\)) などを計算します。引数:
y:ArrayLike観測値。
y_fit:ArrayLikeモデルによる適合値。
p:intモデルのパラメータ数。Adjusted \(R^2\) の計算に使用されます。
戻り値:
Dict[str, float]計算されたスコアの辞書。キーは "N", "p", "RSS", "RMSE", "R2", "adj_R2"。
gaussian_log_likelihood_from_rss(RSS, N, sigma2=None)
ガウス残差を仮定した対数尤度を計算します。
動作: 残差平方和 (RSS) とデータ点数 (N) に基づいて、ガウス分布に従う残差の対数尤度を計算します。
引数:
RSS:float残差平方和。
N:intデータ点の数。
sigma2:Optional[float]既知のノイズ分散。
Noneの場合、MLE (\(\text{RSS}/\text{N}\)) を用います。
戻り値:
float計算された対数尤度。
information_criteria(result, *, k=None)
AIC, AICc, BIC を計算します。
動作:
LeastSquaresResultオブジェクトから得られる情報(データ点数 N、パラメータ数 p、RSS)を用いて、赤池情報量基準 (AIC)、修正赤池情報量基準 (AICc)、ベイズ情報量基準 (BIC) を計算します。引数:
result:LeastSquaresResult線形最小二乗の結果オブジェクト。
k:Optional[int]情報量基準の計算に使用するパラメータの数。
Noneの場合、result.pが使用されます。
戻り値:
Dict[str, float]計算された情報量基準の辞書。キーは "logL", "AIC", "AICc", "BIC"。
normalized_weights_from_scores(scores, *, lower_is_better=True)
AIC/BIC や -log evidence から相対重みを計算します。
動作: 一連のモデルスコアから、各モデルの相対的な重み (例: AIC/BIC重み) を計算します。これは、スコアの最小値または最大値を基準に指数関数的に変換することで行われます。
引数:
scores:ArrayLike各モデルのスコアの配列。
lower_is_better:boolTrueの場合、スコアが小さいほど良いモデルとみなします (AIC, BICなど)。Falseの場合、スコアが大きいほど良いモデルとみなします (log evidenceなど)。
戻り値:
np.ndarray正規化された相対重みの配列。
select_models_by_ic(candidates, y, *, labels=None, criterion='BIC', weights=None, rcond=None)
候補設計行列を AIC/AICc/BIC で比較します。
動作: 複数の候補設計行列をそれぞれ線形最小二乗法で適合させ、指定された情報量基準 (
AIC,AICc, またはBIC) に基づいてモデルを比較し、最も良いモデルを選択します。引数:
candidates:CandidateType比較する候補設計行列のシーケンスまたは辞書。
y:ArrayLike観測値。
labels:Optional[Sequence[str]]候補モデルに割り当てるラベルのシーケンス。
criterion:strモデル選択に使用する情報量基準 ("AIC", "AICC", または "BIC")。デフォルトは "BIC"。
weights:Optional[ArrayLike]linear_lsqに渡す重み。
rcond:Optional[float]linear_lsqに渡すrcond。
戻り値:
ModelSelectionResultモデル選択の結果を格納したデータクラス。
polynomial_candidates(x, max_order, *, min_order=0)
\(0 \dots \text{max_order}\) の多項式候補を辞書として作成します。
動作: 指定された
min_orderからmax_orderまでの各次数について、polynomial_design_matrixを使用して設計行列を生成し、辞書形式で返します。引数:
x:ArrayLike説明変数。
max_order:int含める多項式の最大次数。
min_order:int含める多項式の最小次数。デフォルトは0。
戻り値:
Dict[str, np.ndarray]各次数の多項式設計行列を格納した辞書。キーは "polyN" の形式。
basis_index_candidates(x, models)
任意の多項式べき指数リストから候補設計行列を作成します。
動作: 指定されたべき指数リストのシーケンスに基づき、
polynomial_design_matrixを使用して設計行列を生成し、辞書形式で返します。引数:
x:ArrayLike説明変数。
models:Sequence[Sequence[int]]各要素が使用するべき指数を指定するシーケンスのシーケンス (例:
[[0, 1], [0, 2]])。
戻り値:
Dict[str, np.ndarray]各べき指数モデルの設計行列を格納した辞書。キーはべき指数リストの文字列表現。
Bayesian linear regression / evidence
ベイズ線形回帰とモデルevidenceの計算のための関数。
bayesian_posterior(X, y, *, alpha=1.0, beta=1.0, mu0=None, prior_cov=None)
線形ガウスモデルの事後平均と共分散を計算します。
動作: 事前分布がガウス分布、ノイズがガウス分布と仮定された線形モデルの事後分布の平均と共分散を解析的に計算します。
引数:
X:np.ndarray設計行列。
y:ArrayLike観測値。
alpha:float事前分布の精度パラメータ (係数分散の逆数に比例)。デフォルトは1.0。
beta:floatノイズの精度パラメータ (ノイズ分散の逆数)。デフォルトは1.0。
mu0:Optional[ArrayLike]事前平均回帰係数ベクトル。
Noneの場合、ゼロベクトルが使用されます。
prior_cov:Optional[np.ndarray]事前共分散行列。
Noneの場合、単位行列が使用されます。
戻り値:
BayesianLinearResultベイズ線形回帰の事後分布パラメータを格納したデータクラス。
bayesian_log_evidence(X, y, *, alpha=1.0, beta=1.0, mu0=None, prior_cov=None)
線形ガウスモデルの対数evidenceを計算します。
動作: モデルの対数evidence (周辺尤度) を計算します。これはモデル選択において、異なるモデルを比較するための指標として使用されます。
引数:
X:np.ndarray設計行列。
y:ArrayLike観測値。
alpha:float事前分布の精度パラメータ。デフォルトは1.0。
beta:floatノイズの精度パラメータ。デフォルトは1.0。
mu0:Optional[ArrayLike]事前平均回帰係数ベクトル。
prior_cov:Optional[np.ndarray]事前共分散行列。
戻り値:
float計算された対数evidence。
empirical_bayes_linear(X, y, *, alpha0=1.0, beta0=None, sigma_noise0=1.0, mu0=None, max_iter=100, tol=1.0e-4, eps=1.0e-300)
MacKay型のevidence maximizationで alpha と beta を推定します。
動作: 経験的ベイズ (Type-II Maximum Likelihood) を使用して、モデルの超パラメータ
alpha(事前分布の精度) とbeta(ノイズの精度) をデータから推定します。これにより、回帰係数のスパース性やノイズレベルを自動的に決定できます。現在の実装は等方性事前分布 (beta_vec ~ N(mu0, I/alpha)) 用です。引数:
X:np.ndarray設計行列。
y:ArrayLike観測値。
alpha0:floatalphaの初期値。デフォルトは1.0。
beta0:Optional[float]betaの初期値。Noneの場合、1.0 / (sigma_noise0 ** 2)が使用されます。
sigma_noise0:floatノイズの標準偏差の初期値。デフォルトは1.0。
mu0:Optional[ArrayLike]事前平均回帰係数ベクトル。
Noneの場合、ゼロベクトルが使用されます。
max_iter:int最大イテレーション数。デフォルトは100。
tol:float収束判定のための許容誤差。デフォルトは1.0e-4。
eps:floatゼロ除算を避けるための小さな値。デフォルトは1.0e-300。
戻り値:
BayesianLinearResult推定された超パラメータと事後分布を格納したデータクラス。
predict_bayesian(X_new, result, *, include_noise=True)
ベイズ線形回帰結果から予測平均・分散を計算します。
動作:
BayesianLinearResultを使用して、新しいデータ点X_newにおける予測平均と、信頼区間 (パラメータの不確実性) および予測区間 (パラメータとノイズの不確実性) を計算します。引数:
X_new:np.ndarray予測を行いたい新しいデータ点の設計行列。
result:BayesianLinearResultbayesian_posteriorまたはempirical_bayes_linearの結果オブジェクト。
include_noise:boolTrueの場合、予測誤差にノイズの分散を含めて予測区間を計算します。Falseの場合、パラメータ誤差のみに基づいた信頼区間のみを計算します。
戻り値:
PredictionResult予測結果を格納したデータクラス。
posterior_model_probabilities(log_evidences, log_priors=None)
対数evidenceと対数事前確率から事後モデル確率を計算します。
動作: 複数のモデルの対数evidenceと、それぞれのモデルの対数事前確率 (もしあれば) を組み合わせて、各モデルの事後確率を計算します。
引数:
log_evidences:ArrayLike各モデルの対数evidenceの配列。
log_priors:Optional[ArrayLike]各モデルの対数事前確率の配列。
Noneの場合、全てのモデルに等しい事前確率が仮定されます。
戻り値:
np.ndarray正規化された事後モデル確率の配列。
select_models_by_evidence(candidates, y, *, labels=None, alpha=1.0, beta=1.0, log_priors=None)
候補設計行列を Bayesian log evidence で比較します。
動作: 複数の候補設計行列をそれぞれベイズ線形回帰で適合させ、対数evidenceに基づいてモデルを比較し、事後モデル確率から最も良いモデルを選択します。
引数:
candidates:CandidateType比較する候補設計行列のシーケンスまたは辞書。
y:ArrayLike観測値。
labels:Optional[Sequence[str]]候補モデルに割り当てるラベルのシーケンス。
alpha:floatbayesian_log_evidenceに渡す事前分布の精度パラメータ。デフォルトは1.0。
beta:floatbayesian_log_evidenceに渡すノイズの精度パラメータ。デフォルトは1.0。
log_priors:Optional[ArrayLike]posterior_model_probabilitiesに渡す対数事前確率。
戻り値:
ModelSelectionResultモデル選択の結果を格納したデータクラス。
Optional sparse model selection
疎な基底選択のための関数。scikit-learn に依存します。
ard_select(X, y, *, threshold=1.0e-3, fit_intercept=False, **kwargs)
sklearn の ARDRegression を使った疎な基底選択。
動作: Automatic Relevance Determination (ARD) 回帰を用いて、与えられた基底のうち重要性の低いものを自動的にゼロに近い値にすることで、疎なモデルを選択します。
scikit-learnがインストールされている必要があります。引数:
X:np.ndarray設計行列。
y:ArrayLike観測値。
threshold:float係数の絶対値がこの値よりも小さい基底は選択されないと判断されます。デフォルトは1.0e-3。
fit_intercept:bool切片項をモデルに含めるかどうか。デフォルトは
False。
**kwargs:sklearn.linear_model.ARDRegressionに渡される追加のキーワード引数。
戻り値:
Dict[str, object]ARD回帰モデルオブジェクト、推定された係数、切片、選択された基底のインデックス、スコア、
alpha、lambdaを含む辞書。
main scriptとして実行したときの動作
tklsq.py には、if __name__ == "__main__": ブロックが定義されていません。
したがって、このファイルをPythonインタープリタで直接実行 (python tklsq.py) しても、特定のスクリプトとしての動作は実行されません。このライブラリは、他のプログラムからインポートされて利用されることを想定した、機能を提供するモジュールです。