hanalyze-design-0.2.0.1: src/Hanalyze/Design/Diagnostics.hs
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE BangPatterns #-}
-- |
-- Module : Hanalyze.Design.Diagnostics
-- Description : DoE 設計診断 (Alias Matrix / VIF / D-A-G-I efficiency の一括算出)
-- Copyright : (c) 2026 Aelysce Project (Toshiaki Honda)
-- License : BSD-3-Clause
--
-- [日本語]: DoE 設計診断: Alias Matrix / VIF / D-A-G-I efficiency。
--
-- 設計行列 @X@ (n × p) に対して、 multicollinearity と最適性指標を一括算出。
--
-- === efficiency 指標
--
-- - D-efficiency = (|XᵀX| / n^p)^{1/p}
-- - A-efficiency = p / trace((XᵀX/n)⁻¹)
-- - G-efficiency = p / max_i (n · x_iᵀ (XᵀX)⁻¹ x_i)
-- - I-efficiency = 1 / (n · trace((XᵀX)⁻¹ · M))、
-- M = 1/n · XᵀX (= self-moment 近似版)
--
-- === VIF
--
-- 各列 j について、 「j 以外の列で j を回帰」 した R² を用いて
-- @VIF_j = 1 / (1 − R²_j)@。 切片を含む X を想定し、 切片列 (= 全 1) は
-- スキップ。
--
-- === Alias Matrix
--
-- @A = (XᵀX)⁻¹ Xᵀ Z@、 ここで Z は 「設計に入れていない交絡項」 のモデル
-- 行列。 ここでは Z を Optional 引数として取り、 未指定の場合は
-- アライアス対象が無いとして空行列を返す。
--
-- [English]: DoE design diagnostics: Alias Matrix \/ VIF \/ D-A-G-I
-- efficiency.
--
-- Computes multicollinearity and optimality metrics in bulk for the design
-- matrix @X@ (n × p).
--
-- === Efficiency metrics
--
-- - D-efficiency = (|XᵀX| / n^p)^{1/p}
-- - A-efficiency = p / trace((XᵀX/n)⁻¹)
-- - G-efficiency = p / max_i (n · x_iᵀ (XᵀX)⁻¹ x_i)
-- - I-efficiency = 1 / (n · trace((XᵀX)⁻¹ · M)),
-- M = 1/n · XᵀX (a self-moment approximation)
--
-- === VIF
--
-- For each column j, using the R² of "regressing j on all columns other
-- than j", @VIF_j = 1 / (1 − R²_j)@. Assumes X includes an intercept, and
-- the intercept column (all 1s) is skipped.
--
-- === Alias Matrix
--
-- @A = (XᵀX)⁻¹ Xᵀ Z@, where Z is the model matrix of "confounding terms
-- not included in the design". Z is taken here as an optional argument;
-- when unspecified, an empty matrix is returned as there is nothing to
-- alias.
module Hanalyze.Design.Diagnostics
( DesignDiagnostics (..)
, diagnostics
, diagnosticsWithAlias
, vifVector
, aliasMatrix
) where
import qualified Numeric.LinearAlgebra as LA
-- ===========================================================================
-- 型
-- ===========================================================================
data DesignDiagnostics = DesignDiagnostics
{ ddVIF :: !(LA.Vector Double)
, ddDEff :: !Double
, ddAEff :: !Double
, ddGEff :: !Double
, ddIEff :: !Double
, ddAliasMatrix :: !(LA.Matrix Double)
} deriving (Show)
-- ===========================================================================
-- 公開 API
-- ===========================================================================
-- | [日本語]: Alias を含めない簡易版 (Z = 空)。
-- [English]: A simplified version without Alias (Z = empty).
diagnostics :: LA.Matrix Double -> DesignDiagnostics
diagnostics x =
let dd = computeDiagnostics x
in dd { ddAliasMatrix = LA.fromLists [[]] }
-- | [日本語]: Z (交絡対象モデル行列) 込みの完全版。 Z の行数は X と一致する必要がある。
-- [English]: The full version, including Z (the model matrix of
-- confounding terms). Z's row count must match X's.
diagnosticsWithAlias :: LA.Matrix Double -> LA.Matrix Double -> DesignDiagnostics
diagnosticsWithAlias x z =
let dd = computeDiagnostics x
a = aliasMatrix x z
in dd { ddAliasMatrix = a }
-- | [日本語]: VIF を各列について返す (全 1 列は VIF = 1)。
-- [English]: Returns the VIF for each column (an all-1s column has
-- VIF = 1).
vifVector :: LA.Matrix Double -> LA.Vector Double
vifVector x =
let p = LA.cols x
in LA.fromList [ vifForCol x j | j <- [0 .. p - 1] ]
-- | [日本語]: Alias matrix A = (XᵀX)⁻¹ Xᵀ Z。
-- [English]: The alias matrix A = (XᵀX)⁻¹ Xᵀ Z.
aliasMatrix :: LA.Matrix Double -> LA.Matrix Double -> LA.Matrix Double
aliasMatrix x z =
let xtx = LA.tr x LA.<> x
d = LA.det xtx
in if abs d < 1e-12
then LA.fromLists [[]]
else LA.inv xtx LA.<> LA.tr x LA.<> z
-- ===========================================================================
-- 内部
-- ===========================================================================
computeDiagnostics :: LA.Matrix Double -> DesignDiagnostics
computeDiagnostics x =
let n = LA.rows x
p = LA.cols x
nD = fromIntegral n :: Double
pD = fromIntegral p :: Double
xtx = LA.tr x LA.<> x
d = LA.det xtx
singular = abs d < 1e-12
inv = if singular then LA.ident p else LA.inv xtx
-- D-efficiency: (|XᵀX| / n^p)^{1/p} (clamped to ≥ 0)
dEff = if singular || d <= 0 then 0
else (d / (nD ** pD)) ** (1 / pD)
-- A-efficiency: p / trace((XᵀX/n)⁻¹) = p · n / trace((XᵀX)⁻¹)... wait
-- (XᵀX / n)⁻¹ = n · (XᵀX)⁻¹ なので trace((XᵀX/n)⁻¹) = n · trace((XᵀX)⁻¹)
-- → A-eff = p / (n · trace((XᵀX)⁻¹))
trInv = sum [ inv `LA.atIndex` (i, i) | i <- [0 .. p - 1] ]
aEff = if singular || trInv == 0 then 0
else pD / (nD * trInv)
-- G-efficiency: p / max_i (n · h_ii)、 h_ii = x_iᵀ (XᵀX)⁻¹ x_i
hMax = if singular then 1
else maximum
[ let xi = LA.flatten (x LA.? [i])
v = inv LA.#> xi
in xi `LA.dot` v
| i <- [0 .. n - 1] ]
gEff = if hMax == 0 then 0 else pD / (nD * hMax)
-- I-efficiency 近似 (self-moment)
iEff = if singular then 0
else
let m = LA.scale (1 / nD) xtx
t = LA.sumElements (LA.takeDiag (inv LA.<> m))
in if t == 0 then 0 else 1 / (nD * t)
in DesignDiagnostics
{ ddVIF = vifVector x
, ddDEff = dEff
, ddAEff = aEff
, ddGEff = gEff
, ddIEff = iEff
, ddAliasMatrix = LA.fromLists [[]]
}
-- | [日本語]: 列 j の VIF。 全 1 列 (切片) は 1 を返す。
-- [English]: The VIF of column j. Returns 1 for an all-1s column
-- (intercept).
vifForCol :: LA.Matrix Double -> Int -> Double
vifForCol x j =
let col = LA.flatten (x LA.¿ [j])
isConst = let c0 = LA.atIndex col 0
in LA.maxElement (LA.cmap (\v -> abs (v - c0)) col) < 1e-12
in if isConst then 1
else
let p = LA.cols x
others = [ k | k <- [0 .. p - 1], k /= j ]
xOthers = x LA.¿ others
yj = col
xtx = LA.tr xOthers LA.<> xOthers
d = LA.det xtx
in if abs d < 1e-12 then 1 / 0
else
let beta = LA.inv xtx LA.#> (LA.tr xOthers LA.#> yj)
yhat = xOthers LA.#> beta
yBar = LA.sumElements yj / fromIntegral (LA.size yj)
ssR = LA.sumElements ((yj - yhat) ^ (2 :: Int))
ssT = LA.sumElements ((yj - LA.scalar yBar) ^ (2 :: Int))
r2 = if ssT == 0 then 0 else 1 - ssR / ssT
in if r2 >= 1 then 1 / 0 else 1 / (1 - r2)