diff --git a/Statistics/Sample/Internal.hs b/Statistics/Sample/Internal.hs index 53c9a97a..a9ad7eac 100644 --- a/Statistics/Sample/Internal.hs +++ b/Statistics/Sample/Internal.hs @@ -14,9 +14,10 @@ module Statistics.Sample.Internal ( robustSumVar , sum + , sumF ) where -import Numeric.Sum (kbn, sumVector) +import qualified Numeric.Sum as Sum import Prelude hiding (sum) import Statistics.Function (square) import qualified Data.Vector.Generic as G @@ -26,5 +27,9 @@ robustSumVar m = sum . G.map (square . subtract m) {-# INLINE robustSumVar #-} sum :: (G.Vector v Double) => v Double -> Double -sum = sumVector kbn +sum = Sum.sumVector Sum.kbn {-# INLINE sum #-} + +sumF :: Foldable f => f Double -> Double +sumF = Sum.sum Sum.kbn +{-# INLINE sumF #-} diff --git a/Statistics/Test/Bartlett.hs b/Statistics/Test/Bartlett.hs index 1443bb7c..4dc9ad06 100644 --- a/Statistics/Test/Bartlett.hs +++ b/Statistics/Test/Bartlett.hs @@ -1,35 +1,56 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-| Module : Statistics.Test.Bartlett Description : Bartlett's test for homogeneity of variances. -Copyright : (c) Praneya Kumar, 2025 +Copyright : (c) Praneya Kumar, Alexey Khudyakov, 2025 License : BSD-3-Clause -Implements Bartlett's test to check if multiple groups have equal variances. -Assesses equality of variances assuming normal distribution, sensitive to non-normality. +Bartlett's test is used to check that multiple groups of observations +come from distributions with equal variances. This test assumes that +samples come from normal distribution. If this is not the case it may +simple test for non-normality and Levene's ("Statistics.Test.Levene") +is preferred + +>>> import qualified Data.Vector.Unboxed as VU +>>> import Statistics.Test.Bartlett +>>> :{ +let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99] + b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05] + c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98] +in bartlettTest [a,b,c] +:} +Right (Test {testSignificance = mkPValue 1.1254782518843598e-5, testStatistics = 22.789434813726768, testDistribution = chiSquared 2}) + -} module Statistics.Test.Bartlett ( bartlettTest, module Statistics.Distribution.ChiSquared ) where -import qualified Data.Vector.Generic as G -import qualified Data.Vector.Unboxed as U -import Statistics.Distribution (cumulative) +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as VU +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Primitive as VP +#if MIN_VERSION_vector(0,13,2) +import qualified Data.Vector.Strict as VV +#endif + +import Statistics.Distribution (complCumulative) import Statistics.Distribution.ChiSquared (chiSquared, ChiSquared(..)) import Statistics.Sample (varianceUnbiased) import Statistics.Types (mkPValue) import Statistics.Test.Types (Test(..)) --- | Perform Bartlett's test for equal variances. --- The input is a list of vectors, where each vector represents a group of observations. --- Returns Either an error message or a Test ChiSquared containing the test statistic and p-value. -bartlettTest :: [U.Vector Double] -> Either String (Test ChiSquared) +-- | Perform Bartlett's test for equal variances. The input is a list +-- of vectors, where each vector represents a group of observations. +bartlettTest :: VG.Vector v Double => [v Double] -> Either String (Test ChiSquared) bartlettTest groups - | length groups < 2 = Left "At least two groups are required for Bartlett's test." - | any ((< 2) . G.length) groups = Left "Each group must have at least two observations." - | any (<= 0) groupVariances = Left "All groups must have positive variance." - | otherwise = Right $ Test + | length groups < 2 = Left "At least two groups are required for Bartlett's test." + | any ((< 2) . VG.length) groups = Left "Each group must have at least two observations." + | any ((<= 0) . var) groupVariances = Left "All groups must have positive variance." + | otherwise = Right Test { testSignificance = pValue , testStatistics = tStatistic , testDistribution = chiDist @@ -37,56 +58,42 @@ bartlettTest groups where -- Number of groups k = length groups - -- Sample sizes for each group - ni = map G.length groups - ni' = map fromIntegral ni - + ni = map (fromIntegral . VG.length) groups -- Total number of observations across all groups - nTotal = sum ni - - -- Variance for each group (unbiased estimate) - groupVariances = map varianceUnbiased groups - - -- Pooled variance calculation - sumWeightedVars = sum [ (n - 1) * v | (n, v) <- zip ni' groupVariances ] - pooledVariance = sumWeightedVars / fromIntegral (nTotal - k) - + n_tot = sum $ fromIntegral . VG.length <$> groups + -- Variance estimates + groupVariances = toVar <$> groups + sumWeightedVars = sum [ (n - 1) * v | Var{sampleN=n, var=v} <- groupVariances ] + pooledVariance = sumWeightedVars / fromIntegral (n_tot - k) -- Numerator of Bartlett's statistic numerator = - fromIntegral (nTotal - k) * log pooledVariance - - sum [ (n - 1) * log v | (n, v) <- zip ni' groupVariances ] - + fromIntegral (n_tot - k) * log pooledVariance - + sum [ (n - 1) * log v | Var{sampleN=n, var=v} <- groupVariances ] -- Denominator correction term - sumReciprocals = sum [1 / (n - 1) | n <- ni'] + sumReciprocals = sum [1 / (n - 1) | n <- ni] denomCorrection = - 1 + (sumReciprocals - 1 / fromIntegral (nTotal - k)) / (3 * (fromIntegral k - 1)) + 1 + (sumReciprocals - 1 / fromIntegral (n_tot - k)) / (3 * (fromIntegral k - 1)) - -- Test statistic T + -- Test statistic and test distrubution tStatistic = max 0 $ numerator / denomCorrection - - -- Degrees of freedom and chi-squared distribution - df = k - 1 - chiDist = chiSquared df - pValue = mkPValue $ 1 - cumulative chiDist tStatistic - - --- Example usage: --- import qualified Data.Vector.Unboxed as U --- import Statistics.Test.Bartlett - --- main :: IO () --- main = do --- let a = U.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99] --- b = U.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05] --- c = U.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98] - --- case bartlettTest [a,b,c] of --- Left err -> putStrLn $ "Error: " ++ err --- Right test -> do --- putStrLn $ "Bartlett's Test Statistic: " ++ show (testStatistics test) --- putStrLn $ "P-Value: " ++ show (testSignificance test) - --- Sample Output --- Bartlett's Test Statistic: ~32 --- P-Value: ~1e-5 + chiDist = chiSquared (k - 1) + pValue = mkPValue $ complCumulative chiDist tStatistic +{-# SPECIALIZE bartlettTest :: [V.Vector Double] -> Either String (Test ChiSquared) #-} +{-# SPECIALIZE bartlettTest :: [VU.Vector Double] -> Either String (Test ChiSquared) #-} +{-# SPECIALIZE bartlettTest :: [VS.Vector Double] -> Either String (Test ChiSquared) #-} +{-# SPECIALIZE bartlettTest :: [VP.Vector Double] -> Either String (Test ChiSquared) #-} +#if MIN_VERSION_vector(0,13,2) +{-# SPECIALIZE bartlettTest :: [VV.Vector Double] -> Either String (Test ChiSquared) #-} +#endif + +-- Estimate of variance +data Var = Var + { sampleN :: !Double -- ^ N of elements + , var :: !Double -- ^ Sample variance + } + +toVar :: VG.Vector v Double => v Double -> Var +toVar xs = Var { sampleN = fromIntegral $ VG.length xs + , var = varianceUnbiased xs + } diff --git a/Statistics/Test/Levene.hs b/Statistics/Test/Levene.hs index ffa545c4..3fdc7c34 100644 --- a/Statistics/Test/Levene.hs +++ b/Statistics/Test/Levene.hs @@ -1,152 +1,153 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} - {-| Module : Statistics.Test.Levene Description : Levene's test for homogeneity of variances. -Copyright : (c) Praneya Kumar, 2025 +Copyright : (c) Praneya Kumar, Alexey Khudyakov, 2025 License : BSD-3-Clause -Implements Levene's test to check if multiple groups have equal variances. -Assesses equality of variances, robust to non-normality, and versatile with mean or median centering. +Levene's test used to check whether samples have equal variance. Null +hypothesis is all samples are from distributions with same variance +(homoscedacity). Test is robust to non-normality, and versatile with +mean or median centering. + +>>> import qualified Data.Vector.Unboxed as VU +>>> import Statistics.Test.Levene +>>> :{ +let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99] + b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05] + c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98] +in levenesTest Median [a, b, c] +:} +Right (Test {testSignificance = mkPValue 2.4315059672496814e-3, testStatistics = 7.584952754501659, testDistribution = fDistributionReal 2.0 27.0}) -} module Statistics.Test.Levene ( Center(..), levenesTest ) where -import qualified Data.Vector as V -import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Algorithms.Merge as VA -import Statistics.Distribution (cumulative) +import Control.Monad +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed as VU +import qualified Data.Vector.Generic as VG +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Primitive as VP +#if MIN_VERSION_vector(0,13,2) +import qualified Data.Vector.Strict as VV +#endif +import Statistics.Distribution (complCumulative) import Statistics.Distribution.FDistribution (fDistribution, FDistribution) -import Statistics.Types (mkPValue) +import Statistics.Types (mkPValue) import Statistics.Test.Types (Test(..)) -import qualified Statistics.Sample as Sample -import Control.Exception (assert) +import Statistics.Function (gsort) +import Statistics.Sample (mean) --- | Center calculation method -data Center = - Mean -- ^ Use arithmetic mean - | Median -- ^ Use median - | Trimmed Double -- ^ Trimmed mean with given proportion to cut from each end - deriving (Eq, Show) +import qualified Statistics.Sample.Internal as IS +import Statistics.Quantile --- | Trim data from both ends with error handling and performance optimization -trimboth :: (Ord a, Fractional a) - => V.Vector a - -> Double - -> Either String (V.Vector a) -trimboth vec prop - | prop < 0 || prop > 1 = Left "Proportion must be between 0 and 1" - | V.null vec = Right vec - | otherwise = do - let sorted = V.modify VA.sort vec - n = V.length sorted - lowerCut = floor (prop * fromIntegral n) - upperCut = n - lowerCut - assert (upperCut >= lowerCut) $ - Right $ V.slice lowerCut (upperCut - lowerCut) sorted --- | Calculate median using pre-sorted vector -vectorMedian :: (Fractional a, Ord a) - => V.Vector a - -> Either String a -vectorMedian vec - | V.null vec = Left "Empty vector in median calculation" - | otherwise = Right $ - if odd len - then sorted V.! mid - else (sorted V.! (mid - 1) + sorted V.! mid) / 2 - where - sorted = V.modify VA.sort vec - len = V.length sorted - mid = len `div` 2 +-- | Center calculation method +data Center + = Mean -- ^ Use arithmetic mean + | Median -- ^ Use median + | Trimmed !Double -- ^ Trimmed mean with given proportion to cut from each end + deriving (Eq, Show) -- | Main Levene's test function with full error handling -levenesTest :: Double -- ^ Significance level (alpha) - -> Center -- ^ Centering method - -> [V.Vector Double] -- ^ Input samples - -> Either String (Test FDistribution) -levenesTest alpha center samples - | alpha < 0 || alpha > 1 = Left "Significance level must be between 0 and 1" +levenesTest + :: (VG.Vector v Double) + => Center -- ^ Centering method + -> [v Double] -- ^ Input samples + -> Either String (Test FDistribution) +{-# INLINABLE levenesTest #-} +levenesTest center samples | length samples < 2 = Left "At least two samples required" + -- NOTE: We don't have nice way of computing mean of a list! | otherwise = do - processed <- mapM processSample samples - let (deviationsList, niList) = unzip processed - deviations = V.fromList deviationsList -- V.Vector (U.Vector Double) - ni = V.fromList niList -- V.Vector Int - zbari = V.map Sample.mean deviations -- V.Vector Double - k = V.length deviations - n = V.sum ni - zbar = V.sum (V.zipWith (\z n' -> z * fromIntegral n') zbari ni) / fromIntegral n - - -- Numerator: Sum over (ni * (zbari - zbar)^2) - numerator = V.sum $ V.zipWith (\n z -> fromIntegral n * (z - zbar) ** 2) ni zbari - - -- Denominator: Sum over sum((dev_ij - zbari)^2) - denominator = V.sum $ V.zipWith (\dev z -> U.sum (U.map (\x -> (x - z) ** 2) dev)) deviations zbari - + let residuals = computeResiduals center <$> samples + -- Average of all Z + let n_tot = sum $ VG.length . vecZ <$> residuals -- Total number of samples + let zbar = IS.sumF [ meanZ z * sampleN z + | z <- residuals] + / fromIntegral n_tot + -- Numerator: Sum over (ni * (Z[i] - Z)^2) + let numerator = IS.sumF [ sampleN z * sqr (meanZ z - zbar) + | z <- residuals] + -- Denominator: Sum over Σ((dev_ij - zbari)^2) + let denominator = IS.sumF + [ IS.sum $ VU.map (sqr . subtract (meanZ z)) (vecZ z) + | z <- residuals + ] -- Handle division by zero and invalid values - if denominator <= 0 || isNaN denominator || isInfinite denominator - then Left "Invalid denominator in W-statistic calculation" - else do - let wStat = (fromIntegral (n - k) / fromIntegral (k - 1)) * (numerator / denominator) - df1 = k - 1 - df2 = n - k - fDist = fDistribution df1 df2 - pVal = mkPValue $ 1 - cumulative fDist wStat - - -- Validate distribution parameters - if df1 < 1 || df2 < 1 - then Left "Invalid degrees of freedom" - else Right $ Test - { testStatistics = wStat - , testSignificance = pVal - , testDistribution = fDist - } + when (denominator <= 0 || isNaN denominator || isInfinite denominator) + $ Left "Invalid denominator in W-statistic calculation" + let wStat = (fromIntegral (n_tot - k) / fromIntegral (k - 1)) * (numerator / denominator) + fDist = fDistribution (k - 1) (n_tot - k) + Right Test { testStatistics = wStat + , testSignificance = mkPValue $ complCumulative fDist wStat + , testDistribution = fDist + } where - -- Process samples with error handling and optimized sorting - processSample vec = case center of - Mean -> do - let dev = V.map (abs . subtract (Sample.mean vec)) vec - return (U.convert dev, V.length vec) - - Median -> do - sortedVec <- Right $ V.modify VA.sort vec - m <- vectorMedian sortedVec - let dev = V.map (abs . subtract m) sortedVec - return (U.convert dev, V.length vec) - - Trimmed p -> do - trimmed_for_center_calculation <- trimboth vec p - let robust_center = Sample.mean trimmed_for_center_calculation - -- Calculate deviations for ALL ORIGINAL points from the robust_center - deviations_from_robust_center = V.map (abs . subtract robust_center) vec -- Use 'vec' (original data) - -- Return deviations and the ORIGINAL sample size - return (U.convert deviations_from_robust_center, V.length vec) -- Use 'V.length vec' - - --- Example usage: --- import qualified Data.Vector as V --- import LevenesTest (Center(..), levenesTest) --- import Statistics.Test.Types (testStatistics, testSignificance) --- import Statistics.Types (pValue) - --- main :: IO () --- main = do --- let a = V.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99] --- b = V.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05] --- c = V.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98] - --- case levenesTest (Trimmed 0.05) [a, b, c] of --- Left err -> putStrLn $ "Error: " ++ err --- Right test -> do --- putStrLn $ "Levene's W Statistic: " ++ show (testStatistics test) --- putStrLn $ "P-Value: " ++ show (pValue (testSignificance test)) --- putStrLn $ "Reject null hypothesis at α=0.05: " ++ show (testSignificance test < 0.05) + k = length samples -- Number of groups +{-# SPECIALIZE levenesTest :: Center -> [V.Vector Double] -> Either String (Test FDistribution) #-} +{-# SPECIALIZE levenesTest :: Center -> [VU.Vector Double] -> Either String (Test FDistribution) #-} +{-# SPECIALIZE levenesTest :: Center -> [VS.Vector Double] -> Either String (Test FDistribution) #-} +{-# SPECIALIZE levenesTest :: Center -> [VP.Vector Double] -> Either String (Test FDistribution) #-} +#if MIN_VERSION_vector(0,13,2) +{-# SPECIALIZE levenesTest :: Center -> [VV.Vector Double] -> Either String (Test FDistribution) #-} +#endif + +---------------------------------------------------------------- +-- Implementation +---------------------------------------------------------------- +-- | Trim data from both ends with error handling and performance optimization +trimboth :: (Ord a, Fractional a, VG.Vector v a) + => v a + -> Double + -> v a +{-# INLINE trimboth #-} +trimboth vec p + | p < 0 || p >= 0.5 = error "Statistics.Test.Levene: trimming: proportion must be between 0 and 0.5" + | VG.null vec = vec + | otherwise = VG.slice lowerCut (upperCut - lowerCut) sorted + where + n = VG.length vec + sorted = gsort vec + lowerCut = ceiling $ p * fromIntegral n + upperCut = n - lowerCut + +data Residuals = Residuals + { sampleN :: !Double + , meanZ :: !Double + , vecZ :: !(VU.Vector Double) + } + +computeResiduals + :: VG.Vector v Double + => Center + -> v Double + -> Residuals +{-# INLINE computeResiduals #-} +computeResiduals method xs = case method of + Mean -> + let c = mean xs + zs = VU.map (\x -> abs (x - c)) $ VU.convert xs + in makeR zs + Median -> + let c = median medianUnbiased xs + zs = VU.map (\x -> abs (x - c)) $ VU.convert xs + in makeR zs + Trimmed p -> + let trimmed = trimboth xs p + c = mean trimmed + zs = VU.map (\x -> abs (x - c)) $ VU.convert trimmed + in makeR zs + where + makeR zs = Residuals { sampleN = fromIntegral $ VU.length zs + , meanZ = mean zs + , vecZ = zs + } --- Sample Output --- Levene's W Statistic: 7.905 --- P-Value: 0.002 --- Reject null hypothesis at α=0.05: True +sqr :: Double -> Double +sqr x = x * x diff --git a/tests/Tests/Parametric.hs b/tests/Tests/Parametric.hs index 05f04ab9..6e5af8fe 100644 --- a/tests/Tests/Parametric.hs +++ b/tests/Tests/Parametric.hs @@ -15,7 +15,7 @@ import Statistics.Test.Bartlett tests :: Tst.TestTree -tests = testGroup "Parametric tests" (studentTTests ++ bartlettTests ++ leveneTests) +tests = testGroup "Parametric tests" [studentTTests, bartlettTests, leveneTests] -- 2 samples x 20 obs data -- @@ -84,8 +84,8 @@ testTTest name pVal test = Significant ] -studentTTests :: [Tst.TestTree] -studentTTests = concat +studentTTests :: Tst.TestTree +studentTTests = testGroup "StudentT test" $ concat [ -- R: t.test(sample1, sample2, alt="two.sided", var.equal=T) testTTest "two-sample t-test SamplesDiffer Student" (mkPValue 0.03410) (fromJust $ studentTTest SamplesDiffer sample1 sample2) @@ -112,45 +112,113 @@ studentTTests = concat -- Bartlett's Test ------------------------------------------------------------ -bartlettTests :: [TestTree] -bartlettTests = - let +bartlettTests :: TestTree +bartlettTests = testGroup "Bartlett's test" + [ testCase "a,b,c" $ testBartlettTest [a,b,c] 1.8027132567760222 0.40601846976301237 + , testCase "a,b" $ testBartlettTest [a,b] 0.005221063776321886 0.9423974408021293 + , testCase "a,c" $ testBartlettTest [a,c] 1.1531619271845452 0.2828882244527482 + , testCase "a,a" $ testBartlettTest [a,a] 0.0 1.0 + ] + where a = U.fromList [9.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99] b = U.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 9.36, 9.18, 8.67, 9.05] c = U.fromList [8.95, 8.12, 8.95, 8.85, 8.03, 8.84, 8.07, 8.98, 8.86, 8.98] - result = bartlettTest [a, b, c] - in case result of - Left err -> [testCase "Bartlett test - failed" (assertBool err False)] - Right test -> - [ testApproxEqual "Bartlett's Chi-Square Statistic" 1e-3 (testStatistics test) 1.802 - , testApproxEqual "Bartlett's P-Value" 1e-3 (pValue $ testSignificance test) 0.406 - , testEquality "Reject null hypothesis at alpha = 0.05" - (isSignificant (mkPValue 0.05) test) NotSignificant - ] + +testBartlettTest + :: [U.Vector Double] + -> Double + -> Double + -> IO () +testBartlettTest samples w p = do + r <- case bartlettTest samples of + Left _ -> error "Bartlett's test failed" + Right r -> pure r + approxEqual "W" 1e-9 (testStatistics r) w + approxEqual "p" 1e-9 (pValue $ testSignificance r) p ------------------------------------------------------------ -- Levene's Test (Trimmed Mean) ------------------------------------------------------------ --- Local helper for approximate equality -testApproxEqual :: String -> Double -> Double -> Double -> TestTree -testApproxEqual name epsilon actual expected = - testCase name $ - let diff = abs (actual - expected) - in assertBool (name ++ ": expected ≈ " ++ show expected ++ ", got " ++ show actual) (diff < epsilon) - -leveneTests :: [TestTree] -leveneTests = - let +leveneTests :: TestTree +leveneTests = testGroup "Levene test" + -- Statistics' value and p-values are computed using + [ testCase "a,b,c Mean" $ testLeveneTest [a,b,c] Mean 7.905194483442054 0.001983795817472731 + , testCase "a,b Mean" $ testLeveneTest [a,b] Mean 8.83873787256358 0.008149720958328811 + , testCase "a,a Mean" $ testLeveneTest [a,a] Mean 0.0 1.0 + , testCase "a,b,c Median" $ testLeveneTest [a,b,c] Median 7.584952754501659 0.002431505967249681 + , testCase "a,b Median" $ testLeveneTest [a,b] Median 8.461374333228711 0.009364737715584399 + , testCase "aL,bL Mean" $ testLeveneTest [aL,bL] Mean 5.84424549939465 0.01653410652558999 + , testCase "aL,bL Trimmed" $ testLeveneTest [aL,bL] (Trimmed 0.05) 8.368311226366314 0.004294953946529551 + ] + where a = V.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99] b = V.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05] c = V.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98] - result = levenesTest 0.05 (Trimmed 0.05) [a, b, c] - in case result of - Left err -> [testCase "Levene trimmed mean test - failed" (assertBool err False)] - Right test -> - [ testApproxEqual "Levene's W Statistic (Trimmed 0.05)" 1e-3 (testStatistics test) 7.905 - , testApproxEqual "Levene's P-Value (Trimmed 0.05)" 1e-3 (pValue $ testSignificance test) 0.002 - , testEquality "Reject null hypothesis at alpha = 0.05" - (isSignificant (mkPValue 0.05) test) Significant - ] + -- Large samples for testing trimmed + aL = V.fromList [ + -0.18919252, -1.62837673, 5.21332355, -0.00962043, -0.28417847, + -0.88128233, 1.49698436, 6.1780359 , -1.22301348, 3.34598245, + 5.33227264, -0.88732069, 0.14487346, 2.61060215, 4.22033907, + 2.53139215, -0.72131061, 0.53063607, -0.60510374, -0.73230842, + 1.54037043, -2.81103963, 3.40763063, 0.49005324, 2.13085513, + 5.68650547, 4.16397279, -0.17325097, 1.12664972, 4.23297516, + 4.15943436, -1.01452078, 2.40391646, 0.83019962, 0.29665879, + -3.83031046, -1.98576933, 1.5356527 , 1.30773365, 0.292818 , + 2.45877828, 1.06482289, -0.63241873, 1.58465379, 1.96577614, + 2.25791943, 4.13769848, -2.38595767, -0.65801423, -2.54007791, + 3.17428087, 4.32096964, 0.92240335, -2.38101319, 1.35692587, + 1.48279101, -0.04438309, 0.50296642, 2.08261495, 1.33181215, + -1.95427198, 4.95406809, 1.51294898, -2.68536129, -0.2441218 , + 2.41142613, 4.71051493, 2.66618697, 1.12668301, -0.25732583, + 1.25021838, -1.27523641, 5.01638744, 3.38864442, 0.17979744, + -0.88481645, 3.89346357, -0.51512217, -1.60542888, 0.88378679, + -2.12962732, -1.35989539, 5.09215112, -1.37442481, 0.83578405, + 0.13829571, 1.25171481, 3.60552158, -3.24051591, -0.44301834, + 0.78253445, 1.76098254, 1.79677434, -0.19010505, 3.07640466, + 3.02853882, 1.24849063, 4.84505382, 6.82274999, 2.24063474] + bL = V.fromList [ + 2.15584101, -2.74876744, -0.82231894, 1.97518087, 2.59280595, + 1.28703417, 2.40450278, 1.9761031 , 2.35186598, 1.15611047, + 2.26709318, 1.2832138 , -2.1486074 , 0.27563011, -0.51816861, + 0.89658424, 3.27069545, 1.72846646, 3.84454277, 5.58301459, + -0.40878188, 3.41602853, 1.1281526 , 0.9665913 , 0.76567084, + 1.69522855, 1.69133014, 0.70529264, 2.65243202, -1.0088019 , + -0.62431026, 3.76667396, 3.66225181, 0.73217579, 0.04478736, + 0.4169833 , 0.77065631, -1.31484093, 1.23858618, -0.08339456, + 3.14154286, 1.84358218, -0.53511423, -3.4919477 , 0.24076997, + 3.59381684, 1.99497806, 2.95499775, 1.67157731, 0.0214764 , + 3.32161612, -2.64762427, 0.06486472, 0.19653897, 1.34954235, + 1.18568747, -0.54434597, -3.35544223, 1.41933109, 0.95100195, + 2.7182116 , 1.1334068 , -0.95297806, -0.05421818, 1.42248799, + -3.96201277, -3.21309254, -0.21209211, 0.9689551 , 0.13526401, + -0.88656198, 0.41331783, -3.18766064, 4.34948246, 1.35656384, + 0.41920101, -0.46578994, 1.55181583, 2.43937014, 2.49040644, + 4.10505494, 1.68856296, 1.31503895, 0.41123368, 0.73242999, + 0.2804349 , -1.83494592, -0.31073195, 2.61185513, 2.91645094, + 1.26097638, 2.64197134, 3.88931972, 0.03783002, 2.55209729, + 3.46869549, 0.96348003, 2.27658242, 2.7613171 , -0.1372434 ] + + +testLeveneTest + :: [V.Vector Double] + -> Center + -> Double + -> Double + -> IO () +testLeveneTest samples center w p = do + r <- case levenesTest center samples of + Left _ -> error "Levene's test failed" + Right r -> pure r + approxEqual "W" 1e-9 (testStatistics r) w + approxEqual "p" 1e-9 (pValue $ testSignificance r) p + + +---------------------------------------------------------------- + +approxEqual :: String -> Double -> Double -> Double -> IO () +approxEqual name epsilon actual expected = + assertBool (name ++ ": expected ≈ " ++ show expected ++ ", got " ++ show actual) + (diff < epsilon) + where + diff = abs (actual - expected)