hanalyze-0.2.0.0: src/Hanalyze/Design/RSM.hs
{-# LANGUAGE OverloadedStrings #-}
-- |
-- Module : Hanalyze.Design.RSM
-- Description : 応答曲面法 (RSM) — CCD/Box-Behnken 計画・二次モデル fit・極値の解析解
-- Copyright : (c) 2026 Aelysce Project (Toshiaki Honda)
-- License : BSD-3-Clause
--
-- Response Surface Methodology (RSM).
--
-- * 'centralComposite' — central composite design (CCD): @2^k@ factorial
-- + axial points + center points.
-- * 'boxBehnken' — Box-Behnken design: @k@ three-level factors
-- without axial points.
-- * 'quadraticDesign' — design matrix for the quadratic model
-- (intercept + main + squared + interaction terms).
-- * 'fitQuadratic' — fit the quadratic regression by least squares.
-- * 'optimumPoint' — analytically solve for the extremum (max / min)
-- from the fit.
module Hanalyze.Design.RSM
( CCDType (..)
, centralComposite
, centralCompositeRotatable
, boxBehnken
, quadraticDesign
, quadraticTermNames
, QuadFit (..)
, fitQuadratic
, optimumPoint
, quadBMatrix
, canonicalAnalysis
) where
import Data.List (sortOn)
import Data.Text (Text)
import qualified Data.Text as T
import qualified Numeric.LinearAlgebra as LA
import Hanalyze.Design.Factorial (twoLevelFactorial)
-- ---------------------------------------------------------------------------
-- 中心複合計画 (CCD)
-- ---------------------------------------------------------------------------
-- | Central composite design (CCD) type.
data CCDType
= CCC Double -- ^ Circumscribed: axial distance @α@ (the rotatable
-- choice is @(2^k)^{1/4}@).
| CCF -- ^ Face-centered: @α = 1@ (axial points sit on the cube faces).
| CCI Double -- ^ Inscribed: @α = 1@, factorial part scaled by @1/α@.
deriving (Show, Eq)
-- | Central composite design.
--
-- Composition:
--
-- * @2^k@ factorial part: every @±1@ combination (@2^k@ rows).
-- * @2k@ axial points: @(±α, 0, …, 0)@ for each factor.
-- * @nC@ center points at @(0, …, 0)@.
--
-- @centralComposite k ccdType nC@: @k@ factors and @nC@ centre points.
centralComposite :: Int -> CCDType -> Int -> [[Double]]
centralComposite k ccdType nC =
let factorial = case ccdType of
CCI alpha ->
[[v / alpha | v <- row] | row <- twoLevelFactorial k]
_ -> twoLevelFactorial k
alpha = case ccdType of
CCC a -> a
CCF -> 1.0
CCI _ -> 1.0
axial = concat
[ [ replicate i 0 ++ [-alpha] ++ replicate (k - 1 - i) 0
, replicate i 0 ++ [ alpha] ++ replicate (k - 1 - i) 0
]
| i <- [0 .. k - 1] ]
center = replicate nC (replicate k 0)
in factorial ++ axial ++ center
-- | Rotatable CCD with @α = (2^k)^{1/4}@.
centralCompositeRotatable :: Int -> Int -> [[Double]]
centralCompositeRotatable k nC =
let alpha = (fromIntegral (2 ^ k :: Int) :: Double) ** 0.25
in centralComposite k (CCC alpha) nC
-- ---------------------------------------------------------------------------
-- Box-Behnken 計画
-- ---------------------------------------------------------------------------
-- | Box-Behnken design for @k = 3, 4, 5@. Returns @nC@ additional
-- centre points.
--
-- * @k = 3@: 12 corner points + @nC@ centre points.
-- * @k = 4@: 24 corner points + @nC@ centre points.
-- * @k = 5@: 40 corner points + @nC@ centre points.
boxBehnken :: Int -> Int -> [[Double]]
boxBehnken k nC
| k == 3 = bb3 ++ centers
| k == 4 = bb4 ++ centers
| k == 5 = bb5 ++ centers
| otherwise = error
("boxBehnken: only k = 3, 4, 5 supported (got k = "
++ show k ++ ")")
where
centers = replicate nC (replicate k 0)
-- 因子ペア (i, j) (i < j) の二水準組合せで「他は 0」
pairs n = [(i, j) | i <- [0 .. n - 1], j <- [i + 1 .. n - 1]]
pairBlock n (i, j) =
[ [ if x == i then s1
else if x == j then s2
else 0
| x <- [0 .. n - 1] ]
| s1 <- [-1, 1], s2 <- [-1, 1] ]
bb3 = concatMap (pairBlock 3) (pairs 3)
bb4 = concatMap (pairBlock 4) (pairs 4)
bb5 = concatMap (pairBlock 5) (pairs 5)
-- ---------------------------------------------------------------------------
-- 二次モデル
-- ---------------------------------------------------------------------------
-- | Build the design matrix for a quadratic model.
--
-- Each row @[x_1, …, x_k]@ expands to
-- @[1, x_1, …, x_k, x_1², …, x_k², x_1 x_2, x_1 x_3, …, x_{k-1} x_k]@
-- (intercept, main effects, squared terms, upper-triangle interactions).
--
-- Number of columns: @1 + 2k + k(k-1)/2@.
quadraticDesign :: [[Double]] -> LA.Matrix Double
quadraticDesign rows =
let k = if null rows then 0 else length (head rows)
expand row =
let mainE = row
sqE = [x * x | x <- row]
interE = [(row !! i) * (row !! j)
| i <- [0 .. k - 1], j <- [i + 1 .. k - 1]]
in 1 : mainE ++ sqE ++ interE
in LA.fromLists (map expand rows)
-- | Column names for the quadratic-model design (e.g.
-- @[\"b0\", \"x1\", \"x2\", \"x1^2\", \"x2^2\", \"x1*x2\"]@).
quadraticTermNames :: Int -> [Text]
quadraticTermNames k =
["b0"]
++ [T.pack ("x" ++ show i) | i <- [1 .. k]]
++ [T.pack ("x" ++ show i ++ "^2") | i <- [1 .. k]]
++ [T.pack ("x" ++ show i ++ "*x" ++ show j)
| i <- [1 .. k], j <- [i + 1 .. k]]
-- | Quadratic-model fit result.
data QuadFit = QuadFit
{ qfK :: Int -- ^ Number of factors @k@.
, qfBeta :: LA.Vector Double -- ^ Coefficient vector
-- @[b₀, β_main, β_sq, β_int]@.
, qfYHat :: LA.Vector Double -- ^ Fitted values.
, qfR2 :: Double -- ^ R².
} deriving (Show)
-- | Fit a quadratic model by least squares.
fitQuadratic :: [[Double]] -> [Double] -> QuadFit
fitQuadratic xs ys =
let k = if null xs then 0 else length (head xs)
x = quadraticDesign xs
y = LA.fromList ys
beta = LA.flatten (x LA.<\> LA.asColumn y)
yHat = x LA.#> beta
gm = LA.sumElements y / fromIntegral (LA.size y)
ssT = LA.sumElements ((y - LA.scalar gm) ^ (2 :: Int))
ssR = LA.sumElements ((y - yHat) ^ (2 :: Int))
r2 = if ssT == 0 then 0 else 1 - ssR / ssT
in QuadFit k beta yHat r2
-- | Solve analytically for the extremum (saddle / max / min) of the
-- fitted quadratic model.
--
-- Writing @ŷ = b₀ + bᵀx + xᵀ B x@, set @∂ŷ/∂x = 0@ to obtain
-- x* = −½ B⁻¹ b。固有値の符号で性質を判定。
--
-- 戻り値: (x*, predicted_y, eigenvalues)
-- eigenvalues 全部 < 0 → 極大
-- eigenvalues 全部 > 0 → 極小
-- 混在 → 鞍点
optimumPoint :: QuadFit -> ([Double], Double, [Double])
optimumPoint fit =
let k = qfK fit
beta = LA.toList (qfBeta fit)
b0 = head beta
bMain = take k (drop 1 beta)
bSq = take k (drop (1 + k) beta)
bInt = drop (1 + 2 * k) beta
bMat = quadBMatrix fit
bVec = LA.fromList bMain
xStar = LA.toList (LA.scale (-0.5) (LA.inv bMat LA.#> bVec))
yStar = b0
+ sum (zipWith (*) bMain xStar)
+ sum (zipWith (\b x -> b * x * x) bSq xStar)
+ sum [ (bInt !! pairIndex k i j) * (xStar !! i) * (xStar !! j)
| i <- [0 .. k - 1], j <- [i + 1 .. k - 1] ]
eigs = LA.toList (fst (LA.eigSH (LA.sym bMat)))
in (xStar, yStar, eigs)
where
-- (i, j) ペア (i < j) の β_int 配列内のインデックス
pairIndex n i j = sum [n - 1 - p | p <- [0 .. i - 1]] + (j - i - 1)
-- | 二次モデルの @B@ 行列 (@ŷ = b₀ + bᵀx + xᵀ B x@ の 2 次係数)。 対角は β_sq、
-- 非対角は β_int/2 で対称化。 canonical 解析 / 停留点計算の共通部品。
quadBMatrix :: QuadFit -> LA.Matrix Double
quadBMatrix fit =
let k = qfK fit
beta = LA.toList (qfBeta fit)
bSq = take k (drop (1 + k) beta)
bInt = drop (1 + 2 * k) beta
in LA.fromLists
[ [ if i == j then bSq !! i
else
let (lo, hi) = if i < j then (i, j) else (j, i)
idx = pairIndex k lo hi
in (bInt !! idx) / 2
| j <- [0 .. k - 1] ]
| i <- [0 .. k - 1] ]
where pairIndex n i j = sum [n - 1 - p | p <- [0 .. i - 1]] + (j - i - 1)
-- | Canonical 解析。 @B@ 行列の固有分解を返す (固有値, 固有ベクトル) のペア列。
-- 固有値の符号で応答曲面の性質が読める (全負=極大 / 全正=極小 / 混在=鞍点)、
-- 固有ベクトルは canonical 軸 (停留点周りで応答が最も急/緩に動く coded 方向)。
-- ペアは固有値の昇順。 単位はモデルを当てた座標系 (通常 coded)。
canonicalAnalysis :: QuadFit -> [(Double, [Double])]
canonicalAnalysis fit =
let (vals, vecs) = LA.eigSH (LA.sym (quadBMatrix fit))
pairs = zip (LA.toList vals) (map LA.toList (LA.toColumns vecs))
in sortOn fst pairs