すべてのカテゴリ

Pythonで科学計算やデータ解析をしていると、χ²計算、尤度計算、MCMC、回帰分析などで共分散行列を扱う場面があります。このとき、よくある実装としてC^{-1}を明示的に作り、r.T @ inv(C) @ rのように計算するコードがあります。数式としては正しく見えますが、数値計算の実装としてはあまり推奨されません。理由は、逆行列を明示的に作ることで、計算コストや数値誤差が増えやすくなるためです。より安定した方法は、共分散行列 C を Cholesky 分解し、C = L L.Tとして、残差ベクトル r に対してL w = rを解く方法です。その上で、χ² = w.T @ wとして計算します。この方法は、明示的に C^{-1} を作らずに、同じ共分散重み付きχ²を評価できます。Pythonでは、たとえば以下のように書けます。mport numpy as npfrom scipy.linalg import solve_triangulardef eval_whitened_chi2(L, y_obs, y_pred):    r = np.asarray(y_obs) - np.asarray(y_pred)    w = solve_triangular(L, r, lower=True, check_finite=False)    return float(np.dot(w, w))ここで L は、あらかじめL = np.linalg.cholesky(C)として計算しておきます。この方法の利点は以下です。・明示的な逆行列を作らない  ・数値的に安定しやすい  ・同じ共