hanalyze-0.1.0.0: test/Spec.hs
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeApplications #-}
module Main where
import Test.Hspec
import Hanalyze.Model.GLMM
import Hanalyze.Model.GLM (Family (..), LinkFn (..))
import qualified Data.Vector as V
import qualified Data.Text as T
import qualified Numeric.LinearAlgebra as LA
import Data.List (sort)
import qualified DataFrame as DX
import qualified Hanalyze.Design.Orthogonal as OA
import qualified Hanalyze.Design.Quality as Quality
import qualified Hanalyze.Design.Taguchi as TG
import qualified Hanalyze.Model.LM as LM
import qualified Hanalyze.Model.LM.Diagnostics as LMD
import qualified Hanalyze.DataIO.Preprocess as Pp
import qualified Hanalyze.DataIO.Log as Log
import qualified Hanalyze.DataIO.CSV as CSV
import qualified Hanalyze.DataIO.Convert as Conv
import qualified Hanalyze.DataIO.Health as Health
import qualified Hanalyze.DataIO.Clean as Clean
import qualified Hanalyze.DataIO.Convert as Conv2
import qualified Hanalyze.Stat.Standardize as Std
import qualified Hanalyze.Stat.NumberFormat as NF
import qualified Hanalyze.Stat.Interpolate as Interp
import qualified Hanalyze.Stat.AdaptiveGrid as AG
import qualified Hanalyze.Stat.KernelDist as KD
import qualified Hanalyze.Stat.Cholesky as Chol
import qualified Hanalyze.Stat.QuasiRandom as QR
import qualified Hanalyze.Stat.Test as ST
import qualified Hanalyze.Stat.ClassMetrics as CM
import qualified Hanalyze.Stat.CV as CV
import qualified Hanalyze.Stat.MultipleTesting as MT
import qualified Hanalyze.Stat.Bootstrap as Boot
import qualified Hanalyze.Stat.Effect as Eff
import qualified Hanalyze.Stat.Interpret as Interp
import qualified Hanalyze.Model.PCA as PCA
import qualified Hanalyze.Model.Cluster as Cl
import qualified Hanalyze.Model.DecisionTree as DT
import qualified Hanalyze.Model.TimeSeries as TS
import qualified Hanalyze.Model.Survival as Surv
import qualified Hanalyze.DataIO.Reshape as Reshape
import qualified Hanalyze.Optim.NSGA as NSGA
import qualified System.Random.MWC as MWC
import qualified Hanalyze.Model.Kernel as Kn
import qualified Hanalyze.Model.GP as GP
import qualified Hanalyze.Model.GPRobust as GPR
import qualified Hanalyze.Viz.ReportBuilder as RB
import qualified Data.ByteString as BS
import System.IO.Temp (withSystemTempFile)
import System.IO (hPutStr, hClose)
import qualified Hanalyze.Model.GP as GP
import qualified Hanalyze.Model.GPRobust as GPR
import qualified Hanalyze.Model.RFF as RFF
import qualified Hanalyze.Model.Regularized as Reg
import qualified Hanalyze.Model.Spline as Sp
import qualified Hanalyze.Model.Kernel as K
import qualified Hanalyze.Model.Core as Core
import qualified Hanalyze.Model.GLM as GLM
import qualified Hanalyze.Optim.NelderMead as NM
import qualified Hanalyze.Optim.LBFGS as LBFGS
import qualified Hanalyze.Optim.LineSearch as LS
import qualified Hanalyze.Optim.DifferentialEvolution as DE
import qualified Hanalyze.Optim.CMAES as CMAES
import qualified Hanalyze.Optim.CMAESFull as CMAESF
import qualified Hanalyze.Optim.SimulatedAnnealing as SA
import qualified Hanalyze.Optim.ParticleSwarm as PSO
import qualified Hanalyze.Optim.Constrained as Con
import qualified Hanalyze.Optim.BayesOpt as BO
import qualified Hanalyze.Optim.Common as OC
import qualified System.Random.MWC as MWC
main :: IO ()
main = hspec $ do
describe "Hanalyze.Stat.KernelDist" $ do
let xs = LA.fromLists [[0, 0], [3, 4], [1, 1]] :: LA.Matrix Double
d = KD.pairwiseSqDist xs
-- naive reference
naive m =
let rows = LA.toRows m
n = length rows
in (n LA.>< n)
[ let r = rows !! i - rows !! j in r `LA.dot` r
| i <- [0 .. n - 1], j <- [0 .. n - 1] ]
ref = naive xs
it "returns an n x n matrix with zero diagonal" $ do
LA.rows d `shouldBe` 3
LA.cols d `shouldBe` 3
LA.toList (LA.takeDiag d) `shouldBe` [0, 0, 0]
it "matches the naive reference within 1e-9" $
LA.norm_Inf (d - ref) < 1e-9 `shouldBe` True
it "pairwiseSqDistXY matches reference for cross-matrix" $ do
let ys = LA.fromLists [[0, 0], [1, 0]] :: LA.Matrix Double
dXY = KD.pairwiseSqDistXY xs ys
rxs = LA.toRows xs
rys = LA.toRows ys
ref' = (LA.rows xs LA.>< LA.rows ys)
[ let r = rxs !! i - rys !! j in r `LA.dot` r
| i <- [0 .. LA.rows xs - 1]
, j <- [0 .. LA.rows ys - 1] ]
LA.norm_Inf (dXY - ref') < 1e-9 `shouldBe` True
describe "Hanalyze.Optim.NSGA building blocks" $ do
let mkSol obj = NSGA.Solution
{ NSGA.solDecision = obj -- decision unused here
, NSGA.solObjectives = obj
, NSGA.solViolation = 0
}
s00 = mkSol [0.0, 0.0] -- dominated by no one
s11 = mkSol [1.0, 1.0] -- dominated by s00
s05 = mkSol [0.5, 0.5]
sa = mkSol [0.0, 1.0] -- non-comparable with sb
sb = mkSol [1.0, 0.0]
it "dominates: zero-violation pair uses paretoDominates" $ do
NSGA.dominates s00 s11 `shouldBe` True
NSGA.dominates s11 s00 `shouldBe` False
NSGA.dominates sa sb `shouldBe` False
NSGA.dominates sb sa `shouldBe` False
it "nonDominatedSort: 3-point chain returns 3 fronts" $ do
-- s00 ≻ s05 ≻ s11
let fronts = NSGA.nonDominatedSort [s11, s05, s00]
length fronts `shouldBe` 3
length (head fronts) `shouldBe` 1 -- F_1 = {s00}
it "crowdingDistance: 2-point front keeps both (≤2 → ∞)" $
length (NSGA.crowdingDistance [sa, sb]) `shouldBe` 2
it "crowdingDistance: 3-point front returns all 3 (∞ endpoints + 1 mid)" $ do
let f3 = [mkSol [0.0, 1.0], mkSol [0.5, 0.5], mkSol [1.0, 0.0]]
sorted = NSGA.crowdingDistance f3
length sorted `shouldBe` 3
it "dominationMatrix matches per-pair dominates on a 3-individual pop" $ do
let s00 = NSGA.Solution [0, 0] [0.0, 0.0] 0 -- best
s11 = NSGA.Solution [1, 1] [1.0, 1.0] 0 -- worst
s05 = NSGA.Solution [0.5, 0.5] [0.5, 0.5] 0 -- middle
pop = [s00, s05, s11]
pm = NSGA.fromSolutions pop
mDom = NSGA.dominationMatrix pm
n = length pop
ref = LA.fromLists
[ [ if i == j then 0
else if NSGA.dominates (pop !! i) (pop !! j) then 1
else if NSGA.dominates (pop !! j) (pop !! i) then -1
else 0
| j <- [0 .. n - 1] ]
| i <- [0 .. n - 1] ]
:: LA.Matrix Double
LA.norm_Inf (mDom - ref) `shouldBe` 0
it "polynomialMutation: respects bounds and pMut=0 keeps x" $ do
gen <- MWC.createSystemRandom
let bs = [(0, 1), (0, 1), (0, 1)] :: [(Double, Double)]
x' <- NSGA.polynomialMutation 20 0 bs [0.5, 0.5, 0.5] gen
x' `shouldBe` [0.5, 0.5, 0.5]
y' <- NSGA.polynomialMutation 20 1.0 bs [0.5, 0.5, 0.5] gen
all (\(z, (lo, hi)) -> z >= lo && z <= hi) (zip y' bs)
`shouldBe` True
describe "Hanalyze.Stat.QuasiRandom (Halton)" $ do
it "haltonSequence 10 1 lies in [0, 1) and is a permutation" $ do
let pts = QR.haltonSequence 10 1
length pts `shouldBe` 10
all (\[u] -> u >= 0 && u < 1) pts `shouldBe` True
it "haltonSequence 16 2 covers a 4x4 grid better than random" $ do
-- For n = 16, d = 2 the Halton points should cover all 16
-- 0.25-bins; an iid uniform usually misses some.
let pts = QR.haltonSequence 16 2
binIdx p = (floor (4 * head p) :: Int,
floor (4 * (p !! 1)) :: Int)
unique = length (foldr (\b acc -> if b `elem` acc then acc
else b : acc) [] (map binIdx pts))
unique `shouldSatisfy` (>= 14) -- Halton normally hits ≥ 14 / 16
it "haltonSequenceIn rescales into the supplied bounds" $ do
let bs = [(-2, 2), (10, 20)] :: [(Double, Double)]
pts = QR.haltonSequenceIn 5 bs
all (\[a, b] -> a >= -2 && a < 2
&& b >= 10 && b < 20) pts `shouldBe` True
it "lhsSamples 10 3 lies in [0,1)^3 and fills every per-dim cell" $ do
gen <- MWC.createSystemRandom
pts <- QR.lhsSamples 10 3 gen
length pts `shouldBe` 10
all (\xs -> all (\u -> u >= 0 && u < 1) xs) pts `shouldBe` True
-- For each dim, the 10 points should occupy 10 distinct cells
-- (= floor(10 * u) is a permutation of [0..9]).
let cellsAlong k = map (\xs -> floor (10 * (xs !! k)) :: Int) pts
ok k = sort (cellsAlong k) == [0 .. 9]
ok 0 `shouldBe` True
ok 1 `shouldBe` True
ok 2 `shouldBe` True
it "lhsSamplesIn rescales into bounds" $ do
gen <- MWC.createSystemRandom
let bs = [(-1, 1), (5, 10)] :: [(Double, Double)]
pts <- QR.lhsSamplesIn 8 bs gen
all (\[a, b] -> a >= -1 && a < 1
&& b >= 5 && b < 10) pts `shouldBe` True
describe "Hanalyze.Stat.Cholesky" $ do
let aSPD = LA.fromLists [[4, 2, 1], [2, 5, 3], [1, 3, 6]]
:: LA.Matrix Double
b = LA.asColumn (LA.fromList [1.0, 2.0, 3.0])
it "cholSolve agrees with LA.<\\> on a 3x3 SPD system (1e-9)" $ do
let xC = Chol.cholSolve aSPD b
xR = aSPD LA.<\> b
case xC of
Nothing -> expectationFailure "cholSolve returned Nothing on SPD"
Just xc -> LA.norm_Inf (xc - xR) < 1e-9 `shouldBe` True
it "cholSolveJitter falls back gracefully on a singular matrix" $ do
let aSing = LA.fromLists [[1, 0, 0], [0, 0, 0], [0, 0, 1]]
:: LA.Matrix Double
bSing = LA.asColumn (LA.fromList [1.0, 0.0, 1.0])
x = Chol.cholSolveJitter aSing bSing
LA.rows x `shouldBe` 3 -- did not crash; whatever LSQ gives is fine
it "cholFactor returns Just for SPD and Nothing for non-SPD" $ do
Chol.cholFactor aSPD `shouldSatisfy` (\m -> case m of
Just _ -> True
Nothing -> False)
let aNeg = LA.fromLists [[1, 2], [2, 1]] :: LA.Matrix Double
-- eigenvalues 3 and -1 → not SPD
Chol.cholFactor aNeg `shouldBe` Nothing
describe "Hanalyze.Model.Kernel multi-input (MV)" $ do
-- 2D regression target: y = sin(x1) + 0.5 cos(x2)
let n = 60
h = 0.5
lam = 1e-4
f x1 x2 = sin x1 + 0.5 * cos x2
xs = LA.fromLists
[ [ fromIntegral i / 10
, fromIntegral (n - i) / 10
]
| i <- [0 .. n - 1] ]
ys = LA.asColumn $ LA.fromList
[ f (xs `LA.atIndex` (i, 0)) (xs `LA.atIndex` (i, 1))
| i <- [0 .. n - 1] ]
fit = Kn.kernelRidgeMV Kn.Gaussian h lam xs ys
yhat = Kn.fittedKernelRidgeMV fit
ssErr = LA.sumElements ((ys - yhat) ** 2)
ssTot = let muY = LA.sumElements ys / fromIntegral n
in LA.sumElements ((ys - LA.konst muY (n, 1)) ** 2)
r2 = 1 - ssErr / ssTot
it "achieves R² > 0.95 on a 2D smooth target" $
r2 > 0.95 `shouldBe` True
it "predict at training points equals fitted" $ do
let p = Kn.predictKernelRidgeMV fit xs
LA.norm_Inf (p - yhat) < 1e-9 `shouldBe` True
it "gramMatrixMV matches kernelFromSqDist by element" $ do
let xS = LA.fromLists [[0.0, 0.0], [1.0, 0.0], [0.5, 0.5]] :: LA.Matrix Double
gMV = Kn.gramMatrixMV Kn.Gaussian 1.0 xS
rs = LA.toRows xS
ref = LA.fromLists
[ [ let d = rs !! i - rs !! j
s = (d `LA.dot` d) / (1.0 * 1.0)
in Kn.kernelFromSqDist Kn.Gaussian s
| j <- [0 .. 2] ]
| i <- [0 .. 2] ]
LA.norm_Inf (gMV - ref) < 1e-12 `shouldBe` True
it "Hanalyze.Model.GP MV: 1D input matches legacy 1D fitGP within 1e-6" $ do
let xL = [fromIntegral i / 5 | i <- [0 .. 19 :: Int]]
yL = map sin xL
tL = [0.5, 1.5, 2.5]
mdl = GP.GPModel GP.RBF (GP.GPParams 1.0 1.0 0.05 1.0 Nothing)
legacy = GP.fitGP mdl xL yL tL
xMV = LA.fromLists (map (:[]) xL) :: LA.Matrix Double
yMV = LA.fromList yL
tMV = LA.fromLists (map (:[]) tL) :: LA.Matrix Double
mv = GP.fitGPMV mdl xMV yMV tMV
dMu = LA.norm_Inf
(GP.gpmvMean mv - LA.fromList (GP.gpMean legacy))
dVr = LA.norm_Inf
(GP.gpmvVar mv - LA.fromList (GP.gpVar legacy))
(dMu < 1e-6) `shouldBe` True
(dVr < 1e-6) `shouldBe` True
it "Hanalyze.Model.GP MV: 2D RBF reaches R² > 0.95 with optimized HP" $ do
let nN = 50
gx = [(fromIntegral i / 10, fromIntegral (nN - i) / 10)
| i <- [0 .. nN - 1 :: Int]]
ftn (x1, x2) = sin x1 + 0.5 * cos x2
xMV = LA.fromLists [ [a, b] | (a, b) <- gx ] :: LA.Matrix Double
yMV = LA.fromList [ ftn p | p <- gx ]
p0 = GP.GPParams 1.0 1.0 0.01 1.0 Nothing
po = GP.optimizeGPMV GP.RBF xMV yMV p0
mdl = GP.GPModel GP.RBF po
res = GP.fitGPMV mdl xMV yMV xMV
mu = GP.gpmvMean res
y = yMV
ss = LA.sumElements ((y - mu) ** 2)
mY = LA.sumElements y / fromIntegral nN
st = LA.sumElements ((y - LA.konst mY nN) ** 2)
r2 = 1 - ss / st
(r2 > 0.95) `shouldBe` True
it "Hanalyze.Model.GPRobust MV: 1D input matches legacy fitGPRobust" $ do
let xL = [fromIntegral i / 5 | i <- [0 .. 14 :: Int]]
yL = map sin xL
tL = [0.5, 1.5, 2.5]
ker = GP.RBF
ps = GP.GPParams 1.0 1.0 0.05 1.0 Nothing
lik = GPR.RGaussian 0.1
legFit = GPR.fitGPRobust ker ps lik xL yL
legacy = GPR.predictGPRobust legFit tL
legM = LA.fromList (map fst legacy)
xMV = LA.fromLists (map (:[]) xL) :: LA.Matrix Double
yMV = LA.fromList yL
tMV = LA.fromLists (map (:[]) tL) :: LA.Matrix Double
mvFit = GPR.fitGPRobustMV ker ps lik xMV yMV
(mvM, _) = GPR.predictGPRobustMV mvFit tMV
LA.norm_Inf (mvM - legM) < 1e-6 `shouldBe` True
it "MV gramMatrix on a single-column input agrees with kernelFromSqDist" $ do
let xs1 = LA.fromLists [[fromIntegral i / 5] | i <- [0 .. 19 :: Int]]
:: LA.Matrix Double
gMV2 = Kn.gramMatrixMV Kn.Gaussian 0.4 xs1
ref = LA.fromLists
[ [ let xi = xs1 `LA.atIndex` (i, 0)
xj = xs1 `LA.atIndex` (j, 0)
d = xi - xj
in Kn.kernelFromSqDist Kn.Gaussian (d * d / (0.4 * 0.4))
| j <- [0 .. 19] ]
| i <- [0 .. 19] ]
LA.norm_Inf (gMV2 - ref) < 1e-12 `shouldBe` True
describe "Hanalyze.Model.GLMM" $ do
-- Dataset: 3 groups × 4 obs, strong between-group signal, weak within-group noise.
-- True: β₀≈5, β₁≈0, u_A≈2, u_B≈0, u_C≈-2, σ²_u≈4, σ²≈small → ICC≈high.
let df = DX.fromNamedColumns
[ ("x", DX.fromList ([1,2,3,4, 1,2,3,4, 1,2,3,4] :: [Double]))
, ("y", DX.fromList ([7.1,6.9,7.0,7.0, 5.0,4.9,5.1,5.0, 3.0,2.9,3.1,3.0] :: [Double]))
, ("group", DX.fromList (["A","A","A","A","B","B","B","B","C","C","C","C"] :: [T.Text])) ]
res = fitLMEDataFrame [("x", 1)] "group" "y" df
it "returns Just for valid input" $
res `shouldSatisfy` (\r -> case r of { Just _ -> True; Nothing -> False })
it "ICC is in [0, 1]" $
maybe (expectationFailure "expected Just") (\r ->
glmmICC r `shouldSatisfy` (\v -> v >= 0 && v <= 1)) res
it "ICC is high for strongly grouped data" $
maybe (expectationFailure "expected Just") (\r ->
glmmICC r `shouldSatisfy` (> 0.9)) res
it "random variance is positive" $
maybe (expectationFailure "expected Just") (\r ->
glmmRandVar r `shouldSatisfy` (> 0)) res
it "residual variance is positive" $
maybe (expectationFailure "expected Just") (\r ->
glmmResidVar r `shouldSatisfy` (> 0)) res
it "BLUP count equals number of groups" $
maybe (expectationFailure "expected Just") (\r ->
V.length (glmmBLUPs r) `shouldBe` 3) res
it "group labels are sorted" $
case res of
Just r -> glmmGroups r `shouldBe` V.fromList ["A","B","C"]
Nothing -> expectationFailure "expected Just"
it "returns Nothing for missing column" $
fitLMEDataFrame [("x", 1)] "group" "missing" df
`shouldSatisfy` (\r -> case r of { Nothing -> True; Just _ -> False })
describe "Hanalyze.Model.GLMM (non-Gaussian)" $ do
-- Binomial GLMM: 3 hospitals, binary outcome (treatment success)
-- Strong hospital effect; within each hospital, dose → higher success rate.
-- True: u_A ≈ +1, u_B ≈ 0, u_C ≈ -1 (on logit scale)
let dfBin = DX.fromNamedColumns
[ ("dose", DX.fromList ([1,2,3,4,5, 1,2,3,4,5, 1,2,3,4,5] :: [Double]))
, ("success", DX.fromList ([1,1,1,1,1, 1,1,0,1,0, 0,0,0,1,0] :: [Double]))
, ("hospital", DX.fromList (["A","A","A","A","A","B","B","B","B","B","C","C","C","C","C"] :: [T.Text])) ]
resBin = fitGLMMDataFrame Binomial Logit [("dose", 1)] "hospital" "success" dfBin
it "Binomial GLMM returns Just" $
resBin `shouldSatisfy` (\r -> case r of { Just _ -> True; Nothing -> False })
it "Binomial ICC in [0, 1]" $
maybe (expectationFailure "expected Just") (\r ->
glmmICC r `shouldSatisfy` (\v -> v >= 0 && v <= 1)) resBin
it "Binomial σ²_u is positive" $
maybe (expectationFailure "expected Just") (\r ->
glmmRandVar r `shouldSatisfy` (> 0)) resBin
-- Poisson GLMM: 3 regions, count outcome (events per month)
-- True: β₀ on log scale ≈ 2 (≈7 events baseline), u differs by region.
let dfPois = DX.fromNamedColumns
[ ("time", DX.fromList ([1,2,3,4,5, 1,2,3,4,5, 1,2,3,4,5] :: [Double]))
, ("count", DX.fromList ([15,18,22,20,25, 7,9,8,10,11, 2,3,2,4,3] :: [Double]))
, ("region", DX.fromList (["X","X","X","X","X","Y","Y","Y","Y","Y","Z","Z","Z","Z","Z"] :: [T.Text])) ]
resPois = fitGLMMDataFrame Poisson Log [("time", 1)] "region" "count" dfPois
it "Poisson GLMM returns Just" $
resPois `shouldSatisfy` (\r -> case r of { Just _ -> True; Nothing -> False })
it "Poisson σ²_u is positive" $
maybe (expectationFailure "expected Just") (\r ->
glmmRandVar r `shouldSatisfy` (> 0)) resPois
it "Poisson ICC in [0, 1]" $
maybe (expectationFailure "expected Just") (\r ->
glmmICC r `shouldSatisfy` (\v -> v >= 0 && v <= 1)) resPois
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Design.Orthogonal" $ do
it "L4 has 4 runs and 3 columns" $ do
OA.oaRuns OA.l4 `shouldBe` 4
OA.oaFactors OA.l4 `shouldBe` 3
length (OA.oaTable OA.l4) `shouldBe` 4
it "L8 has 8 runs and 7 columns" $ do
OA.oaRuns OA.l8 `shouldBe` 8
OA.oaFactors OA.l8 `shouldBe` 7
it "L9 has 9 runs and 4 columns at 3 levels each" $ do
OA.oaRuns OA.l9 `shouldBe` 9
OA.oaFactors OA.l9 `shouldBe` 4
OA.oaLevels OA.l9 `shouldBe` [3, 3, 3, 3]
it "L18 has 18 runs and 8 columns (1×2 + 7×3)" $ do
OA.oaRuns OA.l18 `shouldBe` 18
OA.oaLevels OA.l18 `shouldBe` 2 : replicate 7 3
it "L8 columns are balanced (each level appears 4 times)" $ do
let table = OA.oaTable OA.l8
colJ j = [ row !! j | row <- table ]
mapM_ (\j -> do
let cs = colJ j
length (filter (== 1) cs) `shouldBe` 4
length (filter (== 2) cs) `shouldBe` 4) [0 .. 6]
it "L8 column pairs are pairwise orthogonal" $ do
let table = OA.oaTable OA.l8
colJ j = [ row !! j | row <- table ]
pairCount j1 j2 a b =
length (filter id (zipWith (\x y -> x == a && y == b)
(colJ j1) (colJ j2)))
-- For 2-level orthogonality: each pair (1,1)/(1,2)/(2,1)/(2,2) must appear equally
mapM_ (\(j1, j2) -> do
pairCount j1 j2 1 1 `shouldBe` 2
pairCount j1 j2 1 2 `shouldBe` 2
pairCount j1 j2 2 1 `shouldBe` 2
pairCount j1 j2 2 2 `shouldBe` 2)
[(0,1),(0,2),(0,3),(1,2),(1,3),(2,3)]
it "lookupOA finds standard arrays case-insensitively" $ do
OA.oaName <$> OA.lookupOA "L9" `shouldBe` Just "L9(3^4)"
OA.oaName <$> OA.lookupOA "l9" `shouldBe` Just "L9(3^4)"
OA.oaName <$> OA.lookupOA "L99" `shouldBe` Nothing
it "assignFactors fills levels correctly for L4" $ do
let specs = [ OA.FactorSpec "A" [OA.LText "lo", OA.LText "hi"]
, OA.FactorSpec "B" [OA.LNumeric 0, OA.LNumeric 1]
]
case OA.assignFactors OA.l4 specs of
Right ad -> do
length (OA.adRows ad) `shouldBe` 4
map length (OA.adRows ad) `shouldBe` [2, 2, 2, 2]
Left e -> expectationFailure (show e)
it "assignFactors rejects too many factors" $ do
let specs = replicate 5 (OA.FactorSpec "X" [OA.LNumeric 1, OA.LNumeric 2])
OA.assignFactors OA.l4 specs `shouldSatisfy`
\r -> case r of { Left _ -> True; Right _ -> False }
it "assignFactors rejects level count mismatch" $ do
let specs = [ OA.FactorSpec "A" [OA.LText "x"] ] -- only 1 level, L4 needs 2
OA.assignFactors OA.l4 specs `shouldSatisfy`
\r -> case r of { Left _ -> True; Right _ -> False }
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Design.Taguchi" $ do
it "smaller-the-better SN: lower y → higher η" $ do
let etaSmall = TG.snRatio TG.SmallerBetter [0.5, 0.5, 0.5]
etaLarge = TG.snRatio TG.SmallerBetter [5.0, 5.0, 5.0]
etaSmall `shouldSatisfy` (> etaLarge)
it "larger-the-better SN: higher y → higher η" $ do
let etaLarge = TG.snRatio TG.LargerBetter [10, 10, 10]
etaSmall = TG.snRatio TG.LargerBetter [1, 1, 1]
etaLarge `shouldSatisfy` (> etaSmall)
it "nominal-the-best SN: high mean / low var → high η" $ do
let highSN = TG.snRatio TG.NominalBest [10, 10.01, 9.99, 10]
lowSN = TG.snRatio TG.NominalBest [1, 4, 7, 10]
highSN `shouldSatisfy` (> lowSN)
it "nominal-target SN: closer to target → higher η" $ do
let closer = TG.snRatio (TG.NominalBestTarget 5) [4.9, 5.0, 5.1]
farther = TG.snRatio (TG.NominalBestTarget 5) [3, 5, 7]
closer `shouldSatisfy` (> farther)
it "snRatio on empty list is 0" $
TG.snRatio TG.SmallerBetter [] `shouldBe` 0
it "snRatioRows produces same length as input" $ do
let yMatrix = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]
length (TG.snRatioRows TG.SmallerBetter yMatrix) `shouldBe` 3
it "analyzeSN gives one FactorEffect per assigned factor" $ do
let specs = [ OA.FactorSpec "A" [OA.LText "lo", OA.LText "hi"]
, OA.FactorSpec "B" [OA.LText "lo", OA.LText "hi"]
, OA.FactorSpec "C" [OA.LText "lo", OA.LText "hi"]
]
Right ad = OA.assignFactors OA.l4 specs
fes = TG.analyzeSN ad [10, 20, 30, 40]
length fes `shouldBe` 3
map TG.feFactor fes `shouldBe` ["A", "B", "C"]
mapM_ (\fe -> length (TG.feSNByLevel fe) `shouldBe` 2) fes
it "optimalLevels picks max-SN level per factor" $ do
let specs = [ OA.FactorSpec "A" [OA.LText "lo", OA.LText "hi"]
, OA.FactorSpec "B" [OA.LText "lo", OA.LText "hi"]
, OA.FactorSpec "C" [OA.LText "lo", OA.LText "hi"]
]
Right ad = OA.assignFactors OA.l4 specs
-- L4 row 1 ("hi" for A,B,C) has SN=100; others 0
sns = [0, 0, 0, 100]
fes = TG.analyzeSN ad sns
opts = TG.optimalLevels fes
length opts `shouldBe` 3
-- Row 4 has all "hi" by L4 structure (2,2,1) — verify each factor's best
mapM_ (\(_, _, eta) -> eta `shouldSatisfy` (>= 0)) opts
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.DataIO.Preprocess" $ do
let dfNA = DX.fromNamedColumns
[ ("group", DX.fromList (["A","B","A","B","C"] :: [T.Text]))
, ("x", DX.fromList (["1","NA","3","","5"] :: [T.Text]))
, ("y", DX.fromList ([10, 20, 30, 40, 50] :: [Double]))
]
it "isNAString detects standard NA strings" $ do
Pp.isNAString "NA" `shouldBe` True
Pp.isNAString "N/A" `shouldBe` True
Pp.isNAString "null" `shouldBe` True
Pp.isNAString "" `shouldBe` True
Pp.isNAString " " `shouldBe` True
Pp.isNAString "valid" `shouldBe` False
it "countMissing counts NAs in Text columns; numeric is 0" $ do
let counts = Pp.countMissing dfNA
lookup "x" counts `shouldBe` Just 2
lookup "y" counts `shouldBe` Just 0
lookup "group" counts `shouldBe` Just 0
it "dropMissingRows removes rows with NA in target columns" $ do
let df' = Pp.dropMissingRows ["x"] dfNA
(n, _) = DX.dimensions df'
n `shouldBe` 3 -- only rows with x ∈ {"1","3","5"} remain
it "imputeMean converts Text/NA column to Double with mean fill" $ do
case Pp.imputeMean "x" dfNA of
Just df' -> do
let xs = DX.columnAsList (DX.col @Double "x") df'
length xs `shouldBe` 5
-- mean of [1, 3, 5] = 3
(xs !! 1) `shouldBe` 3.0 -- was "NA"
(xs !! 3) `shouldBe` 3.0 -- was ""
Nothing -> expectationFailure "imputeMean failed"
it "selectColumns retains only listed columns" $ do
let df' = Pp.selectColumns ["y", "group"] dfNA
DX.columnNames df' `shouldMatchList` ["y", "group"]
it "filterRowsByNumeric filters numeric column" $ do
let df' = Pp.filterRowsByNumeric "y" (>= 30) dfNA
(n, _) = DX.dimensions df'
n `shouldBe` 3
it "mapNumeric applies a unary function" $ do
let df' = Pp.mapNumeric "y" (* 2) dfNA
xs = DX.columnAsList (DX.col @Double "y") df'
xs `shouldBe` [20, 40, 60, 80, 100]
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.DataIO.Preprocess (groupBy)" $ do
let dfGrp = DX.fromNamedColumns
[ ("group", DX.fromList (["A","B","A","B","A","C"] :: [T.Text]))
, ("y", DX.fromList ([1, 4, 3, 6, 5, 10] :: [Double]))
]
it "groupByMean computes per-group mean" $ do
case Pp.groupByMean "group" "y" dfGrp of
Just df' -> do
let (n, _) = DX.dimensions df'
n `shouldBe` 3
let gs = DX.columnAsList (DX.col @T.Text "group") df'
vs = DX.columnAsList (DX.col @Double "y") df'
pairs = zip gs vs
lookup "A" pairs `shouldBe` Just 3.0 -- (1+3+5)/3
lookup "B" pairs `shouldBe` Just 5.0 -- (4+6)/2
lookup "C" pairs `shouldBe` Just 10.0
Nothing -> expectationFailure "groupByMean failed"
it "groupBySum computes per-group sum" $ do
case Pp.groupBySum "group" "y" dfGrp of
Just df' -> do
let gs = DX.columnAsList (DX.col @T.Text "group") df'
vs = DX.columnAsList (DX.col @Double "y") df'
pairs = zip gs vs
lookup "A" pairs `shouldBe` Just 9.0
lookup "B" pairs `shouldBe` Just 10.0
Nothing -> expectationFailure "groupBySum failed"
it "groupByCount counts rows per group" $ do
case Pp.groupByCount "group" dfGrp of
Just df' -> do
let gs = DX.columnAsList (DX.col @T.Text "group") df'
vs = DX.columnAsList (DX.col @Double "count") df'
pairs = zip gs vs
lookup "A" pairs `shouldBe` Just 3.0
lookup "B" pairs `shouldBe` Just 2.0
lookup "C" pairs `shouldBe` Just 1.0
Nothing -> expectationFailure "groupByCount failed"
it "groupByMin/Max return correct extremes" $ do
case Pp.groupByMin "group" "y" dfGrp of
Just dfMin -> do
let gs = DX.columnAsList (DX.col @T.Text "group") dfMin
vs = DX.columnAsList (DX.col @Double "y") dfMin
pairs = zip gs vs
lookup "A" pairs `shouldBe` Just 1.0
lookup "B" pairs `shouldBe` Just 4.0
Nothing -> expectationFailure "groupByMin failed"
case Pp.groupByMax "group" "y" dfGrp of
Just dfMax -> do
let gs = DX.columnAsList (DX.col @T.Text "group") dfMax
vs = DX.columnAsList (DX.col @Double "y") dfMax
pairs = zip gs vs
lookup "A" pairs `shouldBe` Just 5.0
lookup "B" pairs `shouldBe` Just 6.0
Nothing -> expectationFailure "groupByMax failed"
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Stat.NumberFormat" $ do
it "0 → '0.00'" $ NF.fmtNum 0 `shouldBe` "0.00"
it "中域 (0.01..999) は固定小数点 2 桁" $ do
NF.fmtNum 0.91 `shouldBe` "0.91"
NF.fmtNum 12.34 `shouldBe` "12.34"
NF.fmtNum 998.7 `shouldBe` "998.70"
it "巨大値は指数表記" $ do
NF.fmtNum 1.10e13 `shouldBe` "1.10E+13"
NF.fmtNum 1234.5 `shouldBe` "1.23E+3"
it "極小値は指数表記" $ do
NF.fmtNum 3.057e-24 `shouldBe` "3.06E-24"
NF.fmtNum 0.0099 `shouldBe` "9.90E-3"
it "負の値" $ do
NF.fmtNum (-12.34) `shouldBe` "-12.34"
NF.fmtNum (-1.5e10) `shouldBe` "-1.50E+10"
it "非有限値" $ do
NF.fmtNum (0/0) `shouldBe` "NaN"
NF.fmtNum (1/0) `shouldBe` "+Inf"
NF.fmtNum (-1/0) `shouldBe` "-Inf"
describe "Hanalyze.Stat.Standardize" $ do
let xMat = LA.fromLists [[1, 100], [2, 200], [3, 300], [4, 400], [5, 500]] :: LA.Matrix Double
s = Std.fitStandardizer xMat
it "fit: 各列の μ が一致" $
Std.stMu s `shouldSatisfy`
(\ms -> length ms == 2
&& abs (ms !! 0 - 3) < 1e-9
&& abs (ms !! 1 - 300) < 1e-9)
it "fit: 各列の σ が不偏分散の平方根 (n-1 正規化)" $
Std.stSd s `shouldSatisfy`
(\ss -> length ss == 2
&& abs (ss !! 0 - sqrt 2.5) < 1e-9
&& abs (ss !! 1 - sqrt 25000) < 1e-9)
it "apply 後は各列 mean≈0, std≈1" $ do
let x' = Std.applyStandardizer s xMat
c0 = LA.toColumns x' !! 0
c1 = LA.toColumns x' !! 1
mn v = LA.sumElements v / fromIntegral (LA.size v)
abs (mn c0) `shouldSatisfy` (< 1e-9)
abs (mn c1) `shouldSatisfy` (< 1e-9)
it "unapply で元の値に戻る" $ do
let x' = Std.applyStandardizer s xMat
x'' = Std.unapplyStandardizer s x'
d = LA.norm_2 (xMat - x'') :: Double
d `shouldSatisfy` (< 1e-9)
it "定数列 (std=0) は std=1 にフォールバック (中央化のみ)" $ do
let constMat = LA.fromLists [[7, 1], [7, 2], [7, 3]] :: LA.Matrix Double
s2 = Std.fitStandardizer constMat
abs (Std.stSd s2 !! 0 - 1.0) `shouldSatisfy` (< 1e-12)
let x' = Std.applyStandardizer s2 constMat
c0 = LA.toColumns x' !! 0
abs (LA.sumElements c0) `shouldSatisfy` (< 1e-9)
describe "Hanalyze.Stat.Interpolate" $ do
let pts = [(0,0), (1,1), (2,4), (3,9), (4,16)] -- y = x^2
it "Linear: 観測点で原値を厳密に再現" $ do
let f = Interp.interp1d Interp.Linear pts
mapM_ (\(x, y) -> abs (f x - y) `shouldSatisfy` (< 1e-12)) pts
it "NaturalSpline: 観測点で原値を厳密に再現、中間点で線形より精度高い" $ do
let fl = Interp.interp1d Interp.Linear pts
fs = Interp.interp1d Interp.NaturalSpline pts
true x = x * x
xMid = 1.5
errL = abs (fl xMid - true xMid)
errS = abs (fs xMid - true xMid)
mapM_ (\(x, y) -> abs (fs x - y) `shouldSatisfy` (< 1e-9)) pts
errS `shouldSatisfy` (< errL)
it "PCHIP: 単調データ ([0,1,2,4,8,16]) で出力も単調" $ do
let mp = [(0,0), (1,1), (2,2), (3,4), (4,8), (5,16)]
fp = Interp.interp1d Interp.PCHIP mp
ys = map fp [0, 0.1 .. 5]
and (zipWith (<=) ys (tail ys)) `shouldBe` True
it "PCHIP: 観測点で原値を厳密に再現" $ do
let fp = Interp.interp1d Interp.PCHIP pts
mapM_ (\(x, y) -> abs (fp x - y) `shouldSatisfy` (< 1e-9)) pts
describe "Hanalyze.Stat.AdaptiveGrid" $ do
it "uniformGrid: 端点を含み等間隔" $ do
let g = AG.uniformGrid 5 0 4
g `shouldBe` [0, 1, 2, 3, 4]
it "Adaptive: 急激な変化付近に grid が集中する (step 関数)" $ do
-- y は z=0..1 でほぼ平坦、z=1..2 で急激に変化、z=2..3 で再び平坦
let pts1 = [(z, if z < 1 then 0 else if z < 2 then 10*(z-1) else 10) | z <- [0, 0.1 .. 3]]
spec = (AG.defaultGridSpec 30) { AG.gsKind = AG.Adaptive }
g = AG.makeGrid [pts1] (0, 3) spec
-- 中央領域 [1, 2] にある grid 点数 vs 端領域 [0,1] の grid 点数
midN = length (filter (\z -> z >= 1 && z <= 2) g)
edgeN = length (filter (\z -> z < 1) g)
length g `shouldBe` 30
midN `shouldSatisfy` (> edgeN)
it "Uniform: 端点 + 等間隔" $ do
let g = AG.makeGrid [] (0, 1) ((AG.defaultGridSpec 6) { AG.gsKind = AG.Uniform })
length g `shouldBe` 6
head g `shouldBe` 0
last g `shouldBe` 1
it "N < 10 で adaptive 指定でも uniform にフォールバック" $ do
let pts1 = [(z, sin z) | z <- [0, 0.1 .. 3]]
spec = (AG.defaultGridSpec 5) { AG.gsKind = AG.Adaptive }
g = AG.makeGrid [pts1] (0, 3) spec
gU = AG.uniformGrid 5 0 3
g `shouldBe` gU
describe "Hanalyze.Model.RFF (multivariate, Phase B-RFF)" $ do
it "logMarginalLikRBFMV: 既知 ℓ で最大化される (合成データで)" $ do
-- y = sin(x) (1D) で ℓ をスキャンし、データの z-score 後の長さスケールに
-- 近い値で marg-lik が最大になることを確認。
let xs = [0.0, 0.3 .. 6.0]
ys = map sin xs
xMat = LA.fromLists [[x] | x <- xs]
yV = LA.fromList ys
ells = [0.05, 0.2, 0.5, 1.0, 2.0, 5.0]
mliks = [ RFF.logMarginalLikRBFMV xMat yV ell 1.0 0.05 | ell <- ells ]
best = snd (maximum (zip mliks ells))
best `shouldSatisfy` (\b -> b >= 0.2 && b <= 2.0)
it "loocvRFFRidgeMV: λ → ∞ で残差ベース LOOCV が増える、適度な λ で最小" $ do
let xs = [0.0, 0.3 .. 6.0]
ys = map sin xs
xMat = LA.fromLists [[x] | x <- xs]
yV = LA.fromList ys
gen <- MWC.createSystemRandom
feats <- RFF.sampleRFFRBFMV 1 64 0.5 1.0 gen
let lamSmall = RFF.loocvRFFRidgeMV feats xMat yV 1e-2
lamHuge = RFF.loocvRFFRidgeMV feats xMat yV 1e6
lamSmall `shouldSatisfy` (< lamHuge)
it "gridSearchLOOCVRBFMV: ℓ/λ を自動探索して LOOCV が小さくなる" $ do
let xs = [0.0, 0.5 .. 10.0]
ys = [ sin (x/2) | x <- xs ]
xMat = LA.fromLists [[x] | x <- xs]
yV = LA.fromList ys
gen <- MWC.createSystemRandom
res <- RFF.gridSearchLOOCVRBFMV 1 100 xMat yV (Just (4, 8)) gen
RFF.lcLOOCV res `shouldSatisfy` (< 1.0)
RFF.lcEll res `shouldSatisfy` (> 0)
it "maximizeMarginalLikRBFMV: 雑音ありデータで mlik が改善する" $ do
let xs = [0.0, 0.5 .. 10.0]
ys = [ sin (x/2) + 0.05 * (fromIntegral i / 21) - 0.025
| (i, x) <- zip [0::Int ..] xs ]
xMat = LA.fromLists [[x] | x <- xs]
yV = LA.fromList ys
res = RFF.maximizeMarginalLikRBFMV xMat yV (Just (8, 4, 4))
-- 最適 mlik > 任意の "ヘンな" 値 (ℓ=100, σ_n=10) より高い
weak = RFF.logMarginalLikRBFMV xMat yV 100 1.0 10.0
RFF.mlLogMlik res `shouldSatisfy` (> weak)
RFF.mlEll res `shouldSatisfy` (> 0)
RFF.mlSigmaN res `shouldSatisfy` (> 0)
it "rffRidgeMV: y = x1 * t を完全にフィット" $ do
let xs = [(x1, t) | x1 <- [1, 2, 3, 5, 7], t <- [1..10]]
xss = [[x1, t] | (x1, t) <- xs]
ys = [x1 * t | (x1, t) <- xs]
xMat = LA.fromLists xss
gen <- MWC.createSystemRandom
feats <- RFF.sampleRFFRBFMV 2 256 1.0 1.0 gen
let fit = RFF.rffRidgeMV feats xMat ys 0.001
yhat = RFF.predictRFFRidgeMV fit xMat
rmse = sqrt (sum (zipWith (\a b -> (a-b)*(a-b)) ys yhat)
/ fromIntegral (length ys))
rmse `shouldSatisfy` (< 1.0)
describe "Hanalyze.Model.RFF" $ do
it "feature matrix has correct shape" $ do
gen <- MWC.createSystemRandom
feats <- RFF.sampleRFFRBF 50 1.0 1.0 gen
RFF.rffDim feats `shouldBe` 50
let phi = RFF.rffFeatures feats [0.0, 1.0, 2.0]
-- phi is n × D = 3 × 50
V.length (V.fromList [0::Int]) `shouldBe` 1 -- placeholder for typing
-- We can't easily check matrix shape without hmatrix import here,
-- so just ensure the function doesn't crash.
length (RFF.rffOmegas feats) `shouldSatisfy` (== 50)
let _ = phi
return ()
it "RFF Ridge fits y ≈ x reasonably" $ do
gen <- MWC.createSystemRandom
feats <- RFF.sampleRFFRBF 100 1.0 1.0 gen
let xs = [0.0, 0.1 .. 1.0]
ys = map (\x -> 2 * x + 0.5) xs
fit = RFF.rffRidge feats xs ys 0.001
yhat = RFF.predictRFFRidge fit xs
rmse = sqrt (sum [ (y - yh) ^ (2 :: Int)
| (y, yh) <- zip ys yhat ]
/ fromIntegral (length ys))
rmse `shouldSatisfy` (< 0.5)
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Model.GPRobust" $ do
it "Cauchy GP is more accurate than Gaussian GP under outliers" $ do
let trueF x = sin x
xs = [0.0, 0.5 .. 6.0]
cleanY = map trueF xs
-- Inject outlier at index 5
ys = zipWith (\i y -> if i == 5 then y + 5 else y)
[0::Int ..] cleanY
hp = GP.GPParams 1.0 1.0 0.05 1.0 Nothing
gpRes = GP.fitGP (GP.GPModel GP.RBF hp) xs ys xs
gaussRMSE = sqrt (sum [ (a - b) ^ (2::Int)
| (a, b) <- zip cleanY (GP.gpMean gpRes) ]
/ fromIntegral (length xs))
cauchyFit = GPR.fitGPRobust GP.RBF hp (GPR.RCauchy 0.5) xs ys
cauchyPred = GPR.predictGPRobust cauchyFit xs
cauchyRMSE = sqrt (sum [ (a - b) ^ (2::Int)
| (a, (b, _)) <- zip cleanY cauchyPred ]
/ fromIntegral (length xs))
cauchyRMSE `shouldSatisfy` (< gaussRMSE)
it "IRLS converges in finite iterations" $ do
let xs = [0.0, 1.0, 2.0, 3.0, 4.0]
ys = [0.1, 1.05, 1.95, 2.9, 4.1]
hp = GP.GPParams 1.0 1.0 0.1 1.0 Nothing
fit = GPR.fitGPRobust GP.RBF hp (GPR.RStudentT 4 0.5) xs ys
GPR.rgpIters fit `shouldSatisfy` (\n -> n > 0 && n <= 50)
describe "Hanalyze.DataIO.Log" $ do
it "Monoid: noLog <> r == r" $ do
let r = Log.logReport (Log.mkWarn "W001" "msg" Nothing)
Log.entries (Log.noLog <> r) `shouldBe` Log.entries r
Log.entries (r <> Log.noLog) `shouldBe` Log.entries r
it "addEntry appends" $ do
let r0 = Log.logReport (Log.mkInfo "I001" "first" Nothing)
r1 = Log.addEntry (Log.mkWarn "W001" "second" (Just "ヒント")) r0
length (Log.entries r1) `shouldBe` 2
Log.lgSev (last (Log.entries r1)) `shouldBe` Log.Warn
Log.lgHint (last (Log.entries r1)) `shouldBe` Just "ヒント"
it "hasErrors / hasWarnings detect severity" $ do
let rW = Log.logReport (Log.mkWarn "W" "w" Nothing)
rE = Log.logReport (Log.mkErr "E" "e" Nothing)
Log.hasWarnings rW `shouldBe` True
Log.hasErrors rW `shouldBe` False
Log.hasErrors (rW <> rE) `shouldBe` True
it "severityCount counts each level" $ do
let r = Log.logReport (Log.mkInfo "I" "i" Nothing)
<> Log.logReport (Log.mkWarn "W1" "w" Nothing)
<> Log.logReport (Log.mkWarn "W2" "w" Nothing)
<> Log.logReport (Log.mkErr "E" "e" Nothing)
Log.severityCount Log.Info r `shouldBe` 1
Log.severityCount Log.Warn r `shouldBe` 2
Log.severityCount Log.Err r `shouldBe` 1
it "prettyEntry: includes code, message, hint" $ do
let s = Log.prettyEntry (Log.mkWarn "W042" "壊れている" (Just "助言"))
T.isInfixOf "[WARN]" s `shouldBe` True
T.isInfixOf "W042" s `shouldBe` True
T.isInfixOf "壊れている" s `shouldBe` True
T.isInfixOf "助言" s `shouldBe` True
describe "Hanalyze.DataIO.CSV.loadAutoSafe" $ do
it "Empty file → Left, no exception" $
withSystemTempFile "ha-empty.csv" $ \fp h -> do
hPutStr h ""
hClose h
r <- CSV.loadAutoSafe fp
case r of
Left msg -> T.isInfixOf "Empty" (T.pack msg) `shouldBe` True
Right _ -> expectationFailure "expected Left for empty file"
it "Header-only file → Left" $
withSystemTempFile "ha-hdr.csv" $ \fp h -> do
hPutStr h "x,y,z\n"
hClose h
r <- CSV.loadAutoSafe fp
case r of
Left msg -> T.isInfixOf "header" (T.pack msg) `shouldBe` True
Right _ -> expectationFailure "expected Left for header-only file"
it "Valid CSV → Right with empty log by default" $
withSystemTempFile "ha-ok.csv" $ \fp h -> do
hPutStr h "x,y\n1,2\n3,4\n"
hClose h
r <- CSV.loadAutoSafe fp
case r of
Left msg -> expectationFailure ("unexpected Left: " ++ msg)
Right (_, lg) -> Log.entries lg `shouldBe` []
describe "Hanalyze.DataIO.Convert deep-eval" $ do
it "getMaybeTextVec on text column with mixed NA strings returns Just" $
withSystemTempFile "ha-na.csv" $ \fp h -> do
-- 複数 NA 表現を混ぜると Hackage は Maybe Text 列として保持する。
-- ヘッダ判定で n/a / null は欠損扱い → null bitmap が立つ。
hPutStr h "id,score\n1,A\n2,n/a\n3,null\n4,B\n5,-\n"
hClose h
r <- CSV.loadAutoSafe fp
case r of
Right (df, _) -> case Conv.getMaybeTextVec "score" df of
Just v -> length (V.toList v) `shouldBe` 5
Nothing -> expectationFailure "getMaybeTextVec returned Nothing"
Left msg -> expectationFailure ("load failed: " ++ msg)
it "getDoubleVec returns Nothing without crashing on NA-mixed numeric column" $
withSystemTempFile "ha-na2.csv" $ \fp h -> do
hPutStr h "id,score\n1,85\n2,NA\n3,92\n"
hClose h
r <- CSV.loadAutoSafe fp
case r of
Right (df, _) -> Conv.getDoubleVec "score" df `shouldBe` Nothing
Left msg -> expectationFailure ("load failed: " ++ msg)
describe "Hanalyze.DataIO.Health" $ do
it "W001: ヘッダ無し疑い (列名が全て数値)" $ do
let df = DX.insertColumn "1.0" (DX.fromList ([2.0, 4.0] :: [Double]))
$ DX.insertColumn "2.0" (DX.fromList ([4.1, 8.0] :: [Double]))
$ DX.empty
codes = map Log.lgCode (Log.entries (Health.detectHeaderless df))
codes `shouldContain` ["W001"]
it "W001 は通常ヘッダでは発火しない" $ do
let df = DX.insertColumn "x" (DX.fromList ([1.0, 2.0] :: [Double]))
$ DX.empty
Log.entries (Health.detectHeaderless df) `shouldBe` []
it "W002: コメント行 (# 始まり) を検出" $ do
let preview = "# header comment\n# more comment\nx,y\n1,2\n"
codes = map Log.lgCode (Log.entries (Health.detectCommentLines preview))
codes `shouldContain` ["W002"]
it "W005: 1 列 DataFrame + プレビューにタブ → delimiter ミスマッチ" $ do
let df = DX.insertColumn "x\ty" (DX.fromList ([1.0] :: [Double]))
$ DX.empty
preview = "x\ty\n1\t2\n3\t4\n"
codes = map Log.lgCode
(Log.entries (Health.detectDelimiterMismatch preview df))
codes `shouldContain` ["W005"]
it "W008: 通貨記号付き列を検出" $ do
let df = DX.insertColumn "price"
(DX.fromList (["$1,234.56", "$2,500.00", "$3,000.00", "$4,000"] :: [T.Text]))
$ DX.empty
codes = map Log.lgCode (Log.entries (Health.detectThousandsCurrency df))
codes `shouldContain` ["W008"]
-- BS インポートを使う何かのスモーク (未使用 warning 防止)
it "preview is non-empty for typical use" $
BS.length "x,y\n1,2" `shouldSatisfy` (> 0)
describe "Hanalyze.DataIO.CSV.loadAutoSafeWith" $ do
it "--no-header: 先頭行をデータ行として扱い col0... を生成" $
withSystemTempFile "ha-noh.csv" $ \fp h -> do
hPutStr h "1,2\n3,4\n5,6\n"
hClose h
r <- CSV.loadAutoSafeWith
(CSV.defaultLoadOpts { CSV.loNoHeader = True }) fp
case r of
Left e -> expectationFailure ("unexpected Left: " ++ e)
Right (df, lg) -> do
let cols = DX.columnNames df
cols `shouldBe` ["col0", "col1"]
map Log.lgCode (Log.entries lg) `shouldContain` ["I012"]
it "--skip 2: 先頭 2 行を skip" $
withSystemTempFile "ha-skip.csv" $ \fp h -> do
hPutStr h "# c1\n# c2\nx,y\n1,2\n3,4\n"
hClose h
r <- CSV.loadAutoSafeWith
(CSV.defaultLoadOpts { CSV.loSkip = 2 }) fp
case r of
Left e -> expectationFailure ("unexpected Left: " ++ e)
Right (df, _) -> DX.columnNames df `shouldBe` ["x", "y"]
it "sniff: ヘッダ無し CSV を自動推論で col0... に変える" $
withSystemTempFile "ha-sniff-noh.csv" $ \fp h -> do
hPutStr h "1.0,2.0\n3.0,4.0\n"
hClose h
r <- CSV.loadAutoSafeWith CSV.defaultLoadOpts fp
case r of
Left e -> expectationFailure ("unexpected Left: " ++ e)
Right (df, lg) -> do
DX.columnNames df `shouldBe` ["col0", "col1"]
map Log.lgCode (Log.entries lg) `shouldContain` ["I013"]
it "sniff: コメント行 # を skip 推論" $
withSystemTempFile "ha-sniff-skip.csv" $ \fp h -> do
hPutStr h "# comment 1\n# comment 2\nx,y\n1,2\n3,4\n"
hClose h
r <- CSV.loadAutoSafeWith CSV.defaultLoadOpts fp
case r of
Left e -> expectationFailure ("unexpected Left: " ++ e)
Right (df, _) -> DX.columnNames df `shouldBe` ["x", "y"]
it "sniff: セミコロン区切りを自動検出" $
withSystemTempFile "ha-sniff-semi.csv" $ \fp h -> do
hPutStr h "a;b;c\n1;2;3\n4;5;6\n"
hClose h
r <- CSV.loadAutoSafeWith CSV.defaultLoadOpts fp
case r of
Left e -> expectationFailure ("unexpected Left: " ++ e)
Right (df, _) -> DX.columnNames df `shouldBe` ["a", "b", "c"]
it "sniff: --no-sniff で自動推論を切れる" $
withSystemTempFile "ha-no-sniff.csv" $ \fp h -> do
hPutStr h "1.0,2.0\n3.0,4.0\n"
hClose h
r <- CSV.loadAutoSafeWith
(CSV.defaultLoadOpts { CSV.loSniff = False }) fp
case r of
Left e -> expectationFailure ("unexpected Left: " ++ e)
Right (df, lg) -> do
-- ヘッダ無しの自動修復は走らないので col0 にはならない
DX.columnNames df `shouldBe` ["1.0", "2.0"]
-- 代わりに W001 が出る
map Log.lgCode (Log.entries lg) `shouldContain` ["W001"]
it "Clean.stripUnitsCol: 12.3kg → 12.3" $ do
let df0 = DX.insertColumn "w"
(DX.fromList (["12.3kg", "11.5cm", "10kg"] :: [T.Text]))
$ DX.empty
(df1, lg) = Clean.applyRule Clean.StripUnits "w" df0
map Log.lgCode (Log.entries lg) `shouldContain` ["I100"]
case Conv2.getDoubleVec "w" df1 of
Just v -> V.toList v `shouldBe` [12.3, 11.5, 10.0]
Nothing -> expectationFailure "expected numeric column"
it "Clean.parseCurrencyCol: $1,234.56 → 1234.56" $ do
let df0 = DX.insertColumn "p"
(DX.fromList (["$1,234.56", "$2,500.00"] :: [T.Text]))
$ DX.empty
(df1, _) = Clean.applyRule Clean.ParseCurrency "p" df0
case Conv2.getDoubleVec "p" df1 of
Just v -> V.toList v `shouldBe` [1234.56, 2500.0]
Nothing -> expectationFailure "expected numeric column"
it "Clean.coerceNumericCol: 混在パターンを最大限拾う" $ do
let df0 = DX.insertColumn "x"
(DX.fromList (["12.3", "12.3kg", "$1,000"] :: [T.Text]))
$ DX.empty
(df1, _) = Clean.applyRule Clean.CoerceNumeric "x" df0
case Conv2.getDoubleVec "x" df1 of
Just v -> V.toList v `shouldBe` [12.3, 12.3, 1000.0]
Nothing -> expectationFailure "expected all-success column"
it "Preprocess.meltLonger: wide → long、NA セルは除外、列名を Double に parse" $ do
let df0 = DX.insertColumn "id" (DX.fromList (["a", "b"] :: [T.Text]))
$ DX.insertColumn "1" (DX.fromList ([Just 10.0, Nothing] :: [Maybe Double]))
$ DX.insertColumn "2" (DX.fromList ([Just 20.0, Just 30.0] :: [Maybe Double]))
$ DX.insertColumn "3" (DX.fromList ([Nothing, Just 60.0] :: [Maybe Double]))
$ DX.empty
df1 = Pp.meltLonger ["id"] ["1", "2", "3"] "t" "y" True df0
(nrows, ncols) = DX.dimensions df1
nrows `shouldBe` 4 -- a,1=10; a,2=20; b,2=30; b,3=60
ncols `shouldBe` 3 -- id, t, y
DX.columnNames df1 `shouldMatchList` ["id", "t", "y"]
case Conv2.getDoubleVec "y" df1 of
Just v -> sort (V.toList v) `shouldBe` [10, 20, 30, 60]
Nothing -> expectationFailure "expected y as numeric"
case Conv2.getDoubleVec "t" df1 of
Just v -> sort (V.toList v) `shouldBe` [1, 2, 2, 3]
Nothing -> expectationFailure "expected t parsed as numeric"
it "Preprocess.regridLong: ZIntersection モードで全 id が共通範囲に収まる" $ do
-- id=a: z=0..3, id=b: z=1..4 → intersection は (1, 3)
let df0 = DX.insertColumn "id" (DX.fromList (["a","a","a","a","b","b","b","b"] :: [T.Text]))
$ DX.insertColumn "z" (DX.fromList ([0,1,2,3,1,2,3,4] :: [Double]))
$ DX.insertColumn "y" (DX.fromList ([0,1,4,9,1,4,9,16] :: [Double]))
$ DX.empty
opts = Pp.defaultRegridOpts
{ Pp.roN = 5, Pp.roZBoundsMode = Pp.ZIntersection
, Pp.roGridKind = AG.Uniform }
rr = Pp.regridLong "id" "z" "y" opts df0
Pp.rrZMin rr `shouldBe` 1.0
Pp.rrZMax rr `shouldBe` 3.0
length (Pp.rrZGrid rr) `shouldBe` 5
length (Pp.rrIds rr) `shouldBe` 2
it "Preprocess.regridLong: ZUnion モードで [min,max] が和集合になる" $ do
let df0 = DX.insertColumn "id" (DX.fromList (["a","a","b","b"] :: [T.Text]))
$ DX.insertColumn "z" (DX.fromList ([0,2,1,3] :: [Double]))
$ DX.insertColumn "y" (DX.fromList ([0,4,1,9] :: [Double]))
$ DX.empty
opts = Pp.defaultRegridOpts
{ Pp.roN = 4, Pp.roZBoundsMode = Pp.ZUnion
, Pp.roGridKind = AG.Uniform }
rr = Pp.regridLong "id" "z" "y" opts df0
Pp.rrZMin rr `shouldBe` 0.0
Pp.rrZMax rr `shouldBe` 3.0
-- 外挿が記録される (id=a の上端 0..2、共通 0..3 → above=1)
let stat_a = head [s | s <- Pp.rrPerIdStats rr, Pp.piId s == "a"]
Pp.piExtrapAbove stat_a `shouldBe` 1.0
it "Clean.cleanPipeline: 複数列を一括変換" $ do
let df0 = DX.insertColumn "p"
(DX.fromList (["$10", "$20"] :: [T.Text]))
$ DX.insertColumn "w"
(DX.fromList (["1kg", "2kg"] :: [T.Text]))
$ DX.empty
rules = [("p", Clean.ParseCurrency), ("w", Clean.StripUnits)]
(df1, lg) = Clean.cleanPipeline rules df0
codes = map Log.lgCode (Log.entries lg)
codes `shouldContain` ["I101"]
codes `shouldContain` ["I100"]
Conv2.getDoubleVec "p" df1 `shouldSatisfy` \mv ->
case mv of { Just v -> V.toList v == [10, 20]; Nothing -> False }
it "--strict + 警告ありデータ (sniff off) → Left" $
withSystemTempFile "ha-strict.csv" $ \fp h -> do
hPutStr h "1.0,2.0\n3.0,4.0\n" -- ヘッダ無し疑い W001
hClose h
-- sniff を切ると W001 が残るので strict が短絡する
r <- CSV.loadAutoSafeWith
(CSV.defaultLoadOpts { CSV.loStrict = True
, CSV.loSniff = False }) fp
case r of
Left _ -> return ()
Right _ -> expectationFailure "expected Left under --strict --no-sniff"
-- ===========================================================================
-- 多出力 API の q=1 等価性 (M1〜M8)
-- ===========================================================================
describe "Multi-output equivalence (q=1)" $ do
let xs = LA.fromLists [[1,1.0], [1,2.0], [1,3.0], [1,4.0], [1,5.0]] :: LA.Matrix Double
yV = LA.fromList [2.1, 3.9, 6.0, 8.1, 10.0] :: LA.Vector Double
yM = LA.asColumn yV
approx tol a b = abs (a - b) < tol
approxList tol as bs = length as == length bs &&
all (uncurry (approx tol)) (zip as bs)
buildGroupsLocal gvec =
let lbls = V.fromList . sort . foldr (\x acc -> if x `elem` acc then acc else x:acc) [] $ V.toList gvec
qN = V.length lbls
idxFor x = case V.elemIndex x lbls of
Just i -> i
Nothing -> 0
idx = V.map idxFor gvec
sz = V.fromList [ V.length (V.filter (== j) idx) | j <- [0 .. qN - 1] ]
in (lbls, idx, sz)
it "M1 Regularized Ridge: fitRegularized == fitRegularizedMulti col 0" $ do
let single = Reg.fitRegularized (Reg.L2 0.1) xs yV
multi = Reg.fitRegularizedMulti (Reg.L2 0.1) xs yM
extr = Reg.regFitFromMulti 0 multi
approxList 1e-9 (LA.toList (Reg.rfBeta single))
(LA.toList (Reg.rfBeta extr))
`shouldBe` True
it "M1 Regularized Lasso: q=1 一致" $ do
let single = Reg.fitRegularized (Reg.L1 0.05) xs yV
multi = Reg.fitRegularizedMulti (Reg.L1 0.05) xs yM
extr = Reg.regFitFromMulti 0 multi
approxList 1e-9 (LA.toList (Reg.rfBeta single))
(LA.toList (Reg.rfBeta extr))
`shouldBe` True
it "M2 Spline: fitSpline == fitSplineMulti col 0" $ do
let xv = V.fromList [1,2,3,4,5,6,7,8,9,10] :: V.Vector Double
yv = V.fromList (map (\x -> sin (x/2) + 0.01*x) (V.toList xv))
knots = [1,3,5,7,10]
single = Sp.fitSpline (Sp.BSpline 3) knots xv yv
ymat = LA.asColumn (LA.fromList (V.toList yv))
multi = Sp.fitSplineMulti (Sp.BSpline 3) knots xv ymat
colS = LA.toList (Sp.sfBeta single)
colM = LA.toList (LA.flatten (Sp.smfBeta multi LA.¿ [0]))
approxList 1e-9 colS colM `shouldBe` True
it "M3 Kernel Ridge: kernelRidge == kernelRidgeMulti col 0" $ do
let xv = V.fromList [0.0,1,2,3,4,5,6,7,8,9] :: V.Vector Double
yv = V.fromList [0.0, 0.5, 1.0, 1.4, 1.7, 1.9, 2.0, 2.0, 1.95, 1.8]
single = K.kernelRidge K.Gaussian 1.0 0.01 xv yv
ymat = LA.asColumn (LA.fromList (V.toList yv))
multi = K.kernelRidgeMulti K.Gaussian 1.0 0.01 xv ymat
approxList 1e-9 (LA.toList (K.krAlpha single))
(LA.toList (LA.flatten (K.krmAlpha multi LA.¿ [0])))
`shouldBe` True
it "M3 Kernel NW: nwRegression == nwRegressionMulti col 0" $ do
let xv = V.fromList [0.0,1,2,3,4,5,6,7,8,9] :: V.Vector Double
yv = V.fromList [0.1, 0.3, 0.7, 1.0, 1.5, 1.9, 2.0, 1.95, 1.8, 1.5]
xn = V.fromList [0.5, 2.5, 5.5, 8.5]
single = K.nwRegression K.Gaussian 1.0 xv yv xn
ymat = LA.asColumn (LA.fromList (V.toList yv))
multi = K.nwRegressionMulti K.Gaussian 1.0 xv ymat xn
colM = LA.toList (LA.flatten (multi LA.¿ [0]))
approxList 1e-9 (V.toList single) colM `shouldBe` True
it "M4 RFF Ridge: rffRidge == rffRidgeMulti col 0" $ do
gen <- MWC.create
rff <- RFF.sampleRFFRBF 16 1.0 1.0 gen
let xList = [0.0, 1, 2, 3, 4, 5]
yList = [0.1, 0.5, 1.0, 1.4, 1.7, 1.9]
single = RFF.rffRidge rff xList yList 0.01
ymat = LA.asColumn (LA.fromList yList)
multi = RFF.rffRidgeMulti rff xList ymat 0.01
approxList 1e-9 (LA.toList (RFF.rffrWeights single))
(LA.toList (LA.flatten (RFF.rffrmWeights multi LA.¿ [0])))
`shouldBe` True
it "M5 GP: fitGP mean == fitGPMulti col 0" $ do
let model = GP.GPModel GP.RBF GP.defaultGPParams
trX = [0.0, 1, 2, 3, 4, 5]
trY = [0.1, 0.4, 0.9, 1.3, 1.6, 1.8]
tsX = [0.5, 2.5, 4.5]
single = GP.fitGP model trX trY tsX
ymat = LA.asColumn (LA.fromList trY)
(mMat, _) = GP.fitGPMulti model trX ymat tsX
approxList 1e-9 (GP.gpMean single)
(LA.toList (LA.flatten (mMat LA.¿ [0])))
`shouldBe` True
it "M5 GPRobust: fitGPRobust α == fitGPRobustMulti col 0" $ do
let params = GP.defaultGPParams
trX = [0.0, 1, 2, 3, 4, 5]
trY = [0.1, 0.4, 5.0, 1.3, 1.6, 1.8] -- 1 outlier at idx 2
single = GPR.fitGPRobust GP.RBF params (GPR.RStudentT 4 0.3) trX trY
ymat = LA.asColumn (LA.fromList trY)
multi = GPR.fitGPRobustMulti GP.RBF params (GPR.RStudentT 4 0.3) trX ymat
firstFit = head (GPR.rgmFits multi)
approxList 1e-9 (LA.toList (GPR.rgpAlpha single))
(LA.toList (GPR.rgpAlpha firstFit))
`shouldBe` True
it "M6 GLM Gaussian: fitGLM == fitGLMMulti col 0" $ do
let single = GLM.fitGLM GLM.Gaussian xs yV
multi = GLM.fitGLMMulti GLM.Gaussian GLM.Identity xs yM
colM = LA.toList (LA.flatten (GLM.gfmBeta multi LA.¿ [0]))
colS = LA.toList (LA.flatten (Core.coefficients single))
approxList 1e-7 colS colM `shouldBe` True
it "M7 LME: fitLME == fitLMEMulti col 0" $ do
let xMat = LA.fromLists
[[1,1],[1,2],[1,3],[1,4],
[1,1],[1,2],[1,3],[1,4],
[1,1],[1,2],[1,3],[1,4]] :: LA.Matrix Double
y1 = LA.fromList [7.1,6.9,7.0,7.0, 5.0,4.9,5.1,5.0, 3.0,2.9,3.1,3.0]
ym = LA.asColumn y1
gv = V.fromList (["A","A","A","A","B","B","B","B","C","C","C","C"] :: [T.Text])
(lbls, idx, sz) = buildGroupsLocal gv
single = fitLME xMat y1 idx lbls sz
multi = fitLMEMulti xMat ym idx lbls sz
firstM = head (glmmFits multi)
glmmRandVar firstM `shouldSatisfy` approx 1e-9 (glmmRandVar single)
-- ===========================================================================
-- 単目的オプティマイザ (Hanalyze.Optim.NelderMead)
-- ===========================================================================
describe "Hanalyze.Optim.NelderMead" $ do
let l2 :: [Double] -> [Double] -> Double
l2 a b = sqrt (sum (zipWith (\x y -> (x-y)^(2::Int)) a b))
sphere xs = sum [x*x | x <- xs]
rosenbrock [x, y] = (1 - x)^(2::Int) + 100 * (y - x*x)^(2::Int)
rosenbrock _ = error "rosenbrock: 2D only"
it "minimises sphere f(x)=Σx² to ~0 from x0=[3,-2,1]" $ do
r <- NM.runNelderMead sphere [3, -2, 1]
OC.orValue r `shouldSatisfy` (< 1e-6)
it "minimises Rosenbrock 2D to (1,1) within 0.05" $ do
let cfg = NM.defaultNMConfig
{ NM.nmStop = OC.defaultStopCriteria { OC.stMaxIter = 5000 } }
r <- NM.runNelderMeadWith cfg rosenbrock [-1.2, 1.0]
l2 (OC.orBest r) [1, 1] `shouldSatisfy` (< 0.05)
it "Maximize: -sphere has optimum 0 at origin" $ do
let cfg = NM.defaultNMConfig { NM.nmDir = OC.Maximize }
r <- NM.runNelderMeadWith cfg (\xs -> negate (sphere xs)) [3, -2, 1]
-- Maximize: orValue は元尺度 (= negate sphere の最大値、すなわち 0 に近い)
OC.orValue r `shouldSatisfy` (\v -> v > -1e-6)
-- 最良点は原点近傍
l2 (OC.orBest r) [0, 0, 0] `shouldSatisfy` (< 1e-2)
-- ===========================================================================
-- 単目的オプティマイザ (Hanalyze.Optim.LBFGS)
-- ===========================================================================
describe "Hanalyze.Optim.LBFGS" $ do
let l2 :: [Double] -> [Double] -> Double
l2 a b = sqrt (sum (zipWith (\x y -> (x-y)^(2::Int)) a b))
sphere xs = sum [x*x | x <- xs]
sphereGrad xs = [2*x | x <- xs]
rosen [x, y] = (1-x)^(2::Int) + 100*(y - x*x)^(2::Int)
rosen _ = error "rosen: 2D"
rosenGrad [x, y] =
[ -2*(1-x) - 400*x*(y - x*x), 200*(y - x*x) ]
rosenGrad _ = error "rosenGrad: 2D"
it "minimises sphere 5D with analytic grad to ~0" $ do
r <- LBFGS.runLBFGS sphere sphereGrad [3, -2, 1, 0.5, -1.5]
OC.orValue r `shouldSatisfy` (< 1e-8)
it "minimises Rosenbrock 2D within 0.01 of (1,1)" $ do
let cfg = LBFGS.defaultLBFGSConfig
{ LBFGS.lbStop = OC.defaultStopCriteria { OC.stMaxIter = 500 } }
r <- LBFGS.runLBFGSWith cfg rosen rosenGrad [-1.2, 1.0]
l2 (OC.orBest r) [1, 1] `shouldSatisfy` (< 0.01)
it "numeric gradient: sphere 30D converges" $ do
let x0 = take 30 (cycle [1.5, -2.0, 0.5])
r <- LBFGS.runLBFGSNumeric LBFGS.defaultLBFGSConfig sphere x0
OC.orValue r `shouldSatisfy` (< 1e-4)
-- ===========================================================================
-- 1D オプティマイザ (Hanalyze.Optim.LineSearch)
-- ===========================================================================
describe "Hanalyze.Optim.LineSearch" $ do
let parabola [x] = (x - 2.5)^(2::Int) + 1.0
parabola _ = error "1D"
cosBowl [x] = cos x + 0.1 * x * x -- 単峰、最小 ≈ 1.428 付近
cosBowl _ = error "1D"
it "Brent: parabola minimum at x = 2.5" $ do
let r = LS.brent LS.defaultBrentConfig parabola 0 5
abs (head (OC.orBest r) - 2.5) `shouldSatisfy` (< 1e-5)
OC.orValue r `shouldSatisfy` (\v -> abs (v - 1) < 1e-8)
it "Brent: cos x + 0.1 x² minimum (verified by GS)" $ do
let rB = LS.brent LS.defaultBrentConfig cosBowl 0 4
rG = LS.goldenSection OC.Minimize cosBowl 0 4 1e-8 200
abs (head (OC.orBest rB) - head (OC.orBest rG)) `shouldSatisfy` (< 1e-3)
it "GoldenSection: parabola minimum at x = 2.5" $ do
let r = LS.goldenSection OC.Minimize parabola 0 5 1e-7 200
abs (head (OC.orBest r) - 2.5) `shouldSatisfy` (< 1e-3)
it "Kernel.autoBandwidthBrent: 同じ最適 h を grid 法とほぼ一致" $ do
let xs = V.fromList [0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
ys = V.map (\x -> sin x + 0.1 * x) xs
(hG, _) = K.gridSearchBandwidth K.Gaussian xs ys [0.5, 0.8, 1.0, 1.5, 2.0, 3.0]
(hB, _) = K.autoBandwidthBrent K.Gaussian xs ys 0.3 4.0
abs (hB - hG) `shouldSatisfy` (< 1.0) -- グリッドと近い領域
-- ===========================================================================
-- 大域オプティマイザ (Hanalyze.Optim.DifferentialEvolution)
-- ===========================================================================
describe "Hanalyze.Optim.DifferentialEvolution" $ do
let sphere xs = sum [x*x | x <- xs]
rastrigin xs =
10 * fromIntegral (length xs) +
sum [x*x - 10 * cos (2 * pi * x) | x <- xs]
l2 a b = sqrt (sum (zipWith (\x y -> (x-y)^(2::Int)) a b))
it "DE: sphere 5D が原点付近に到達" $ do
gen <- MWC.create
let bs = replicate 5 (-5, 5)
cfg = (DE.defaultDEConfig bs)
{ DE.deStop = OC.defaultStopCriteria { OC.stMaxIter = 200 } }
r <- DE.runDEWith cfg sphere gen
OC.orValue r `shouldSatisfy` (< 1e-3)
it "DE: Rastrigin 3D の大域最小 (原点) を見つける" $ do
gen <- MWC.create
let bs = replicate 3 (-5.12, 5.12)
cfg = (DE.defaultDEConfig bs)
{ DE.deStop = OC.defaultStopCriteria { OC.stMaxIter = 400 } }
r <- DE.runDEWith cfg rastrigin gen
l2 (OC.orBest r) [0, 0, 0] `shouldSatisfy` (< 0.5)
-- ===========================================================================
-- 大域オプティマイザ (Hanalyze.Optim.CMAES、簡易対角版)
-- ===========================================================================
describe "Hanalyze.Optim.CMAES" $ do
let sphere xs = sum [x*x | x <- xs]
l2 a b = sqrt (sum (zipWith (\x y -> (x-y)^(2::Int)) a b))
it "CMA-ES: sphere 5D を最小化、原点近傍に到達" $ do
gen <- MWC.create
let cfg = CMAES.defaultCMAESConfig
{ CMAES.cmStop = OC.defaultStopCriteria { OC.stMaxIter = 200 }
, CMAES.cmSigma0 = 1.0 }
r <- CMAES.runCMAESWith cfg sphere [3, -2, 1, 0.5, -1.5] gen
l2 (OC.orBest r) [0,0,0,0,0] `shouldSatisfy` (< 0.5)
it "CMA-ES Full: sphere 5D で 1e-3 以内に到達" $ do
gen <- MWC.create
let cfg = CMAESF.defaultCMAESFConfig
{ CMAESF.cmfStop = OC.defaultStopCriteria { OC.stMaxIter = 300 }
, CMAESF.cmfSigma0 = 1.0 }
r <- CMAESF.runCMAESFullWith cfg sphere [3, -2, 1, 0.5, -1.5] gen
OC.orValue r `shouldSatisfy` (< 1e-3)
it "CMA-ES Full: Rosenbrock 2D で (1,1) に 0.1 以内" $ do
gen <- MWC.create
let cfg = CMAESF.defaultCMAESFConfig
{ CMAESF.cmfStop = OC.defaultStopCriteria { OC.stMaxIter = 500 }
, CMAESF.cmfSigma0 = 0.5
, CMAESF.cmfLambda = Just 20 }
rosen [x, y] = (1-x)^(2::Int) + 100 * (y - x*x)^(2::Int)
rosen _ = error "2D"
r <- CMAESF.runCMAESFullWith cfg rosen [-1.2, 1.0] gen
l2 (OC.orBest r) [1, 1] `shouldSatisfy` (< 0.1)
-- ===========================================================================
-- メタヒューリスティック (Tier 2: Simulated Annealing, PSO)
-- ===========================================================================
describe "Hanalyze.Optim.SimulatedAnnealing" $ do
it "SA: sphere 5D で sufficient annealing" $ do
gen <- MWC.create
let bs = replicate 5 (-3, 3)
cfg = (SA.defaultSAConfig bs)
{ SA.saStop = OC.defaultStopCriteria { OC.stMaxIter = 5000 }
, SA.saInitTemp = 2.0
, SA.saSchedule = SA.Geometric 0.997 }
sphere xs = sum [x*x | x <- xs]
r <- SA.runSAWith cfg sphere [2, -1.5, 1, 0.5, -0.7] gen
OC.orValue r `shouldSatisfy` (< 0.5)
describe "Hanalyze.Optim.ParticleSwarm" $ do
it "PSO: sphere 5D で原点近傍 (0.5 以内)" $ do
gen <- MWC.create
let bs = replicate 5 (-5, 5)
cfg = (PSO.defaultPSOConfig bs)
{ PSO.psoStop = OC.defaultStopCriteria { OC.stMaxIter = 200 }
, PSO.psoNum = 30 }
sphere xs = sum [x*x | x <- xs]
r <- PSO.runPSOWith cfg sphere gen
OC.orValue r `shouldSatisfy` (< 0.5)
it "PSO: Rastrigin 3D の大域最小に近い" $ do
gen <- MWC.create
let bs = replicate 3 (-5.12, 5.12)
cfg = (PSO.defaultPSOConfig bs)
{ PSO.psoStop = OC.defaultStopCriteria { OC.stMaxIter = 300 }
, PSO.psoNum = 40 }
rastrigin xs =
10 * fromIntegral (length xs) +
sum [x*x - 10 * cos (2 * pi * x) | x <- xs]
r <- PSO.runPSOWith cfg rastrigin gen
OC.orValue r `shouldSatisfy` (< 5.0) -- 大域近傍 (10 程度の局所有り)
-- ===========================================================================
-- 制約付き最適化 (Augmented Lagrangian)
-- ===========================================================================
describe "Hanalyze.Optim.Constrained" $ do
it "Augmented Lagrangian: min x1²+x2² s.t. x1+x2=1 → (0.5, 0.5)" $ do
let f xs = (head xs)^(2::Int) + (xs !! 1)^(2::Int)
cs = Con.ConstraintSet
{ Con.csEq = [\xs -> head xs + xs !! 1 - 1]
, Con.csIneq = []
}
(r, viol) <- Con.runAugmentedLagrangian
Con.defaultConstrainedConfig f cs [0, 0]
let [x1, x2] = OC.orBest r
abs (x1 - 0.5) `shouldSatisfy` (< 0.05)
abs (x2 - 0.5) `shouldSatisfy` (< 0.05)
viol `shouldSatisfy` (< 1e-3)
it "Augmented Lagrangian: 不等式 x ≥ 1 (h(x)=1-x≤0) で min x²" $ do
let f xs = (head xs)^(2::Int)
cs = Con.ConstraintSet
{ Con.csEq = []
, Con.csIneq = [\xs -> 1 - head xs] -- 1 - x ≤ 0 ⇔ x ≥ 1
}
(r, viol) <- Con.runAugmentedLagrangian
Con.defaultConstrainedConfig f cs [0]
let [x] = OC.orBest r
x `shouldSatisfy` (\v -> v >= 1 - 0.01)
viol `shouldSatisfy` (< 1e-3)
it "Penalty method: 等式 x1+x2=1 (簡易版)" $ do
let f xs = (head xs)^(2::Int) + (xs !! 1)^(2::Int)
cs = Con.ConstraintSet
{ Con.csEq = [\xs -> head xs + xs !! 1 - 1]
, Con.csIneq = []
}
(r, _) <- Con.penaltyMethod Con.defaultConstrainedConfig f cs [0, 0]
let [x1, x2] = OC.orBest r
abs (x1 + x2 - 1) `shouldSatisfy` (< 0.05)
it "boxToIneq: bounds [(1,5)] を不等式 2 本に展開、x=3 で両方満たす" $ do
let ineqs = Con.boxToIneq [(1.0, 5.0)]
length ineqs `shouldBe` 2
all (\h -> h [3.0] <= 0) ineqs `shouldBe` True
any (\h -> h [0.0] > 0) ineqs `shouldBe` True -- 下限違反
any (\h -> h [6.0] > 0) ineqs `shouldBe` True -- 上限違反
-- ===========================================================================
-- box 制約 (Bounds) 統一インターフェース (Hanalyze.Optim.Common)
-- ===========================================================================
describe "Hanalyze.Optim.Common box constraints" $ do
it "clipToBounds: 範囲外を反射、範囲内はそのまま" $ do
OC.clipToBounds [(0,10)] [5] `shouldBe` [5]
OC.clipToBounds [(0,10)] [-2] `shouldBe` [2] -- 反射
OC.clipToBounds [(0,10)] [12] `shouldBe` [8]
it "boundsPenalty: 範囲内 0、範囲外で大きい正値" $ do
OC.boundsPenalty (Just [(0,10)]) [5] `shouldBe` 0
OC.boundsPenalty Nothing [-99] `shouldBe` 0
OC.boundsPenalty (Just [(0,10)]) [-1] `shouldSatisfy` (> 1e5)
it "NelderMead: nmBounds で box 内に解が収まる" $ do
let f xs = (head xs - 5)^(2::Int) -- 真の最小は x=5
cfg = NM.defaultNMConfig { NM.nmBounds = Just [(0, 2)] }
r <- NM.runNelderMeadWith cfg f [1.0]
head (OC.orBest r) `shouldSatisfy` (\v -> v >= -0.1 && v <= 2.1)
it "LBFGS: lbBounds で box 内に解が収まる" $ do
let f xs = (head xs - 5)^(2::Int)
cfg = LBFGS.defaultLBFGSConfig { LBFGS.lbBounds = Just [(0, 2)] }
r <- LBFGS.runLBFGSNumeric cfg f [1.0]
head (OC.orBest r) `shouldSatisfy` (\v -> v >= -0.1 && v <= 2.1)
it "CMAES (簡易版): cmBounds で box 内サンプル" $ do
gen <- MWC.create
let f xs = sum [x*x | x <- xs]
cfg = CMAES.defaultCMAESConfig
{ CMAES.cmBounds = Just [(-1, 1), (-1, 1)]
, CMAES.cmStop = (CMAES.cmStop CMAES.defaultCMAESConfig)
{ OC.stMaxIter = 80 }
}
r <- CMAES.runCMAESWith cfg f [0.5, 0.5] gen
all (\v -> v >= -1.05 && v <= 1.05) (OC.orBest r) `shouldBe` True
it "CMAESFull: cmfBounds で box 内サンプル" $ do
gen <- MWC.create
let f xs = sum [x*x | x <- xs]
cfg = CMAESF.defaultCMAESFConfig
{ CMAESF.cmfBounds = Just [(-1, 1), (-1, 1)]
, CMAESF.cmfStop = (CMAESF.cmfStop CMAESF.defaultCMAESFConfig)
{ OC.stMaxIter = 80 }
}
r <- CMAESF.runCMAESFullWith cfg f [0.5, 0.5] gen
all (\v -> v >= -1.05 && v <= 1.05) (OC.orBest r) `shouldBe` True
-- ===========================================================================
-- Bayesian Optimization 内部最適化の差し替え (Hanalyze.Optim.BayesOpt)
-- ===========================================================================
describe "Hanalyze.Optim.BayesOpt (acquisition optimizer swap)" $ do
it "bayesOpt 1D: Brent 内側で簡単な凸関数 (x-1.5)^2 を見つける" $ do
gen <- MWC.create
let cfg = BO.defaultBayesOptConfig
{ BO.boIterations = 8
, BO.boInitPoints = 4
, BO.boGridSize = 32
}
target x = pure ((x - 1.5)^(2::Int) :: Double)
(_, (xb, _)) <- BO.bayesOpt cfg target (0, 3) gen
-- 8 反復では精度は緩めに。3.0 範囲のうち 0.5 以内に収束を期待
abs (xb - 1.5) `shouldSatisfy` (< 0.5)
-- ===========================================================================
-- RFF HP 自動チューニングの DE 版 (Phase O9)
-- ===========================================================================
describe "Hanalyze.Model.RFF DE-based auto-HP" $ do
it "maximizeMarginalLikRBFMV_DE: y = sin x + noise で妥当な ℓ" $ do
gen <- MWC.create
let n = 30
xs = [ fromIntegral i / 5 | i <- [0 .. n - 1] ] :: [Double]
ys = [ sin x + 0.05 * cos (3 * x) | x <- xs ]
xMat = LA.fromColumns [LA.fromList xs]
yVec = LA.fromList ys
r <- RFF.maximizeMarginalLikRBFMV_DE xMat yVec 30 gen
-- ℓ が極端に小さくない (>1e-2) ことだけ確認
RFF.mlEll r `shouldSatisfy` (> 1e-2)
it "gridSearchLOOCVRBFMV_DE: LOOCV が有限値、ℓ が探索範囲内" $ do
gen <- MWC.create
let n = 25
xs = [ fromIntegral i / 4 | i <- [0 .. n - 1] ] :: [Double]
ys = [ x + 0.1 * sin (2 * x) | x <- xs ]
xMat = LA.fromColumns [LA.fromList xs]
yVec = LA.fromList ys
r <- RFF.gridSearchLOOCVRBFMV_DE 1 50 xMat yVec 20 gen
RFF.lcLOOCV r `shouldSatisfy` (\v -> not (isNaN v) && v >= 0)
RFF.lcEll r `shouldSatisfy` (> 1e-3)
-- ===========================================================================
-- Hanalyze.Viz.ReportBuilder.secInterpolation (Phase G4)
-- ===========================================================================
describe "Hanalyze.Viz.ReportBuilder.secInterpolation" $ do
it "defaultInterpReport で renderReport まで通る (smoke)" $
withSystemTempFile "ha-interp.html" $ \fp h -> do
hClose h
let ir = RB.defaultInterpReport "test"
RB.renderReport fp (RB.defaultReportConfig "T") [RB.secInterpolation ir]
out <- readFile fp
length out `shouldSatisfy` (> 100)
it "InterpReport 全フィールド + extra で HTML に主要要素が含まれる" $
withSystemTempFile "ha-interp2.html" $ \fp h -> do
hClose h
let ir = (RB.defaultInterpReport "regrid")
{ RB.irInterpKind = "PCHIP"
, RB.irGridKind = "Adaptive"
, RB.irN = 3
, RB.irPerIdObserved = [("a", [(0, 0), (1, 1)])]
, RB.irPerIdInterpY = [("a", [(0, 0), (0.5, 0.5), (1, 1)])]
, RB.irGrid = [0, 0.5, 1]
, RB.irDensity = [(0, 1), (0.5, 2), (1, 1)]
, RB.irPerIdSummary = [("a", 2, 0, 1, 0, 0, 0)]
, RB.irExtraEnabled = True
, RB.irPerIdYRange = [("a", 0, 1, 0, 1)]
}
RB.renderReport fp (RB.defaultReportConfig "T") [RB.secInterpolation ir]
out <- readFile fp
out `shouldContain` "Parameters"
out `shouldContain` "PCHIP"
-- ===========================================================================
-- Hanalyze.Stat.Test (Phase 1: hypothesis tests)
-- ===========================================================================
describe "Hanalyze.Stat.Test" $ do
it "tTest1Sample: μ₀=0 で同分布なら p > 0.05" $ do
let xs = LA.fromList [0.1, -0.2, 0.3, 0.0, 0.15, -0.1, 0.05, 0.2]
tr = ST.tTest1Sample xs 0 ST.TwoSided
ST.trPValue tr `shouldSatisfy` (> 0.05)
it "tTestWelch: 明らかにずれた 2 群で p < 0.05" $ do
let xs = LA.fromList [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
ys = LA.fromList [10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0]
tr = ST.tTestWelch xs ys ST.TwoSided
ST.trPValue tr `shouldSatisfy` (< 0.05)
it "anovaOneWay: 3 群が異なる平均なら有意" $ do
let g1 = LA.fromList [1.0, 2.0, 3.0]
g2 = LA.fromList [4.0, 5.0, 6.0]
g3 = LA.fromList [7.0, 8.0, 9.0]
tr = ST.anovaOneWay [g1, g2, g3]
ST.trStatistic tr `shouldSatisfy` (> 1)
ST.trPValue tr `shouldSatisfy` (< 0.05)
it "chiSquareIndep: 強い独立で p > 0.05" $ do
let tbl = LA.matrix 2 [25, 25, 25, 25] -- 完全独立
tr = ST.chiSquareIndep tbl
ST.trPValue tr `shouldSatisfy` (> 0.5)
it "chiSquareIndep: 強い従属で p < 0.05" $ do
let tbl = LA.matrix 2 [40, 10, 10, 40] -- 高 chi2
tr = ST.chiSquareIndep tbl
ST.trPValue tr `shouldSatisfy` (< 0.05)
it "leveneTest: 分散が大きく異なれば p < 0.05" $ do
let g1 = LA.fromList [1, 1.1, 0.9, 1.05, 0.95, 1.02, 0.98, 1.03]
g2 = LA.fromList [10, 20, 5, 25, 8, 30, 3, 22] -- much larger var
tr = ST.leveneTest [g1, g2]
ST.trPValue tr `shouldSatisfy` (< 0.05)
it "shapiroWilk: roughly-normal 系列で p > 0.05" $ do
let xs = LA.fromList [-1.5, -0.5, 0.0, 0.3, 0.8, 1.2, -0.3, 0.5, 1.0, -1.0]
tr = ST.shapiroWilk xs
ST.trPValue tr `shouldSatisfy` (> 0.05)
it "mannWhitneyU: 明らかにずれた 2 群で p < 0.05" $ do
let xs = LA.fromList [1, 2, 3, 4, 5, 6]
ys = LA.fromList [11, 12, 13, 14, 15, 16]
tr = ST.mannWhitneyU xs ys ST.TwoSided
ST.trPValue tr `shouldSatisfy` (< 0.05)
it "fisherExact2x2: 強い偏りで p < 0.05" $ do
let tr = ST.fisherExact2x2 ((20, 5), (5, 20)) ST.TwoSided
ST.trPValue tr `shouldSatisfy` (< 0.05)
-- ===========================================================================
-- Hanalyze.Model.PCA (Phase 2)
-- ===========================================================================
describe "Hanalyze.Model.PCA" $ do
it "PCA on rank-1 matrix: 1st component explains ~100% var" $ do
let -- Rank-1: each row = scalar × [1, 2, 3]
ks = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
xs = LA.fromLists [[k, 2 * k, 3 * k] | k <- ks]
r = PCA.pca PCA.Center Nothing xs
head (LA.toList (PCA.pcaExplainedRatio r))
`shouldSatisfy` (> 0.999)
it "pcaTransform + pcaInverse は Center mode で全成分なら復元 ≈ x" $ do
let xs = LA.fromLists [[1, 2, 3], [4, 5, 6], [7, 8, 9], [2, 1, 0]]
r = PCA.pca PCA.Center Nothing xs
scores = PCA.pcaTransform r xs
recon = PCA.pcaInverse r scores
diff = LA.norm_2 (LA.flatten (xs - recon))
diff `shouldSatisfy` (< 1e-9)
it "pcaCumExplained は monotone increasing で max ≤ 1" $ do
let xs = LA.fromLists [[k, 2*k+1, k*k]
| k <- [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]]
r = PCA.pca PCA.Center Nothing xs
cum = LA.toList (PCA.pcaCumExplained r)
last cum `shouldSatisfy` (<= 1.0001)
and (zipWith (<=) cum (tail cum)) `shouldBe` True
it "CenterScale で各列の SD が 1 に正規化される" $ do
let xs = LA.fromLists [[1, 100, 0.1], [2, 200, 0.2],
[3, 300, 0.3], [4, 400, 0.4]]
r = PCA.pca PCA.CenterScale Nothing xs
-- SD は元データの per-col SD で保存される
LA.size (PCA.pcaScale r) `shouldBe` 3
all (> 0) (LA.toList (PCA.pcaScale r)) `shouldBe` True
-- ===========================================================================
-- Hanalyze.Stat.ClassMetrics (Phase 3)
-- ===========================================================================
describe "Hanalyze.Stat.ClassMetrics" $ do
it "confusionMatrix: 完全一致で TP/TN のみ" $ do
let c = CM.confusionMatrix [1, 0, 1, 0, 1] [1, 0, 1, 0, 1]
CM.confTP c `shouldBe` 3
CM.confTN c `shouldBe` 2
CM.confFP c `shouldBe` 0
CM.confFN c `shouldBe` 0
CM.accuracy c `shouldBe` 1.0
CM.f1Score c `shouldBe` 1.0
it "precision/recall/f1: 不均衡な誤分類" $ do
let c = CM.confusionMatrix [1, 1, 1, 0, 0] [1, 1, 0, 1, 0]
-- TP=2, FN=1, FP=1, TN=1
CM.precision c `shouldBe` 2/3 -- 2/(2+1)
CM.recall c `shouldBe` 2/3 -- 2/(2+1)
CM.f1Score c `shouldBe` 2/3
it "AUC: 完全分離で 1.0、ランダムスコアで ~0.5" $ do
let ys = [0, 0, 0, 1, 1, 1]
perfectScores = [0.1, 0.2, 0.3, 0.7, 0.8, 0.9]
aucPerfect = CM.auc ys perfectScores
aucPerfect `shouldBe` 1.0
it "logLoss: 自信ある正解で小、自信ある誤りで大" $ do
let lossGood = CM.logLoss [1, 0, 1, 0] [0.99, 0.01, 0.99, 0.01]
lossBad = CM.logLoss [1, 0, 1, 0] [0.01, 0.99, 0.01, 0.99]
lossGood `shouldSatisfy` (< 0.05)
lossBad `shouldSatisfy` (> 4.0)
it "brierScore: 完全予測で 0、最悪予測で 1" $ do
let bsGood = CM.brierScore [1, 0, 1, 0] [1.0, 0.0, 1.0, 0.0]
bsBad = CM.brierScore [1, 0, 1, 0] [0.0, 1.0, 0.0, 1.0]
bsGood `shouldSatisfy` (< 1e-10)
bsBad `shouldBe` 1.0
it "MCC: 完全一致で 1、ランダムで 0 付近" $ do
let cGood = CM.confusionMatrix [1, 0, 1, 0, 1] [1, 0, 1, 0, 1]
CM.matthewsCorr cGood `shouldBe` 1.0
it "macroF1 (multi-class): 完全分類で 1.0" $ do
let cm = CM.confusionMulti [0, 1, 2, 0, 1, 2] [0, 1, 2, 0, 1, 2]
CM.accuracyMulti cm `shouldBe` 1.0
CM.macroF1 cm `shouldBe` 1.0
-- ===========================================================================
-- Hanalyze.Stat.CV (Phase 4)
-- ===========================================================================
describe "Hanalyze.Stat.CV" $ do
it "kFold(5, 100): 5 fold で全 100 行を test に使用、重複なし" $ do
gen <- MWC.createSystemRandom
folds <- CV.kFold 5 100 gen
length folds `shouldBe` 5
let allTest = concatMap snd folds
length allTest `shouldBe` 100
length (V.toList (V.fromList allTest)) `shouldBe` 100 -- 重複なし
it "kFold: train + test = total samples per fold" $ do
gen <- MWC.createSystemRandom
folds <- CV.kFold 5 100 gen
mapM_ (\(tr, te) -> length tr + length te `shouldBe` 100) folds
it "leaveOneOut(10): 10 folds、test set size 1 each" $ do
folds <- CV.leaveOneOut 10
length folds `shouldBe` 10
mapM_ (\(_, te) -> length te `shouldBe` 1) folds
it "stratifiedKFold(3): クラスバランスがほぼ保持される" $ do
gen <- MWC.createSystemRandom
let labels = replicate 30 0 ++ replicate 30 1 ++ replicate 30 2
folds <- CV.stratifiedKFold 3 labels gen
length folds `shouldBe` 3
-- 各 fold の test set には各クラスの ~10 が含まれる
mapM_ (\(_, te) -> length te `shouldSatisfy` (\n -> n >= 27 && n <= 33)) folds
it "shuffleSplit: 反復回数とテストサイズが正しい" $ do
gen <- MWC.createSystemRandom
folds <- CV.shuffleSplit 5 0.2 100 gen
length folds `shouldBe` 5
mapM_ (\(_, te) -> length te `shouldBe` 20) folds
it "timeSeriesSplit: forward-chaining で過去のみで学習" $ do
let folds = CV.timeSeriesSplit 50 10 100 -- initial=50, step=10, n=100
length folds `shouldBe` 5 -- (100-50)/10 = 5 folds
-- 全 fold で train indices < min(test indices)
mapM_ (\(tr, te) ->
(maximum tr < minimum te) `shouldBe` True) folds
-- ===========================================================================
-- Hanalyze.Model.Cluster (Phase 5)
-- ===========================================================================
describe "Hanalyze.Model.Cluster" $ do
it "kMeans: 2 つの離れたクラスタで正しく分類" $ do
gen <- MWC.createSystemRandom
-- Cluster 1: around (0, 0); cluster 2: around (10, 10)
let xs = LA.fromLists $
[[0.1*x, 0.1*y] | x <- [-3..3], y <- [-3..3]] ++
[[10 + 0.1*x, 10 + 0.1*y] | x <- [-3..3], y <- [-3..3]]
cfg = Cl.defaultKMeansConfig 2
r <- Cl.kMeans cfg xs gen
LA.rows (Cl.kmrCentroids r) `shouldBe` 2
-- 全 49 points がクラスタ 1、49 が クラスタ 2
let labels = Cl.kmrLabels r
(c0, c1) = (length (filter (== 0) labels), length (filter (== 1) labels))
(min c0 c1) `shouldBe` 49
(max c0 c1) `shouldBe` 49
it "silhouette: well-separated clusters で > 0.5" $ do
gen <- MWC.createSystemRandom
let xs = LA.fromLists $
[[0.1*x, 0.1*y] | x <- [-3..3], y <- [-3..3]] ++
[[20 + 0.1*x, 20 + 0.1*y] | x <- [-3..3], y <- [-3..3]]
cfg = Cl.defaultKMeansConfig 2
r <- Cl.kMeans cfg xs gen
let s = Cl.silhouette xs (Cl.kmrLabels r)
s `shouldSatisfy` (> 0.7)
it "kMeans: inertia は monotone non-increasing in iter" $ do
gen <- MWC.createSystemRandom
let xs = LA.fromLists [[fromIntegral i, fromIntegral j]
| i <- [0..9::Int], j <- [0..9::Int]]
cfg = (Cl.defaultKMeansConfig 4) { Cl.kmRestarts = 5 }
r <- Cl.kMeans cfg xs gen
Cl.kmrInertia r `shouldSatisfy` (>= 0)
Cl.kmrConverged r `shouldBe` True
-- ===========================================================================
-- Hanalyze.Stat.MultipleTesting (Phase 6)
-- ===========================================================================
describe "Hanalyze.Stat.MultipleTesting" $ do
it "Bonferroni: p × m, capped at 1" $ do
let ps = [0.01, 0.04, 0.05, 0.10]
adj = MT.bonferroni ps
adj `shouldBe` [0.04, 0.16, 0.20, 0.40]
it "Holm: 単調 + Bonferroni より緩い" $ do
let ps = [0.01, 0.02, 0.03, 0.04]
adj = MT.holm ps
-- 各 adj ≥ 元 p、最初は p × m = 0.04
head adj `shouldBe` 0.04
and (zipWith (<=) ps adj) `shouldBe` True
it "BH: monotonic non-decreasing in sorted p order" $ do
let ps = [0.01, 0.04, 0.03, 0.05]
adj = MT.benjaminiHochberg ps
length adj `shouldBe` 4
all (<= 1.0) adj `shouldBe` True
all (>= 0.0) adj `shouldBe` True
it "BY: BH より conservative" $ do
let ps = [0.01, 0.02, 0.03, 0.04, 0.05]
bh = MT.benjaminiHochberg ps
by = MT.benjaminiYekutieli ps
-- BY は cumulative harmonic factor を掛けるので大きい
and (zipWith (>=) by bh) `shouldBe` True
-- ===========================================================================
-- Hanalyze.Stat.Bootstrap (Phase 7)
-- ===========================================================================
describe "Hanalyze.Stat.Bootstrap" $ do
it "bootstrapCI on N(0,1) sample: 0 が 95% CI 内" $ do
gen <- MWC.createSystemRandom
let xs = LA.fromList [-1.0, -0.5, 0.0, 0.5, 1.0,
-0.3, 0.3, -0.7, 0.7, 0.0]
(lo, hi) <- Boot.bootstrapCI 2000 0.95 Boot.sampleMean xs gen
lo `shouldSatisfy` (< 0.5)
hi `shouldSatisfy` (> -0.5)
it "permutationTest: 異なる平均で p < 0.05" $ do
gen <- MWC.createSystemRandom
let xs = LA.fromList [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
ys = LA.fromList [10.0, 11.0, 12.0, 13.0, 14.0, 15.0]
(_diff, p) <- Boot.permutationTest 2000 xs ys gen
p `shouldSatisfy` (< 0.05)
it "sampleMean / sampleVar / sampleMedian の整合性" $ do
let v = LA.fromList [1.0, 2.0, 3.0, 4.0, 5.0]
Boot.sampleMean v `shouldBe` 3.0
Boot.sampleVar v `shouldBe` 2.5 -- variance of 1..5 (unbiased)
Boot.sampleMedian v `shouldBe` 3.0
-- ===========================================================================
-- Hanalyze.DataIO.Reshape (Phase 8)
-- ===========================================================================
describe "Hanalyze.DataIO.Reshape" $ do
it "lagColumn(1): 先頭が NaN、残りは 1 つずれる" $ do
let df = DX.fromNamedColumns
[("x", DX.fromList [1.0, 2.0, 3.0, 4.0, 5.0 :: Double])]
df' = Reshape.lagColumn 1 "x" "x_lag1" df
-- Lagged column should exist
"x_lag1" `elem` DX.columnNames df' `shouldBe` True
it "leadColumn(1): 末尾が NaN、残りは 1 つ前進" $ do
let df = DX.fromNamedColumns
[("x", DX.fromList [1.0, 2.0, 3.0, 4.0, 5.0 :: Double])]
df' = Reshape.leadColumn 1 "x" "x_lead1" df
"x_lead1" `elem` DX.columnNames df' `shouldBe` True
it "rollingMean(3): 最初 2 つが NaN、3 番目以降は窓内平均" $ do
let df = DX.fromNamedColumns
[("x", DX.fromList [1.0, 2.0, 3.0, 4.0, 5.0 :: Double])]
df' = Reshape.rollingMean 3 "x" "x_rmean3" df
"x_rmean3" `elem` DX.columnNames df' `shouldBe` True
it "oneHot: text 列を indicator 列に展開" $ do
let df = DX.fromNamedColumns
[ ("id", DX.fromList [1, 2, 3, 4, 5 :: Int])
, ("category", DX.fromList ["A", "B", "A", "C", "B" :: T.Text])
]
df' = Reshape.oneHot False "category" df
let cols = DX.columnNames df'
"category" `elem` cols `shouldBe` False -- 元列削除
"category_A" `elem` cols `shouldBe` True
"category_B" `elem` cols `shouldBe` True
"category_C" `elem` cols `shouldBe` True
it "oneHot dropFirst=True: 1 列分減る" $ do
let df = DX.fromNamedColumns
[("c", DX.fromList ["X", "Y", "Z" :: T.Text])]
df' = Reshape.oneHot True "c" df
length (DX.columnNames df') `shouldBe` 2 -- Y, Z (X drop)
-- ===========================================================================
-- Hanalyze.Stat.Effect (Phase 9)
-- ===========================================================================
describe "Hanalyze.Stat.Effect" $ do
it "cohenD: 同じ分布で 0、平均差 1 SD で ~1" $ do
let xs = LA.fromList [1, 2, 3, 4, 5] :: LA.Vector Double
ys = LA.fromList [1, 2, 3, 4, 5] :: LA.Vector Double
Eff.cohenD xs ys `shouldBe` 0.0
let zs = LA.fromList [2, 3, 4, 5, 6] :: LA.Vector Double -- shifted
d2 = Eff.cohenD zs xs
d2 `shouldSatisfy` (> 0.5)
it "hedgesG ≤ |cohenD| (small-sample correction)" $ do
let xs = LA.fromList [1, 2, 3, 4, 5] :: LA.Vector Double
ys = LA.fromList [3, 4, 5, 6, 7] :: LA.Vector Double
d = Eff.cohenD xs ys
g = Eff.hedgesG xs ys
abs g `shouldSatisfy` (<= abs d + 1e-9)
it "eta2: 完全分離で大、同分布で 0" $ do
let g1 = LA.fromList [1, 2, 3] :: LA.Vector Double
g2 = LA.fromList [10, 11, 12] :: LA.Vector Double
g3 = LA.fromList [20, 21, 22] :: LA.Vector Double
Eff.eta2 [g1, g2, g3] `shouldSatisfy` (> 0.9)
it "powerTTest: d = 0 で α 付近 (= ~0.05)" $ do
let p = Eff.powerTTest 30 0.05 0.0
p `shouldSatisfy` (\x -> abs (x - 0.05) < 0.01)
it "powerTTest: 大きい d で高い power" $ do
let p = Eff.powerTTest 30 0.05 0.8
p `shouldSatisfy` (> 0.85)
it "sampleSizeTTest: 小 d ほど大 n が必要" $ do
let n1 = Eff.sampleSizeTTest 0.80 0.05 0.2 -- small effect
n2 = Eff.sampleSizeTTest 0.80 0.05 0.5 -- medium
n3 = Eff.sampleSizeTTest 0.80 0.05 0.8 -- large
n1 `shouldSatisfy` (> n2)
n2 `shouldSatisfy` (> n3)
it "cramerV: 2×2 完全独立で ~0、強従属で大" $ do
Eff.cramerV 0 100 2 2 `shouldBe` 0.0
Eff.cramerV 50 100 2 2 `shouldSatisfy` (> 0.5)
-- ===========================================================================
-- Hanalyze.Model.DecisionTree (Phase 10)
-- ===========================================================================
describe "Hanalyze.Model.DecisionTree" $ do
it "fitDT: 線形分離可能なデータで perfect train accuracy" $ do
-- y = 0 if x[0] < 5, else 1
let xs = [[fromIntegral x] | x <- [1..10::Int]]
ys = [if x < 5 then 0 else 1 | x <- [1..10::Int]]
tree = DT.fitDT DT.defaultDTConfig xs ys
preds = map (DT.predictDT tree) xs
preds `shouldBe` ys
it "fitDT: 2D XOR-like パターン" $ do
let xs = [[0, 0], [0, 1], [1, 0], [1, 1]]
ys = [0, 1, 1, 0] -- XOR
tree = DT.fitDT DT.defaultDTConfig xs ys
preds = map (DT.predictDT tree) xs
preds `shouldBe` ys
it "predictDTProbs: leaf で確率 1.0、混合 leaf で fractional" $ do
let xs = [[1.0], [2.0], [3.0]]
ys = [0, 0, 1]
tree = DT.fitDT DT.defaultDTConfig xs ys
probs = DT.predictDTProbs tree [1.5]
-- x=1.5 should reach a leaf where most samples are class 0
probs `shouldSatisfy` (\m -> length m >= 1)
it "giniImpurity: 純粋クラスで 0、均等で 0.5 (2 クラス)" $ do
DT.giniImpurity [0, 0, 0, 0] `shouldBe` 0.0
DT.giniImpurity [1, 1, 1, 1] `shouldBe` 0.0
DT.giniImpurity [0, 0, 1, 1] `shouldBe` 0.5
it "maxDepth=1: shallow tree、underfit に近い" $ do
let cfg = DT.defaultDTConfig { DT.dtMaxDepth = Just 1 }
xs = [[fromIntegral x, fromIntegral y]
| x <- [1..5::Int], y <- [1..5::Int]]
ys = [if x + y > 5 then 1 else 0
| x <- [1..5::Int], y <- [1..5::Int]]
tree = DT.fitDT cfg xs ys
-- maxDepth=1 → 1 split, root = decision node
case tree of
DT.DNode {} -> True `shouldBe` True
_ -> expectationFailure "Expected DNode at root"
-- ===========================================================================
-- Hanalyze.Model.TimeSeries (Phase 11)
-- ===========================================================================
describe "Hanalyze.Model.TimeSeries" $ do
it "autocorrelation: lag 0 = 1.0" $ do
let y = LA.fromList [1.0, 2.0, 1.5, 2.5, 1.8]
acf = TS.autocorrelation 5 y
LA.atIndex acf 0 `shouldBe` 1.0
it "autocorrelation: 周期 4 の sin 系列で lag 4 が高い" $ do
let y = LA.fromList [sin (2*pi*fromIntegral t / 4) | t <- [0..40::Int]]
acf = TS.autocorrelation 8 y
LA.atIndex acf 4 `shouldSatisfy` (> 0.5)
it "fitAR(1) on quasi-AR(1) data: φ̂ in (0, 1)" $ do
-- Quasi-AR(1) with φ=0.7 + small deterministic perturbation
let phi = 0.7
go _ 0 = []
go prev n =
let yi = phi * prev + 0.01 * sin (fromIntegral n :: Double)
in yi : go yi (n - 1)
ys = take 100 (drop 50 (go 1.0 200)) -- burn-in 50
y = LA.fromList ys
fit = TS.fitAR 1 y
let phiHat = LA.atIndex (TS.arPhi fit) 0
-- AR(1) coefficient should be in (0, 1) for stationary positive AR
phiHat `shouldSatisfy` (\p -> p > 0 && p < 1)
it "forecastAR: h-step forecast 同 size" $ do
let y = LA.fromList [1.0, 1.5, 2.0, 2.5, 3.0, 2.5, 2.0, 1.5, 1.0, 1.5,
2.0, 2.5, 3.0, 2.5, 2.0]
fit = TS.fitAR 2 y
fc = TS.forecastAR fit y 5
LA.size fc `shouldBe` 5
it "differencing: y' length = n - 1" $ do
let y = LA.fromList [1.0, 3.0, 2.0, 5.0, 4.0]
d1 = TS.differencing y
LA.size d1 `shouldBe` 4
LA.toList d1 `shouldBe` [2.0, -1.0, 3.0, -1.0]
it "simpleExpSmoothing: α=1 で原系列、α=0 で初期値固定" $ do
let y = LA.fromList [1.0, 2.0, 3.0, 4.0, 5.0]
s1 = TS.simpleExpSmoothing 1.0 y -- α=1
s0 = TS.simpleExpSmoothing 0.0 y -- α=0
LA.toList s1 `shouldBe` LA.toList y -- exact
all (== 1.0) (LA.toList s0) `shouldBe` True
it "holtWinters: 線形 trend + 周期 4 を再現" $ do
let -- y_t = t + 5 sin(2πt/4)
ys = [ fromIntegral t + 3 * sin (2 * pi * fromIntegral t / 4)
| t <- [0 .. 39 :: Int] ]
y = LA.fromList ys
fit = TS.holtWinters TS.HWAdditive 4 y
fc = TS.hwForecast fit 4
LA.size fc `shouldBe` 4
-- Forecast at t=40,41,42,43 should be roughly 40 + sin pattern
LA.atIndex fc 0 `shouldSatisfy` (> 35)
it "stlDecompose: 周期成分が period 倍で繰り返す" $ do
let ys = [ fromIntegral t / 10 + 2 * sin (2 * pi * fromIntegral t / 4)
| t <- [0 .. 39 :: Int] ]
y = LA.fromList ys
(_trend, seasonal, _resid) = TS.stlDecompose 4 y
LA.size seasonal `shouldBe` LA.size y
-- ===========================================================================
-- Hanalyze.Model.Survival (Phase 12)
-- ===========================================================================
describe "Hanalyze.Model.Survival" $ do
it "kaplanMeier: 全 event observed で S(t) は単調減少" $ do
let samples = [ Surv.SurvSample t Surv.Observed | t <- [1, 2, 3, 4, 5] ]
km = Surv.kaplanMeier samples
length (Surv.kmrTimes km) `shouldBe` 5
let ss = Surv.kmrSurvival km
and (zipWith (>=) ss (tail ss)) `shouldBe` True
it "kaplanMeier: censored data でも非負の生存確率" $ do
let samples = [ Surv.SurvSample 1 Surv.Observed
, Surv.SurvSample 2 Surv.Censored
, Surv.SurvSample 3 Surv.Observed
, Surv.SurvSample 4 Surv.Observed
, Surv.SurvSample 5 Surv.Censored
]
km = Surv.kaplanMeier samples
all (>= 0) (Surv.kmrSurvival km) `shouldBe` True
all (<= 1) (Surv.kmrSurvival km) `shouldBe` True
it "nelsonAalen: 累積ハザードは monotone non-decreasing" $ do
let samples = [ Surv.SurvSample t Surv.Observed | t <- [1, 2, 3, 4, 5] ]
na = Surv.nelsonAalen samples
h = Surv.narCumHazard na
and (zipWith (<=) h (tail h)) `shouldBe` True
it "logRankTest: 同一分布で p > 0.05" $ do
let g1 = [ Surv.SurvSample t Surv.Observed | t <- [1, 2, 3, 4, 5] ]
g2 = [ Surv.SurvSample t Surv.Observed | t <- [1, 2, 3, 4, 5] ]
lr = Surv.logRankTest [g1, g2]
Surv.lrPValue lr `shouldSatisfy` (> 0.05)
it "logRankTest: 異なる分布で p < 0.05" $ do
let g1 = [ Surv.SurvSample t Surv.Observed | t <- [1, 2, 3, 4, 5] ]
g2 = [ Surv.SurvSample t Surv.Observed | t <- [10, 11, 12, 13, 14] ]
lr = Surv.logRankTest [g1, g2]
Surv.lrPValue lr `shouldSatisfy` (< 0.05)
it "coxPH: 共変量と event 時間に強い相関で β > 0" $ do
-- x が大きい個体が早く event を起こす想定の合成データ
let n = 30
xs = [ LA.fromList [fromIntegral i / fromIntegral n :: Double]
| i <- [1 .. n] ]
times = [ 30 - i | i <- [1 .. n] ] -- x 大きいほど early
samples = [ Surv.SurvSample (fromIntegral t) Surv.Observed
| t <- times ]
fit = Surv.coxPH xs samples
LA.atIndex (Surv.coxBeta fit) 0 `shouldSatisfy` (> 0)
-- ===========================================================================
-- Hanalyze.Stat.Interpret (Phase 13)
-- ===========================================================================
describe "Hanalyze.Stat.Interpret" $ do
it "permutationImportance: 重要 feature が高 importance" $ do
gen <- MWC.createSystemRandom
-- y = x_0 のみに依存、x_1 は無関係
let xs = [[fromIntegral i, fromIntegral (i * i)]
| i <- [1 .. 20 :: Int]]
ys = [head row | row <- xs]
predict zs = [head row | row <- zs]
score yt yp =
let mse = sum [(yt !! i - yp !! i) ^ (2 :: Int)
| i <- [0 .. length yt - 1]]
in negate mse -- higher (= 0) is better
cfg = (Interp.defaultPermutationConfig)
{ Interp.pcNRepeats = 5 }
r <- Interp.permutationImportance cfg predict score xs ys gen
let imps = Interp.piMeanImportance r
-- 0番目 feature の importance が高いはず (重要)
head imps `shouldSatisfy` (> 0.5 * (imps !! 1))
it "partialDependence: y = x[0] で PDP が grid に追従" $ do
let xs = [[fromIntegral i, 0.0] | i <- [1..10::Int]]
predict zs = [head row | row <- zs]
grid = [1.0, 5.0, 10.0]
pdp = Interp.partialDependence predict xs 0 grid
-- 各 grid 点での PD は値そのもの (= grid 値)
Interp.pdpMeanPredict pdp `shouldBe` grid
it "icePlot: 全 sample について curve を返す" $ do
let xs = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]
predict zs = [sum row | row <- zs]
grid = [0.0, 1.0, 2.0]
ice = Interp.icePlot predict xs 0 grid
length (Interp.iceCurves ice) `shouldBe` 3
length (Interp.iceFeatureValues ice) `shouldBe` 3
-- iceMean should equal partial dependence
length (Interp.iceMean ice) `shouldBe` 3
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Model.LM.Diagnostics (vs statsmodels OLS)" $ do
let xRaw = [1.0, 2, 3, 4, 5, 6, 7, 8, 9, 10]
yRaw = [3.7, 5.2, 6.8, 8.5, 9.9, 11.3, 13.1, 14.4, 15.8, 17.2]
x = LA.fromColumns
[ LA.konst 1.0 (length xRaw)
, LA.fromList xRaw
]
y = LA.fromList yRaw
fit = LM.fitLMVec x y
approx tol a b = abs (a - b) < tol
approxV tol expected actual =
length expected == LA.size actual &&
and (zipWith (approx tol) expected (LA.toList actual))
it "ciTValue: 95% / df=8 ≈ 2.306" $
LMD.ciTValue 0.95 8 `shouldSatisfy` approx 1e-3 2.306
it "lmStdErrors matches statsmodels (intercept, slope)" $
LMD.lmStdErrors x fit `shouldSatisfy`
approxV 1e-5 [0.09602188, 0.01547533]
it "lmCoefStats t / p match statsmodels" $ do
let cs = LMD.lmCoefStats x fit
length cs `shouldBe` 2
LMD.csTValue (head cs) `shouldSatisfy` approx 1e-3 23.8834
LMD.csTValue (cs !! 1) `shouldSatisfy` approx 1e-3 97.4768
LMD.csPValue (head cs) `shouldSatisfy` approx 1e-9 1.0062e-8
LMD.csPValue (cs !! 1) `shouldSatisfy` approx 1e-13 1.3699e-13
it "lmFStatistic matches statsmodels (F, p, df1, df2)" $ do
let fs = head (LMD.lmFStatistic x fit)
LMD.fsValue fs `shouldSatisfy` approx 1e-1 9501.7193
LMD.fsPValue fs `shouldSatisfy` approx 1e-13 1.3699e-13
LMD.fsDf1 fs `shouldBe` 1
LMD.fsDf2 fs `shouldBe` 8
it "lmInformationCriteria (R lm() convention, k = p + 1 with σ)" $ do
let ic = LMD.lmInformationCriteria fit
LMD.icLogLik ic `shouldSatisfy` approx 1e-5 6.547424
LMD.icAIC ic `shouldSatisfy` approx 1e-5 (-7.094848)
LMD.icBIC ic `shouldSatisfy` approx 1e-5 (-6.187093)
it "hatDiagonal matches statsmodels leverage" $
LMD.hatDiagonal x `shouldSatisfy`
approxV 1e-6
[ 0.34545455, 0.24848485, 0.17575758, 0.12727273, 0.10303030
, 0.10303030, 0.12727273, 0.17575758, 0.24848485, 0.34545455 ]
it "standardizedResiduals match statsmodels (resid_studentized_internal)" $
LMD.standardizedResiduals x fit `shouldSatisfy`
approxV 1e-5
[ -0.89534127, -0.90521501, -0.14722564, 1.31539084, 0.48257654
, -0.33234045, 1.88308584, 0.30394970, -0.57197652, -1.56684722 ]
it "cooksDistance matches statsmodels" $
LMD.cooksDistance x fit `shouldSatisfy`
approxV 1e-5
[ 0.21154283, 0.13546767, 0.00231098, 0.12616429, 0.01337487
, 0.00634342, 0.25856339, 0.00984992, 0.05408646, 0.64784992 ]
it "predictorStdDevs: intercept col=0, x col≈3.0277" $ do
let sds = LMD.predictorStdDevs x
LA.size sds `shouldBe` 2
(LA.toList sds !! 0) `shouldSatisfy` approx 1e-12 0
(LA.toList sds !! 1) `shouldSatisfy` approx 1e-6 3.027650
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Design.Orthogonal.listArraysWithSize" $ do
it "returns one entry per standard array" $
length OA.listArraysWithSize `shouldBe` length OA.standardArrays
it "L9 entry exposes runs / factors / levels" $ do
let l9meta = head [ m | m <- OA.listArraysWithSize
, OA.omName m == OA.oaName OA.l9 ]
OA.omRuns l9meta `shouldBe` 9
OA.omFactors l9meta `shouldBe` 4
OA.omLevels l9meta `shouldBe` [3, 3, 3, 3]
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Design.Taguchi extras" $ do
it "snRatioWithDetails: SmallerBetter on [1, 2, 3]" $ do
let d = TG.snRatioWithDetails TG.SmallerBetter [1.0, 2.0, 3.0]
TG.sdN d `shouldBe` 3
TG.sdMean d `shouldSatisfy` (\v -> abs (v - 2.0) < 1e-12)
TG.sdVariance d `shouldSatisfy` (\v -> abs (v - 1.0) < 1e-12)
-- η = -10 log10((1+4+9)/3) = -10 log10(14/3) ≈ -6.690
TG.sdSN d `shouldSatisfy` (\v -> abs (v - (-6.690)) < 1e-2)
it "factorEffectsTable: contributions sum to 1" $ do
let specs = [ OA.FactorSpec "A" [OA.LText "lo", OA.LText "hi"]
, OA.FactorSpec "B" [OA.LNumeric 0, OA.LNumeric 1]
]
case OA.assignFactors OA.l4 specs of
Right ad -> do
let sns = [10.0, 12.0, 11.0, 13.0]
ext = TG.factorEffectsTable ad sns
length ext `shouldBe` 2
let totalC = sum (map TG.feeContribution ext)
totalC `shouldSatisfy` (\v -> abs (v - 1.0) < 1e-12)
all ((>= 0) . TG.feeRange) ext `shouldBe` True
Left e -> expectationFailure (show e)
-- ─────────────────────────────────────────────────────────────────────
describe "Hanalyze.Design.Quality.processCapability" $ do
it "centred process with σ=1, USL=6, LSL=−6 → Cp ≈ 2.0, Cpk ≈ 2.0" $ do
-- 11-point symmetric sample around 0 with σ=1 (population)
let xs = LA.fromList [-1.5, -1.0, -0.5, 0.5, 1.0, 1.5,
1.5, 1.0, 0.5, -0.5, -1.0, -1.5]
cap = Quality.processCapability (-6) 6 xs
Quality.capCp cap `shouldSatisfy` (> 0)
-- For a centred sample, Cp == Cpk by symmetry.
abs (Quality.capCp cap - Quality.capCpk cap)
`shouldSatisfy` (< 1e-2)
it "shifted process: Cpk < Cp" $ do
let xs = LA.fromList [4.0, 4.5, 4.2, 4.8, 4.3, 4.6, 4.4, 4.7]
cap = Quality.processCapability 0 6 xs
Quality.capCp cap `shouldSatisfy` (> Quality.capCpk cap)
it "processCapabilityUpper: only USL → Cp == Cpk" $ do
let xs = LA.fromList [1.0, 1.2, 0.9, 1.1, 1.05, 0.95]
cap = Quality.processCapabilityUpper 2.0 xs
Quality.capCp cap `shouldBe` Quality.capCpk cap
describe "Hanalyze.Model.GLM diagnostics (request/090-AB)" $ do
let xsG = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] :: [Double]
ysG = [1.1, 2.9, 5.2, 6.8, 9.1, 11.0, 13.2, 14.8, 17.0, 19.1] :: [Double]
xMatG = LA.matrix 2 (concatMap (\x -> [1, x]) xsG)
yVecG = LA.fromList ysG
(frG, sigmaG) = GLM.fitGLMFull GLM.Gaussian GLM.Identity xMatG yVecG
betaG = head (LA.toColumns (Core.coefficients frG))
muVG = head (LA.toColumns (Core.fitted frG))
it "glmPearsonResiduals (Gaussian, V=1) == raw residuals" $ do
let pr = GLM.glmPearsonResiduals GLM.Gaussian yVecG muVG
rr = yVecG - muVG
LA.norm_2 (pr - rr) `shouldSatisfy` (< 1e-9)
it "glmDevianceResiduals (Gaussian) == sign(y-μ)·|y-μ|" $ do
let dr = GLM.glmDevianceResiduals GLM.Gaussian yVecG muVG
ref = LA.fromList
[ signum (y - m) * abs (y - m)
| (y, m) <- zip ysG (LA.toList muVG) ]
LA.norm_2 (dr - ref) `shouldSatisfy` (< 1e-9)
it "glmVariance: Binomial μ(1-μ); Poisson μ" $ do
GLM.glmVariance GLM.Binomial 0.3 `shouldBe` 0.3 * 0.7
GLM.glmVariance GLM.Poisson 4.0 `shouldBe` 4.0
it "predictGlmEtaWithSE: η = xᵀβ, SE > 0" $ do
let xNew = LA.fromList [1, 5.0]
(eta, se) = GLM.predictGlmEtaWithSE betaG sigmaG xNew
eta `shouldSatisfy` (\e -> abs (e - (1 + 2 * 5.0)) < 0.5)
se `shouldSatisfy` (> 0)
se `shouldSatisfy` (< 1)
it "predictGlmMuWithCI (Identity): half-width ≈ 1.96·SE" $ do
let xNew = LA.fromList [1, 4.5]
ci = GLM.predictGlmMuWithCI GLM.Identity 0.95 betaG sigmaG xNew
(_, se) = GLM.predictGlmEtaWithSE betaG sigmaG xNew
halfW = (GLM.gpHi ci - GLM.gpLo ci) / 2
abs (halfW - 1.96 * se) `shouldSatisfy` (< 1e-2)
it "predictGlmMuWithCI (Logit): CI stays in (0,1)" $ do
let xs2 = LA.matrix 2 (concatMap (\i -> [1, fromIntegral (i :: Int)]) [0..9])
ys2 = LA.fromList [0,1,0,1,0,1,0,1,0,1]
(_, sigma2) = GLM.fitGLMFull GLM.Binomial GLM.Logit xs2 ys2
beta2 = LA.fromList [0.0, 0.05]
xNew = LA.fromList [1, 5.0]
ci = GLM.predictGlmMuWithCI GLM.Logit 0.95 beta2 sigma2 xNew
GLM.gpMu ci `shouldSatisfy` (\v -> v > 0 && v < 1)
GLM.gpLo ci `shouldSatisfy` (\v -> v >= 0)
GLM.gpHi ci `shouldSatisfy` (\v -> v <= 1)
GLM.gpLo ci `shouldSatisfy` (< GLM.gpHi ci)
describe "Hanalyze.Model.GLMM SE (request/100)" $ do
-- Same fixture as the GLMM tests above (3 groups × 4 obs).
-- design X = [1, x] over the same 12 rows.
let xMat12 = LA.matrix 2
( concatMap (\v -> [1, v])
[1,2,3,4,1,2,3,4,1,2,3,4 :: Double] )
yVec12 = LA.fromList
[7.1,6.9,7.0,7.0, 5.0,4.9,5.1,5.0, 3.0,2.9,3.1,3.0]
gVec12 = V.fromList
["A","A","A","A","B","B","B","B","C","C","C","C" :: T.Text]
-- Inline group construction (mirrors Hanalyze.Model.GLMM.buildGroups
-- which is currently internal).
gLabels12 = V.fromList ["A", "B", "C"] :: V.Vector T.Text
gIdx12 = V.map
(\g -> case V.elemIndex g gLabels12 of
Just i -> i
Nothing -> 0)
gVec12
gSizes12 = V.fromList [4, 4, 4]
glmmRes = fitLME xMat12 yVec12 gIdx12 gLabels12 gSizes12
idx12 = gIdx12
it "glmmFixedSE: returns one SE per coefficient (length p)" $ do
let ses = glmmFixedSE xMat12 idx12 glmmRes
LA.size ses `shouldBe` LA.cols xMat12
it "glmmFixedSE: all SEs are positive" $ do
let ses = glmmFixedSE xMat12 idx12 glmmRes
mapM_ (\v -> v `shouldSatisfy` (> 0)) (LA.toList ses)
it "glmmFixedSE: σ_u → 0 reduces to OLS SE within tolerance" $ do
-- Force σ²_u to 0 (no random effects → OLS).
let resOLS = glmmRes { glmmRandVar = 0 }
sesG = glmmFixedSE xMat12 idx12 resOLS
-- Reference OLS SE from σ² (XᵀX)⁻¹.
xtx = LA.tr xMat12 LA.<> xMat12
covOLS = LA.scale (glmmResidVar resOLS) (LA.inv xtx)
sesOLS = LA.fromList
[ sqrt (LA.atIndex covOLS (i, i))
| i <- [0 .. LA.cols xMat12 - 1] ]
LA.norm_Inf (sesG - sesOLS) `shouldSatisfy` (< 1e-9)
it "glmmBLUPSE: one entry per group, all positive" $ do
let ses = glmmBLUPSE idx12 glmmRes
V.length ses `shouldBe` V.length (glmmGroups glmmRes)
mapM_ (\v -> v `shouldSatisfy` (> 0)) (V.toList ses)
it "glmmBLUPSE: shrinkage formula (1/σ²_u + n_j/σ²)⁻¹^½" $ do
let ses = glmmBLUPSE idx12 glmmRes
sig2u = glmmRandVar glmmRes
sig2 = glmmResidVar glmmRes
-- Group sizes: A=4, B=4, C=4 (balanced design)
expected = sqrt (1.0 / (1.0 / sig2u + 4.0 / sig2))
mapM_ (\(_, v) -> abs (v - expected) `shouldSatisfy` (< 1e-9))
(zip [0 :: Int ..] (V.toList ses))