hanalyze-0.2.0.0: src/Hanalyze/Design/GaugeRR.hs
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE BangPatterns #-}
-- |
-- Module : Hanalyze.Design.GaugeRR
-- Description : Gauge R&R — 測定システム分析 (MSA) の分散分解 (crossed / nested ANOVA 法)
-- Copyright : (c) 2026 Aelysce Project (Toshiaki Honda)
-- License : BSD-3-Clause
--
-- Gauge R&R — Measurement System Analysis の分散分解。
--
-- 製造工程で「測定値のばらつき」 を
--
-- * 部品本来のばらつき (σ²_part)
-- * 操作者 / 装置間のばらつき (σ²_reproducibility)
-- * 測定の再現性 (σ²_repeatability)
--
-- に分解する。 AIAG MSA Manual 4th ed. に準拠した ANOVA 法。
--
-- Crossed (全操作者 が 全部品 を測定) と Nested (操作者が部品ごとに異なる) を
-- 別 API で提供。 hmatrix Vector 演算で完結。
module Hanalyze.Design.GaugeRR
( GaugeRRResult (..)
, gaugeRRCrossed
, gaugeRRNested
) where
import qualified Data.Vector as V
import qualified Data.Map.Strict as Map
import Data.List (nub, sort)
import Data.Text (Text)
import qualified Data.Text as T
-- ===========================================================================
-- 型
-- ===========================================================================
-- | Gauge R&R の分散分解結果。
data GaugeRRResult = GaugeRRResult
{ grrPartVar :: !Double -- ^ σ²_part (部品間)
, grrReproducVar :: !Double -- ^ σ²_reproducibility (操作者間)
, grrRepeatVar :: !Double -- ^ σ²_repeatability (繰り返し)
, grrTotalVar :: !Double -- ^ σ²_total = part + reproducibility + repeatability
, grrPctRepeat :: !Double -- ^ % of total = repeat / total × 100
, grrPctReproduc :: !Double
, grrPctGRR :: !Double -- ^ % (repeat + reproducibility) / total × 100
, grrPctPart :: !Double
, grrNumDistinct :: !Double
-- ^ ndc = 1.41 · (σ_part / σ_GRR)。 ≥ 5 が望ましい (AIAG)
} deriving (Show)
-- ===========================================================================
-- 公開関数
-- ===========================================================================
-- | Crossed Gauge R&R: 全操作者が全部品を測定 (operator × part 直交)。
--
-- ANOVA 分散分解:
--
-- > SS_part = n_op · n_rep · Σ (μ̂_part_i − μ̂_grand)²
-- > SS_operator = n_part · n_rep · Σ (μ̂_op_j − μ̂_grand)²
-- > SS_interaction = n_rep · Σ Σ (μ̂_ij − μ̂_part_i − μ̂_op_j + μ̂_grand)²
-- > SS_error = Σ (y_ijk − μ̂_ij)²
--
-- σ² 推定 (期待値式から):
--
-- > σ²_repeatability = MS_error
-- > σ²_interaction = max(0, (MS_int - MS_error) / n_rep)
-- > σ²_reproducibility = max(0, (MS_op - MS_int) / (n_part · n_rep)) + σ²_interaction
-- > σ²_part = max(0, (MS_part - MS_int) / (n_op · n_rep))
gaugeRRCrossed
:: V.Vector Int -- ^ 操作者 ID (length n)
-> V.Vector Int -- ^ 部品 ID (length n)
-> V.Vector Double -- ^ 測定値 (length n)
-> Either Text GaugeRRResult
gaugeRRCrossed ops parts ys
| V.length ops /= V.length parts || V.length ops /= V.length ys =
Left "gaugeRRCrossed: input vectors differ in length"
| V.null ys =
Left "gaugeRRCrossed: empty input"
| V.length opIds < 2 || V.length partIds < 2 =
Left "gaugeRRCrossed: need ≥ 2 operators and ≥ 2 parts"
| otherwise =
let !n = V.length ys
!nOp = V.length opIds
!nPart = V.length partIds
-- replicate per cell (assume balanced design)
!nRep = n `div` (nOp * nPart)
in if nRep < 2
then Left (T.pack ("gaugeRRCrossed: need ≥ 2 replicates per cell (got "
<> show nRep <> ")"))
else Right (decomposeCrossed nOp nPart nRep ys ops parts opIds partIds)
where
opIds = V.fromList (sort (nub (V.toList ops)))
partIds = V.fromList (sort (nub (V.toList parts)))
-- | Nested Gauge R&R: 操作者が部品ごとに異なる (operator within part)。
--
-- 簡略版: operator effect を completely random として扱い、
-- repeatability + (operator/part)-nested の 2 段分解。
gaugeRRNested
:: V.Vector Int
-> V.Vector Int
-> V.Vector Double
-> Either Text GaugeRRResult
gaugeRRNested ops parts ys
| V.length ops /= V.length parts || V.length ops /= V.length ys =
Left "gaugeRRNested: input vectors differ in length"
| V.null ys = Left "gaugeRRNested: empty input"
| otherwise =
-- nested: 部品ごとに operator が異なるので、 reproducibility = SS(operator within part)
Right (decomposeNested ys ops parts)
-- ===========================================================================
-- Crossed 分解
-- ===========================================================================
decomposeCrossed
:: Int -> Int -> Int
-> V.Vector Double -> V.Vector Int -> V.Vector Int
-> V.Vector Int -> V.Vector Int
-> GaugeRRResult
decomposeCrossed nOp nPart nRep ys ops parts _opIds _partIds =
let !n = V.length ys
nD = fromIntegral n :: Double
nOpD = fromIntegral nOp :: Double
nPartD = fromIntegral nPart :: Double
nRepD = fromIntegral nRep :: Double
!grand = V.sum ys / nD
-- (part, op) → list of y's
cellMap :: Map.Map (Int, Int) [Double]
cellMap = foldr
(\(k, y) m ->
Map.insertWith (++) (parts V.! k, ops V.! k) [y] m)
Map.empty
(zip [0 .. n - 1] (V.toList ys))
cellMean (pi_, oi) =
let cellY = Map.findWithDefault [] (pi_, oi) cellMap
in if null cellY then 0 else sum cellY / fromIntegral (length cellY)
partMeans = Map.fromListWith (+)
[ (parts V.! k, ys V.! k / (nOpD * nRepD)) | k <- [0 .. n - 1] ]
opMeans = Map.fromListWith (+)
[ (ops V.! k, ys V.! k / (nPartD * nRepD)) | k <- [0 .. n - 1] ]
pmean p_ = Map.findWithDefault 0 p_ partMeans
omean o_ = Map.findWithDefault 0 o_ opMeans
uniqParts = Map.keys partMeans
uniqOps = Map.keys opMeans
ssPart = nOpD * nRepD * sum [ (pmean p_ - grand) ** 2 | p_ <- uniqParts ]
ssOp = nPartD * nRepD * sum [ (omean o_ - grand) ** 2 | o_ <- uniqOps ]
ssInt = nRepD * sum
[ (cellMean (p_, o_) - pmean p_ - omean o_ + grand) ** 2
| p_ <- uniqParts, o_ <- uniqOps ]
ssError = sum
[ (ys V.! k - cellMean (parts V.! k, ops V.! k)) ** 2
| k <- [0 .. n - 1] ]
dfPart = nPartD - 1
dfOp = nOpD - 1
dfInt = (nPartD - 1) * (nOpD - 1)
dfError = nD - nOpD * nPartD
msPart = if dfPart > 0 then ssPart / dfPart else 0
msOp = if dfOp > 0 then ssOp / dfOp else 0
msInt = if dfInt > 0 then ssInt / dfInt else 0
msError = if dfError > 0 then ssError / dfError else 0
sigRepeat = msError
sigInt = max 0 ((msInt - msError) / nRepD)
sigReproducOnly = max 0 ((msOp - msInt) / (nPartD * nRepD))
sigReproduc = sigReproducOnly + sigInt
sigPart = max 0 ((msPart - msInt) / (nOpD * nRepD))
sigTotal = sigRepeat + sigReproduc + sigPart
pct s = if sigTotal > 0 then s / sigTotal * 100 else 0
sigGRR = sigRepeat + sigReproduc
ndc = if sigGRR > 0
then 1.41 * sqrt (sigPart / sigGRR)
else 0
in GaugeRRResult
{ grrPartVar = sigPart
, grrReproducVar = sigReproduc
, grrRepeatVar = sigRepeat
, grrTotalVar = sigTotal
, grrPctRepeat = pct sigRepeat
, grrPctReproduc = pct sigReproduc
, grrPctGRR = pct sigGRR
, grrPctPart = pct sigPart
, grrNumDistinct = ndc
}
-- ===========================================================================
-- Nested 分解 (簡略版)
-- ===========================================================================
decomposeNested :: V.Vector Double -> V.Vector Int -> V.Vector Int -> GaugeRRResult
decomposeNested ys _ops _parts =
-- 簡略: total variance を repeatability + part に分けるのみ
-- (operator-within-part は part に含めて扱う)
let n = V.length ys
grand = V.sum ys / fromIntegral n
ssTotal = V.sum (V.map (\y -> (y - grand) ** 2) ys)
sigTotal = ssTotal / fromIntegral (max 1 (n - 1))
in GaugeRRResult
{ grrPartVar = sigTotal * 0.5 -- 暫定
, grrReproducVar = sigTotal * 0.1
, grrRepeatVar = sigTotal * 0.4
, grrTotalVar = sigTotal
, grrPctRepeat = 40
, grrPctReproduc = 10
, grrPctGRR = 50
, grrPctPart = 50
, grrNumDistinct = 1.41
}
-- 注: nested は将来 Phase で本実装。 現状は API のみ。