diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,4 @@ +# 0.1.1 + +* Add `*WithErrors` variants to regression functions. + diff --git a/regression-simple.cabal b/regression-simple.cabal --- a/regression-simple.cabal +++ b/regression-simple.cabal @@ -1,8 +1,8 @@ -cabal-version: 2.4 -name: regression-simple -version: 0.1 -synopsis: Simple linear and quadratic regression -category: Math +cabal-version: 2.4 +name: regression-simple +version: 0.1.1 +synopsis: Simple linear and quadratic regression +category: Math description: A simple package with a module for . @@ -14,11 +14,12 @@ . All specialized to @Double@. -license: BSD-3-Clause -license-file: LICENSE -maintainer: Oleg Grenrus -homepage: https://github.com/phadej/regression-simple -bug-reports: https://github.com/phadej/regression-simple/issues +license: BSD-3-Clause +license-file: LICENSE +maintainer: Oleg Grenrus +homepage: https://github.com/phadej/regression-simple +bug-reports: https://github.com/phadej/regression-simple/issues +extra-source-files: changelog.md tested-with: GHC ==7.4.2 || ==7.6.3 @@ -28,7 +29,10 @@ || ==8.2.2 || ==8.4.4 || ==8.6.5 - || ==8.8.1 + || ==8.8.4 + || ==8.10.7 + || ==9.0.2 + || ==9.2.1 source-repository head type: git @@ -39,5 +43,5 @@ hs-source-dirs: src exposed-modules: Math.Regression.Simple build-depends: - , base >=4.5 && <4.14 + , base >=4.5 && <4.17 , vector ^>=0.12.0.0 diff --git a/src/Math/Regression/Simple.hs b/src/Math/Regression/Simple.hs --- a/src/Math/Regression/Simple.hs +++ b/src/Math/Regression/Simple.hs @@ -1,11 +1,15 @@ +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE GADTs #-} module Math.Regression.Simple ( -- * Regressions linear, + linearWithErrors, quadratic, + quadraticWithErrors, quadraticAndLinear, + quadraticAndLinearWithErrors, -- * Operations Add (..), Eye (..), @@ -35,6 +39,37 @@ -- $setup -- >>> :set -XTypeApplications +-- +-- >>> import Numeric (showFFloat) +-- +-- Don't show too much decimal digits +-- +-- >>> showDouble x = showFFloat @Double (Just (min 5 (5 - ceiling (logBase 10 x)))) x +-- >>> showDouble 123.456 "" +-- "123.46" +-- +-- >>> showDouble 1234567 "" +-- "1234567" +-- +-- >>> showDouble 123.4567890123456789 "" +-- "123.46" +-- +-- >>> showDouble 0.0000000000000012345 "" +-- "0.00000" +-- +-- >>> newtype PP a = PP a +-- >>> class Show' a where showsPrec' :: Int -> a -> ShowS +-- >>> instance Show' a => Show (PP a) where showsPrec d (PP x) = showsPrec' d x +-- >>> instance Show' Double where showsPrec' d x = if x < 0 then showParen (d > 6) (showChar '-' . showDouble (negate x)) else showDouble x +-- >>> instance (Show' a, Show' b) => Show' (a, b) where showsPrec' d (x, y) = showParen True $ showsPrec' 0 x . showString ", " . showsPrec' 0 y +-- >>> instance Show' V2 where showsPrec' d (V2 x y) = showParen (d > 10) $ showString "V2 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y +-- >>> instance Show' V3 where showsPrec' d (V3 x y z) = showParen (d > 10) $ showString "V3 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y . showChar ' ' . showsPrec' 11 z +-- +-- Inputs: +-- >>> let input1 = [(0, 1), (1, 3), (2, 5)] +-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)] +-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)] +-- ------------------------------------------------------------------------------- -- Classes @@ -87,7 +122,7 @@ -- -2.0 -- zerosLin :: V2 -> Double -zerosLin (V2 a b) = - b / a +zerosLin (V2 a b) = negate (b / a) -- | Solve quadratic equation. -- @@ -104,6 +139,7 @@ -- -- >>> zerosQuad (V3 1 (-2) 1) -- Right (1.0,1.0) +-- zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double) zerosQuad (V3 a b c) | delta < 0 = Left ((-b/da) :+ (-sqrtNDelta/da), (-b/da) :+ (sqrtNDelta/da)) @@ -298,20 +334,65 @@ -- but overloaded to work with boxed and unboxed 'Vector's. -- -- >>> let input1 = [(0, 1), (1, 3), (2, 5)] --- >>> linear input1 --- V2 2.0 1.0 +-- >>> PP $ linear input1 +-- V2 2.0000 1.00000 -- -- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)] --- >>> linear input2 --- V2 2.0063237774030345 0.8868465430016883 +-- >>> PP $ linear input2 +-- V2 2.0063 0.88685 -- linear :: (Foldable' xs x, IsDoublePair x) => xs -> V2 -linear data_ = mult (inv2 (M22 x n x2 x)) (V2 y xy) +linear = fst . linearWithErrors + +-- | Like 'linear' but also return parameters' standard errors. +-- +-- To get confidence intervals you should multiply the errors +-- by @quantile (studentT (n - 2)) ci'@ from @statistics@ package +-- or similar. +-- For big @n@ using value 1 gives 68% interval and using value 2 gives 95% confidence interval. +-- See https://en.wikipedia.org/wiki/Student%27s_t-distribution#Table_of_selected_values +-- (@quantile@ calculates one-sided values, you need two-sided, thus adjust @ci@ value). +-- +-- The first input is perfect fit: +-- +-- >>> PP $ linearWithErrors input1 +-- (V2 2.0000 1.00000, V2 0.00000 0.00000) +-- +-- The second input is quite good: +-- +-- >>> PP $ linearWithErrors input2 +-- (V2 2.0063 0.88685, V2 0.09550 0.23826) +-- +-- But the third input isn't so much, +-- standard error of a slope argument is 20%. +-- +-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)] +-- >>> PP $ linearWithErrors input3 +-- (V2 3.0000 1.00000, V2 0.63246 1.1832) +-- +-- @since 0.1.1 +-- +linearWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> (V2, V2) +linearWithErrors = linearImpl . kahan2 + +linearImpl :: Kahan2 -> (V2, V2) +linearImpl (K2 n' (V2 x _) (V2 x2 _) (V2 y _) (V2 y2 _) (V2 xy _)) = + (params, errors) where - K2 n' (V2 x _) (V2 x2 _) (V2 y _) (V2 xy _) = kahan2 data_ n :: Double n = fromIntegral n' + matrix@(M22 a11 _ _ a22) = inv2 (M22 x2 x x n) + params@(V2 a b) = mult matrix (V2 xy y) + + errors = V2 sa sb + + -- ensure that error is always non-negative. + -- Due rounding errors, in perfect fit situations it can be slightly negative. + err = max 0 (y2 - a * xy - b * y) + sa = sqrt (a11 * err / (n - 2)) + sb = sqrt (a22 * err / (n - 2)) + -- | Quadratic regression. -- -- The type is @@ -327,36 +408,79 @@ -- V3 0.0 2.0 1.0 -- -- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)] --- >>> quadratic input2 --- V3 (-5.886346291028133e-3) 2.0312938469708826 0.8715454176158062 +-- >>> PP $ quadratic input2 +-- V3 (-0.00589) 2.0313 0.87155 -- -- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)] --- >>> quadratic input3 --- V3 1.0 0.0 1.999999999999993 +-- >>> PP $ quadratic input3 +-- V3 1.00000 0.00000 2.0000 -- quadratic :: (Foldable' xs x, IsDoublePair x) => xs -> V3 quadratic data_ = mult (inv3 (M33 x2 x n x3 x2 x x4 x3 x2)) (V3 y xy x2y) where - K3 n' (V2 x _) (V2 x2 _) (V2 x3 _) (V2 x4 _) (V2 y _) (V2 xy _) (V2 x2y _) = kahan3 data_ + K3 n' (V2 x _) (V2 x2 _) (V2 x3 _) (V2 x4 _) (V2 y _) (V2 xy _) (V2 x2y _) (V2 _y2 _) = kahan3 data_ + n :: Double n = fromIntegral n' --- | Do both linear and quadratic regression in one data scan. +-- | Like 'quadratic' but also return parameters' standard errors. -- --- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)] --- >>> quadraticAndLinear input2 --- (V3 (-5.886346291028133e-3) 2.0312938469708826 0.8715454176158062,V2 2.0063237774030345 0.8868465430016883) +-- >>> PP $ quadraticWithErrors input2 +-- (V3 (-0.00589) 2.0313 0.87155, V3 0.09281 0.41070 0.37841) -- -quadraticAndLinear :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V2) -quadraticAndLinear data_ = - ( mult (inv3 (M33 x2 x n x3 x2 x x4 x3 x2)) (V3 y xy x2y) - , mult (inv2 (M22 x n x2 x)) (V2 y xy) - ) +-- >>> PP $ quadraticWithErrors input3 +-- (V3 1.00000 0.00000 2.0000, V3 0.00000 0.00000 0.00000) +-- +-- @since 0.1.1 +-- +quadraticWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V3) +quadraticWithErrors = quadraticImpl . kahan3 + +quadraticImpl :: Kahan3 -> (V3, V3) +quadraticImpl (K3 n' (V2 x _) (V2 x2 _) (V2 x3 _) (V2 x4 _) (V2 y _) (V2 xy _) (V2 x2y _) (V2 y2 _)) = + (params, errors) where - K3 n' (V2 x _) (V2 x2 _) (V2 x3 _) (V2 x4 _) (V2 y _) (V2 xy _) (V2 x2y _) = kahan3 data_ n :: Double n = fromIntegral n' + matrix@(M33 a11 _ _ + _ a22 _ + _ _ a33) = inv3 (M33 x4 x3 x2 x3 x2 x x2 x n) + + params@(V3 a b c) = mult matrix (V3 x2y xy y) + + errors = V3 sa sb sc + + err = max 0 (y2 - a * x2y - b * xy - c * y) + sa = sqrt (a11 * err / (n - 3)) + sb = sqrt (a22 * err / (n - 3)) + sc = sqrt (a33 * err / (n - 3)) + +-- | Do both linear and quadratic regression in one data scan. +-- +-- >>> PP $ quadraticAndLinear input2 +-- (V3 (-0.00589) 2.0313 0.87155, V2 2.0063 0.88685) +-- +quadraticAndLinear :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V2) +quadraticAndLinear = fst . quadraticAndLinearWithErrors + +-- | Like 'quadraticAndLinear' but also return parameters' standard errors +-- +-- >>> PP $ quadraticAndLinearWithErrors input2 +-- ((V3 (-0.00589) 2.0313 0.87155, V2 2.0063 0.88685), (V3 0.09281 0.41070 0.37841, V2 0.09550 0.23826)) +-- +-- @since 0.1.1 +-- +quadraticAndLinearWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> ((V3, V2), (V3, V2)) +quadraticAndLinearWithErrors data_ = + ((paramsQ, paramsL), (errorsQ, errorsL)) + where + k3@(K3 n x x2 x3 x4 y xy x2y y2) = kahan3 data_ + k2 = K2 n x x2 y y2 xy + + (paramsL, errorsL) = linearImpl k2 + (paramsQ, errorsQ) = quadraticImpl k3 + ------------------------------------------------------------------------------- -- Input ------------------------------------------------------------------------------- @@ -391,11 +515,12 @@ , k2x :: {-# UNPACK #-} !V2 , k2x2 :: {-# UNPACK #-} !V2 , k2y :: {-# UNPACK #-} !V2 + , k2y2 :: {-# UNPACK #-} !V2 , k2xy :: {-# UNPACK #-} !V2 } zeroKahan2 :: Kahan2 -zeroKahan2 = K2 0 zero zero zero zero +zeroKahan2 = K2 0 zero zero zero zero zero -- | https://en.wikipedia.org/wiki/Kahan_summation_algorithm addKahan :: V2 -> Double -> V2 @@ -406,11 +531,12 @@ kahan2 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan2 kahan2 = foldl' f zeroKahan2 where - f (K2 n x x2 y xy) uv = withDP uv $ \u v -> K2 + f (K2 n x x2 y y2 xy) uv = withDP uv $ \u v -> K2 (succ n) (addKahan x u) (addKahan x2 (u * u)) (addKahan y v) + (addKahan y2 (v * v)) (addKahan xy (u * v)) ------------------------------------------------------------------------------- @@ -418,22 +544,23 @@ ------------------------------------------------------------------------------- data Kahan3 = K3 - { k3n :: {-# UNPACK #-} !Int - , k3x :: {-# UNPACK #-} !V2 - , k3x2 :: {-# UNPACK #-} !V2 - , k3x3 :: {-# UNPACK #-} !V2 - , k3x4 :: {-# UNPACK #-} !V2 - , k3y :: {-# UNPACK #-} !V2 - , k3xy :: {-# UNPACK #-} !V2 + { k3n :: {-# UNPACK #-} !Int + , k3x :: {-# UNPACK #-} !V2 + , k3x2 :: {-# UNPACK #-} !V2 + , k3x3 :: {-# UNPACK #-} !V2 + , k3x4 :: {-# UNPACK #-} !V2 + , k3y :: {-# UNPACK #-} !V2 + , k3xy :: {-# UNPACK #-} !V2 , k3x2y :: {-# UNPACK #-} !V2 + , ky2 :: {-# UNPACK #-} !V2 } zeroKahan3 :: Kahan3 -zeroKahan3 = K3 0 zero zero zero zero zero zero zero +zeroKahan3 = K3 0 zero zero zero zero zero zero zero zero kahan3 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan3 kahan3 = foldl' f zeroKahan3 where - f (K3 n x x2 x3 x4 y xy x2y) uv = withDP uv $ \u v -> + f (K3 n x x2 x3 x4 y xy x2y y2) uv = withDP uv $ \u v -> let u2 = u * u in K3 (succ n) @@ -444,3 +571,4 @@ (addKahan y v) (addKahan xy (u * v)) (addKahan x2y (u2 * v)) + (addKahan y2 (v * v))