diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,4 +1,11 @@ +# 0.2 + +Large refactoring: + +* Add `*WithWeights` variants +* Move linear algebra minilib into own module +* Expose accumulators + # 0.1.1 * Add `*WithErrors` variants to regression functions. - diff --git a/gnuplot/linear.dat b/gnuplot/linear.dat new file mode 100644 --- /dev/null +++ b/gnuplot/linear.dat @@ -0,0 +1,20 @@ + 1.33343 8.43365 0.20000 1.20000 + 2.46969 11.46688 0.30000 1.40000 + 2.80364 15.59652 0.20000 1.60000 + 4.00046 15.56544 0.30000 1.20000 + 5.18323 20.55600 0.20000 1.40000 + 5.66586 21.16630 0.30000 1.60000 + 6.82541 29.06629 0.20000 1.20000 + 7.64329 33.09135 0.30000 1.40000 + 8.80296 34.41986 0.20000 1.60000 + 9.65792 35.80643 0.30000 1.20000 + 10.63990 39.34492 0.20000 1.40000 + 12.48292 39.14740 0.30000 1.60000 + 12.83448 44.32188 0.20000 1.20000 + 13.93802 45.79094 0.30000 1.40000 + 14.98339 49.75683 0.20000 1.60000 + 16.09785 53.33996 0.30000 1.20000 + 16.74771 55.50356 0.20000 1.40000 + 17.81197 59.54537 0.30000 1.60000 + 18.56691 61.87966 0.20000 1.20000 + 20.64801 65.51835 0.30000 1.40000 diff --git a/gnuplot/quad.dat b/gnuplot/quad.dat new file mode 100644 --- /dev/null +++ b/gnuplot/quad.dat @@ -0,0 +1,20 @@ + 1.23669 1.90145 0.20000 1.20000 + 2.11352 -0.19517 0.30000 1.40000 + 2.98278 -1.34899 0.20000 1.60000 + 3.87135 -6.42140 0.30000 1.20000 + 5.50783 -6.20445 0.20000 1.40000 + 6.74626 -9.14974 0.30000 1.60000 + 7.41109 -12.14659 0.20000 1.20000 + 8.44699 -13.34493 0.30000 1.40000 + 9.34861 -14.56375 0.20000 1.60000 + 10.16446 -15.58626 0.30000 1.20000 + 11.17560 -16.32637 0.20000 1.40000 + 11.85809 -17.63580 0.30000 1.60000 + 12.88762 -17.86456 0.20000 1.20000 + 14.05061 -17.48367 0.30000 1.40000 + 15.15438 -18.61585 0.20000 1.60000 + 16.42014 -17.24635 0.30000 1.20000 + 17.00652 -17.96829 0.20000 1.40000 + 18.07526 -12.69336 0.30000 1.60000 + 18.63128 -18.25398 0.20000 1.20000 + 20.70280 -12.68441 0.30000 1.40000 diff --git a/regression-simple.cabal b/regression-simple.cabal --- a/regression-simple.cabal +++ b/regression-simple.cabal @@ -1,6 +1,6 @@ cabal-version: 2.4 name: regression-simple -version: 0.1.1 +version: 0.2 synopsis: Simple linear and quadratic regression category: Math description: @@ -19,9 +19,15 @@ maintainer: Oleg Grenrus homepage: https://github.com/phadej/regression-simple bug-reports: https://github.com/phadej/regression-simple/issues -extra-source-files: changelog.md +extra-source-files: + changelog.md + gnuplot/linear.dat + gnuplot/quad.dat + tested-with: - GHC ==7.4.2 + GHC ==7.0.4 + || ==7.2.2 + || ==7.4.2 || ==7.6.3 || ==7.8.4 || ==7.10.3 @@ -32,16 +38,60 @@ || ==8.8.4 || ==8.10.7 || ==9.0.2 - || ==9.2.1 + || ==9.2.5 + || ==9.4.4 source-repository head type: git location: https://github.com/phadej/regression-simple library + default-language: Haskell2010 + ghc-options: -Wall + hs-source-dirs: src + exposed-modules: + Math.Regression.Simple + Math.Regression.Simple.LinAlg + Numeric.KBN + + build-depends: + , base >=4.3 && <4.18 + , deepseq + + if !impl(ghc >=8.0) + build-depends: semigroups >=0.18.5 && <0.21 + + x-docspec-extra-packages: math-functions statistics ad + +test-suite generate-test-data default-language: Haskell2010 - hs-source-dirs: src - exposed-modules: Math.Regression.Simple + ghc-options: -Wall + hs-source-dirs: test + type: exitcode-stdio-1.0 + main-is: generate-test-data.hs build-depends: - , base >=4.5 && <4.17 - , vector ^>=0.12.0.0 + , base + , splitmix ^>=0.1.0.4 + +test-suite regression-simple-tests + default-language: Haskell2010 + ghc-options: -Wall + hs-source-dirs: test + type: exitcode-stdio-1.0 + main-is: regression-simple-tests.hs + build-depends: + , base + , regression-simple + + build-depends: + , tasty ^>=1.4.0.1 + , tasty-hunit ^>=0.10.0.3 + + if impl(ghc >=7.4) + build-depends: + , ad >=4.4.1 && <4.6 + , math-functions ^>=0.3.4.2 + , statistics ^>=0.10.2.0 || ^>=0.15.2.0 || ^>=0.16.0.1 + + if !impl(ghc >=8.0) + build-depends: semigroups 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 @@ -2,573 +2,1206 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE GADTs #-} -module Math.Regression.Simple ( - -- * Regressions - linear, - linearWithErrors, - quadratic, - quadraticWithErrors, - quadraticAndLinear, - quadraticAndLinearWithErrors, - -- * Operations - Add (..), - Eye (..), - Mult (..), - Det (..), - Inv (..), - -- * Zeros - zerosLin, - zerosQuad, - optimaQuad, - -- * Two dimensions - V2 (..), - M22 (..), - -- * Three dimensions - V3 (..), - M33 (..), - -- * Auxiliary classes - Foldable' (..), - IsDoublePair (..), - ) where - -import Data.Complex (Complex (..)) - -import qualified Data.List as L -import qualified Data.Vector as V -import qualified Data.Vector.Unboxed as U - --- $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 -------------------------------------------------------------------------------- - --- | Addition -class Add a where - zero :: a - add :: a -> a -> a - --- | Identity -class Eye a where - eye :: a - --- | Multiplication of different things. -class Eye a => Mult a b c | a b -> c where - mult :: a -> b -> c - --- | Determinant -class Eye a => Det a where - det :: a -> Double - --- | Inverse -class Det a => Inv a where - inv :: a -> a - -infixl 6 `add` -infixl 7 `mult` - -instance Eye Double where - eye = 1 - -instance Add Double where - zero = 0 - add = (+) - -instance Det Double where - det = id - -instance Inv Double where - inv = recip - -------------------------------------------------------------------------------- --- Zeros -------------------------------------------------------------------------------- - --- | Solve linear equation. --- --- >>> zerosLin (V2 1 2) --- -2.0 --- -zerosLin :: V2 -> Double -zerosLin (V2 a b) = negate (b / a) - --- | Solve quadratic equation. --- --- >>> zerosQuad (V3 2 0 (-1)) --- Right (-0.7071067811865476,0.7071067811865476) --- --- >>> zerosQuad (V3 2 0 1) --- Left ((-0.0) :+ (-0.7071067811865476),(-0.0) :+ 0.7071067811865476) --- --- Double root is not treated separately: --- --- >>> zerosQuad (V3 1 0 0) --- Right (-0.0,0.0) --- --- >>> 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)) - | otherwise = Right ((- b - sqrtDelta) / da, (-b + sqrtDelta) / da) - where - delta = b*b - 4 * a * c - sqrtDelta = sqrt delta - sqrtNDelta = sqrt (- delta) - da = 2 * a - --- | Find an optima point. --- --- >>> optimaQuad (V3 1 (-2) 0) --- 1.0 --- --- compare to --- --- >>> zerosQuad (V3 1 (-2) 0) --- Right (0.0,2.0) --- -optimaQuad :: V3 -> Double -optimaQuad (V3 a b _) = zerosLin (V2 (2 * a) b) - -------------------------------------------------------------------------------- --- 2 dimensions -------------------------------------------------------------------------------- - --- | 2d vector. Strict pair of 'Double's. --- --- Also used to represent linear polynomial: @V2 a b@ \(= a x + b\). --- -data V2 = V2 !Double !Double - deriving (Eq, Show) - -instance Add V2 where - zero = V2 0 0 - add (V2 x y) (V2 x' y') = V2 (x + x') (y + y') - {-# INLINE zero #-} - {-# INLINE add #-} - -instance Mult Double V2 V2 where - mult k (V2 x y) = V2 (k * x) (k * y) - {-# INLINE mult #-} - --- | 2×2 matrix. -data M22 = M22 !Double !Double !Double !Double - deriving (Eq, Show) - -instance Add M22 where - zero = M22 0 0 0 0 - add (M22 a b c d) (M22 a' b' c' d') = M22 (a + a') (b + b') (c + c') (d + d') - {-# INLINE zero #-} - {-# INLINE add #-} - -instance Eye M22 where - eye = M22 1 0 0 1 - {-# INLINE eye #-} - -instance Det M22 where det = det2 -instance Inv M22 where inv = inv2 - -instance Mult Double M22 M22 where - mult k (M22 a b c d) = M22 (k * a) (k * b) (k * c) (k * d) - {-# INLINE mult #-} - -instance Mult M22 V2 V2 where - mult (M22 a b c d) (V2 u v) = V2 (a * u + b * v) (c * u + d * v) - {-# INLINE mult #-} - --- | >>> M22 1 2 3 4 `mult` eye @M22 --- M22 1.0 2.0 3.0 4.0 -instance Mult M22 M22 M22 where - mult (M22 a b c d) (M22 x y z w) = M22 - (a * x + b * z) (a * y + b * w) - (c * x + d * z) (c * y + d * w) - {-# INLINE mult #-} - -det2 :: M22 -> Double -det2 (M22 a b c d) = a * d - b * c -{-# INLINE det2 #-} - -inv2 :: M22 -> M22 -inv2 m@(M22 a b c d) = M22 - ( d / det) (- b / det) - (- c / det) ( a / det) - where - det = det2 m -{-# INLINE inv2 #-} - -------------------------------------------------------------------------------- --- 3 dimensions -------------------------------------------------------------------------------- - --- | 3d vector. Strict triple of 'Double's. --- --- Also used to represent quadratic polynomial: @V3 a b c@ \(= a x^2 + b x + c\). -data V3 = V3 !Double !Double !Double - deriving (Eq, Show) - -instance Add V3 where - zero = V3 0 0 0 - add (V3 x y z) (V3 x' y' z') = V3 (x + x') (y + y') (z + z') - {-# INLINE zero #-} - {-# INLINE add #-} - -instance Mult Double V3 V3 where - mult k (V3 x y z) = V3 (k * x) (k * y) (k * z) - {-# INLINE mult #-} - --- | 3×3 matrix. -data M33 = M33 - !Double !Double !Double - !Double !Double !Double - !Double !Double !Double - deriving (Eq, Show) - -instance Add M33 where - zero = M33 0 0 0 0 0 0 0 0 0 - - add (M33 a b c d e f g h i) (M33 a' b' c' d' e' f' g' h' i') = M33 - (a + a') (b + b') (c + c') - (d + d') (e + e') (f + f') - (g + g') (h + h') (i + i') - {-# INLINE zero #-} - {-# INLINE add #-} - -instance Eye M33 where - eye = M33 1 0 0 0 1 0 0 0 1 - {-# INLINE eye #-} - -instance Det M33 where det = det3 -instance Inv M33 where inv = inv3 - -instance Mult Double M33 M33 where - mult k (M33 a b c d e f g h i) = M33 - (k * a) (k * b) (k * c) - (k * d) (k * e) (k * f) - (k * g) (k * h) (k * i) - {-# INLINE mult #-} - -instance Mult M33 V3 V3 where - mult (M33 a b c - d e f - g h i) (V3 u v w) = V3 - (a * u + b * v + c * w) - (d * u + e * v + f * w) - (g * u + h * v + i * w) - {-# INLINE mult #-} - --- TODO: instance Mult M33 M33 M33 where - -det3 :: M33 -> Double -det3 (M33 a b c - d e f - g h i) - = a * (e*i-f*h) - d * (b*i-c*h) + g * (b*f-c*e) -{-# INLINE det3 #-} - -inv3 :: M33 -> M33 -inv3 m@(M33 a b c - d e f - g h i) - = M33 a' b' c' - d' e' f' - g' h' i' - where - a' = cofactor e f h i / det - b' = cofactor c b i h / det - c' = cofactor b c e f / det - d' = cofactor f d i g / det - e' = cofactor a c g i / det - f' = cofactor c a f d / det - g' = cofactor d e g h / det - h' = cofactor b a h g / det - i' = cofactor a b d e / det - cofactor q r s t = det2 (M22 q r s t) - det = det3 m -{-# INLINE inv3 #-} - -------------------------------------------------------------------------------- --- Regressions -------------------------------------------------------------------------------- - --- | Linear regression. --- --- The type is --- --- @ --- 'linear' :: [('Double', 'Double')] -> 'V2' --- @ --- --- but overloaded to work with boxed and unboxed 'Vector's. --- --- >>> let input1 = [(0, 1), (1, 3), (2, 5)] --- >>> 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)] --- >>> PP $ linear input2 --- V2 2.0063 0.88685 --- -linear :: (Foldable' xs x, IsDoublePair x) => xs -> V2 -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 - 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 --- --- @ --- 'quadratic' :: [('Double', 'Double')] -> 'V3' --- @ --- --- but overloaded to work with boxed and unboxed 'Vector's. --- --- >>> let input1 = [(0, 1), (1, 3), (2, 5)] --- >>> quadratic input1 --- 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)] --- >>> PP $ quadratic input2 --- V3 (-0.00589) 2.0313 0.87155 --- --- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)] --- >>> 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 _) (V2 _y2 _) = kahan3 data_ - - n :: Double - n = fromIntegral n' - --- | Like 'quadratic' but also return parameters' standard errors. --- --- >>> PP $ quadraticWithErrors input2 --- (V3 (-0.00589) 2.0313 0.87155, V3 0.09281 0.41070 0.37841) --- --- >>> 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 - 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 -------------------------------------------------------------------------------- - --- | Like 'Foldable' but with element in the class definition. -class Foldable' xs x | xs -> x where - foldl' :: (b -> x -> b) -> b -> xs -> b - -instance Foldable' [a] a where foldl' = L.foldl' -instance Foldable' (V.Vector a) a where foldl' = V.foldl' -instance U.Unbox a => Foldable' (U.Vector a) a where foldl' = U.foldl' - --- | Class witnessing that @dp@ has a pair of 'Double's. -class IsDoublePair dp where - withDP :: dp -> (Double -> Double -> r) -> r - makeDP :: Double -> Double -> dp - -instance IsDoublePair V2 where - withDP (V2 x y) k = k x y - makeDP = V2 - -instance (a ~ Double, b ~ Double) => IsDoublePair (a, b) where - withDP ~(x, y) k = k x y - makeDP = (,) - -------------------------------------------------------------------------------- --- Kahan2 -------------------------------------------------------------------------------- - -data Kahan2 = K2 - { k2n :: {-# UNPACK #-} !Int - , k2x :: {-# UNPACK #-} !V2 - , k2x2 :: {-# UNPACK #-} !V2 - , k2y :: {-# UNPACK #-} !V2 - , k2y2 :: {-# UNPACK #-} !V2 - , k2xy :: {-# UNPACK #-} !V2 - } - -zeroKahan2 :: Kahan2 -zeroKahan2 = K2 0 zero zero zero zero zero - --- | https://en.wikipedia.org/wiki/Kahan_summation_algorithm -addKahan :: V2 -> Double -> V2 -addKahan (V2 acc c) i = - let y = i - c - t = acc + y - in V2 t ((t - acc) - y) - -kahan2 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan2 -kahan2 = foldl' f zeroKahan2 where - 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)) - -------------------------------------------------------------------------------- --- Kahan3 -------------------------------------------------------------------------------- - -data Kahan3 = K3 - { 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 zero - -kahan3 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan3 -kahan3 = foldl' f zeroKahan3 where - f (K3 n x x2 x3 x4 y xy x2y y2) uv = withDP uv $ \u v -> - let u2 = u * u - in K3 - (succ n) - (addKahan x u) - (addKahan x2 u2) - (addKahan x3 (u * u2)) - (addKahan x4 (u2 * u2)) - (addKahan y v) - (addKahan xy (u * v)) - (addKahan x2y (u2 * v)) - (addKahan y2 (v * v)) +{-# LANGUAGE RecordWildCards #-} +-- | +-- +-- @regression-simple@ provides (hopefully) simple regression functions. +-- +-- The @'linear' :: Foldable f => (a -> (Double, Double)) -> f a -> 'V2'@ +-- is the simplest one. +-- +-- There are variants with weights, y-errors, and x and y-errors. +-- In addition, package includes Levenberg–Marquardt algorithm implementation +-- to fit arbitrary functions (with one, two or three parameters), +-- as long as you can give their partial derivatives as well (@ad@ package is handy for that). +-- +-- For multiple independent variable ordinary least squares +-- or Levenberg-Marquard with functions with \> 3 parameter you should look elsewhere. +-- +-- Package has been tested to return similar results as @fit@ functionality in @gnuplot@ +-- (L-M doesn't always converge to exactly the same points in parameter space). +-- +module Math.Regression.Simple ( + -- * Linear regression + linear, + linearFit, + linearWithWeights, + linearWithYerrors, + linearWithXYerrors, + -- ** Step-by-step interface + linearFit', + LinRegAcc (..), + zeroLinRegAcc, + addLinReg, + addLinRegW, + -- * Quadratic regression + quadratic, + quadraticFit, + quadraticWithWeights, + quadraticWithYerrors, + quadraticWithXYerrors, + -- ** Step-by-step interface + quadraticFit', + QuadRegAcc (..), + zeroQuadRegAcc, + addQuadReg, + addQuadRegW, + quadRegAccToLin, + -- * Levenberg–Marquardt algorithm + -- ** One parameter + levenbergMarquardt1, + levenbergMarquardt1WithWeights, + levenbergMarquardt1WithYerrors, + levenbergMarquardt1WithXYerrors, + -- ** Two parameters + levenbergMarquardt2, + levenbergMarquardt2WithWeights, + levenbergMarquardt2WithYerrors, + levenbergMarquardt2WithXYerrors, + -- ** Three parameters + levenbergMarquardt3, + levenbergMarquardt3WithWeights, + levenbergMarquardt3WithYerrors, + levenbergMarquardt3WithXYerrors, + -- * Auxiliary types + Fit (..), + V2 (..), + V3 (..), +) where + +import Control.DeepSeq (NFData (..)) + +import qualified Data.Foldable as F +import qualified Data.List.NonEmpty as NE + +import Math.Regression.Simple.LinAlg +import Numeric.KBN + +-- $setup +-- >>> :set -XDeriveFunctor -XDeriveFoldable -XDeriveTraversable +-- +-- >>> import Numeric (showFFloat) +-- >>> import Data.List.NonEmpty (NonEmpty (..)) +-- >>> import qualified Data.List.NonEmpty as NE +-- +-- Don't show too much decimal digits +-- +-- >>> let showDouble x = showFFloat (Just (min 5 (5 - ceiling (logBase 10 x)))) (x :: Double) +-- >>> 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' Int where showsPrec' = showsPrec +-- >>> instance (Show' a, Show' b) => Show' (a, b) where showsPrec' d (x, y) = showParen True $ showsPrec' 0 x . showString ", " . showsPrec' 0 y +-- >>> instance Show' v => Show' (Fit v) where showsPrec' d (Fit p e ndf wssr) = showParen (d > 10) $ showString "Fit " . showsPrec' 11 p . showChar ' ' . showsPrec' 11 e . showChar ' ' . showsPrec' 11 ndf . showChar ' ' . showsPrec' 11 wssr +-- >>> 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 +-- >>> instance Show' a => Show' [a] where showsPrec' _ [] = id; showsPrec' _ [x] = showsPrec' 0 x; showsPrec' _ (x:xs) = showsPrec' 0 x . showChar '\n' . showsPrec' 0 xs +-- >>> instance Show' a => Show' (NonEmpty a) where showsPrec' _ (x :| []) = showsPrec' 0 x; showsPrec' _ (x:|xs) = showsPrec' 0 x . showChar '\n' . showsPrec' 0 xs +-- >>> +-- +-- 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)] +-- +-- >>> let sq z = z * z +-- + +------------------------------------------------------------------------------- +-- Linear +------------------------------------------------------------------------------- + +-- | Linear regression. +-- +-- >>> let input1 = [(0, 1), (1, 3), (2, 5)] +-- >>> PP $ linear id 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)] +-- >>> PP $ linear id input2 +-- V2 2.0063 0.88685 +-- +linear :: F.Foldable f => (a -> (Double, Double)) -> f a -> V2 +linear f = fitParams . linearFit f + +-- | Like 'linear' but returns complete 'Fit'. +-- +-- 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: +-- +-- >>> let fit = linearFit id input1 +-- >>> PP fit +-- Fit (V2 2.0000 1.00000) (V2 0.00000 0.00000) 1 0.00000 +-- +-- The second input is quite good: +-- +-- >>> PP $ linearFit id input2 +-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962 +-- +-- But the third input isn't so much, +-- standard error of a slope parameter is 20%. +-- +-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)] +-- >>> PP $ linearFit id input3 +-- Fit (V2 3.0000 1.00000) (V2 0.63246 1.1832) 2 4.0000 +-- +linearFit :: F.Foldable f => (a -> (Double, Double)) -> f a -> Fit V2 +linearFit f = linearFit' . linRegAcc f + +-- | Weighted linear regression. +-- +-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)] +-- >>> PP $ linearFit id input2 +-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962 +-- +-- >>> let input2w = [(0.1, 1.2, 1), (1.3, 3.1, 1), (1.9, 4.9, 1), (3.0, 7.1, 1/4), (4.1, 9.0, 1/4)] +-- >>> PP $ linearWithWeights id input2w +-- Fit (V2 2.0060 0.86993) (V2 0.12926 0.23696) 3 0.22074 +-- +linearWithWeights :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V2 +linearWithWeights f = linearFit' . linRegAccW f + +-- | Linear regression with y-errors. +-- +-- >>> let input2y = [(0.1, 1.2, 0.12), (1.3, 3.1, 0.31), (1.9, 4.9, 0.49), (3.0, 7.1, 0.71), (4.1, 9.0, 1.9)] +-- >>> let fit = linearWithYerrors id input2y +-- >>> PP fit +-- Fit (V2 1.9104 0.98302) (V2 0.13006 0.10462) 3 2.0930 +-- +-- When we know actual y-errors, we can calculate the Q-value using @statistics@ package: +-- +-- >>> import qualified Statistics.Distribution as S +-- >>> import qualified Statistics.Distribution.ChiSquared as S +-- >>> S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit) +-- 0.446669639443138 +-- +-- or using @math-functions@ +-- +-- >>> import Numeric.SpecFunctions (incompleteGamma) +-- >>> incompleteGamma (fromIntegral (fitNDF fit) / 2) (fitWSSR fit / 2) +-- 0.446669639443138 +-- +-- It is not uncommon to deem acceptable on equal terms any models with, say, Q > 0.001. +-- If Q is too large, too near to 1 is most likely caused by overestimating +-- the y-errors. +-- +-- +linearWithYerrors :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V2 +linearWithYerrors f = linearWithWeights f' where + f' a = case f a of + (x, y, dy) -> (x, y, recip (dy * dy)) + +-- | Iterative linear regression with x and y errors. +-- +-- /Orear, J. (1982). Least squares when both variables have uncertainties. American Journal of Physics, 50(10), 912–916. doi:10.1119\/1.12972/ +-- +-- >>> let input2xy = [(0.1, 1.2, 0.01, 0.12), (1.3, 3.1, 0.13, 0.31), (1.9, 4.9, 0.19, 0.49), (3.0, 7.1, 0.3, 0.71), (4.1, 9.0, 0.41, 1.9)] +-- >>> let fit :| fits = linearWithXYerrors id input2xy +-- +-- First fit is done using 'linearWithYerrors': +-- +-- >>> PP fit +-- Fit (V2 1.9104 0.98302) (V2 0.13006 0.10462) 3 2.0930 +-- +-- After that the effective variance is used to refine the fit, +-- just a few iterations is often enough: +-- +-- >>> PP $ take 3 fits +-- Fit (V2 1.9092 0.99251) (V2 0.12417 0.08412) 3 1.2992 +-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998 +-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998 +-- +linearWithXYerrors + :: F.Foldable f + => (a -> (Double, Double, Double, Double)) -- ^ \(x_i, y_i, \delta x_i, \delta y_i\) + -> f a -- ^ data + -> NE.NonEmpty (Fit V2) +linearWithXYerrors f xs = iterate1 go fit0 where + fit0 = linearWithYerrors (\a -> case f a of (x,y,_,dy) -> (x,y,dy)) xs + go fit = linearWithWeights (\a -> case f a of (x,y,dx,dy) -> (x,y,recip $ sq (param1 * dx) + sq dy)) xs where + V2 param1 _ = fitParams fit + +-- >>> import qualified Numeric.AD.Mode.Reverse.Double as AD +-- >>> data H3 a = H3 a a a deriving (Functor, Foldable, Traversable) +-- >>> let linearF (H3 a b x) = a * x + b +-- >>> let lin' (V2 a b) (x, y, dx, dy) = case AD.grad' linearF (H3 a b x) of (f, H3 da db f') -> (y, f, V2 da db, recip $ sq (f' * dx) + sq dy) +-- +-- >>> PP $ NE.last $ levenbergMarquardt2WithWeights lin' (V2 1 1) input2xy +-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998 + +-- | Calculate linear fit from 'LinRegAcc'. +linearFit' :: LinRegAcc -> Fit V2 +linearFit' LinRegAcc {..} = Fit params errors ndf wssr where + matrix@(SM22 a11 _ a22) = inv (SM22 x2 x w) + 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. + wssr = max 0 (y2 - a * xy - b * y) + ndf = lra_n - 2 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + sb = sqrt (a22 * wssr / ndf') + + w = getKBN lra_w + x = getKBN lra_x + x2 = getKBN lra_x2 + y = getKBN lra_y + xy = getKBN lra_xy + y2 = getKBN lra_y2 + + -- is it useful? + -- r2 = (n * xy - x*y) / (sqrt (n * x2 - x*x) * sqrt (n * y2 - y*y)) + +linRegAcc :: F.Foldable f => (a -> (Double, Double)) -> f a -> LinRegAcc +linRegAcc f = F.foldl' (\acc a -> case f a of (x,y) -> addLinReg acc x y) zeroLinRegAcc + +linRegAccW :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> LinRegAcc +linRegAccW f = F.foldl' (\acc a -> case f a of (x,y,w) -> addLinRegW acc x y w) zeroLinRegAcc + +------------------------------------------------------------------------------- +-- Quadractic +------------------------------------------------------------------------------- + +-- | Quadratic regression. +-- +-- >>> let input1 = [(0, 1), (1, 3), (2, 5)] +-- >>> quadratic id input1 +-- 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)] +-- >>> PP $ quadratic id input2 +-- V3 (-0.00589) 2.0313 0.87155 +-- +-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)] +-- >>> PP $ quadratic id input3 +-- V3 1.00000 0.00000 2.0000 +-- +quadratic :: F.Foldable f => (a -> (Double, Double)) -> f a -> V3 +quadratic f = fitParams . quadraticFit f + +-- | Like 'quadratic' but returns complete 'Fit'. +-- +-- >>> PP $ quadraticFit id input2 +-- Fit (V3 (-0.00589) 2.0313 0.87155) (V3 0.09281 0.41070 0.37841) 2 0.25910 +-- +-- >>> PP $ quadraticFit id input3 +-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000 +-- +quadraticFit :: F.Foldable f => (a -> (Double, Double)) -> f a -> Fit V3 +quadraticFit f = quadraticFit' . quadRegAcc f + +-- | Weighted quadratic regression. +-- +-- >>> let input2w = [(0.1, 1.2, 1), (1.3, 3.1, 1), (1.9, 4.9, 1), (3.0, 7.1, 1/4), (4.1, 9.0, 1/4)] +-- >>> PP $ quadraticWithWeights id input2w +-- Fit (V3 0.02524 1.9144 0.91792) (V3 0.10775 0.42106 0.35207) 2 0.21484 +-- +quadraticWithWeights :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V3 +quadraticWithWeights f = quadraticFit' . quadRegAccW f + +-- | Quadratic regression with y-errors. +-- +-- >>> let input2y = [(0.1, 1.2, 0.12), (1.3, 3.1, 0.31), (1.9, 4.9, 0.49), (3.0, 7.1, 0.71), (4.1, 9.0, 0.9)] +-- >>> PP $ quadraticWithYerrors id input2y +-- Fit (V3 0.08776 1.6667 1.0228) (V3 0.10131 0.31829 0.11917) 2 1.5398 +-- +quadraticWithYerrors :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V3 +quadraticWithYerrors f = quadraticWithWeights f' where + f' a = case f a of + (x, y, dy) -> (x, y, recip (dy * dy)) + +-- | Iterative quadratic regression with x and y errors. +-- +-- /Orear, J. (1982). Least squares when both variables have uncertainties. American Journal of Physics, 50(10), 912–916. doi:10.1119\/1.12972/ +-- +quadraticWithXYerrors + :: F.Foldable f + => (a -> (Double, Double, Double, Double)) -- ^ \(x_i, y_i, \delta x_i, \delta y_i\) + -> f a -- ^ data + -> NE.NonEmpty (Fit V3) +quadraticWithXYerrors f xs = iterate1 go fit0 where + fit0 = quadraticWithYerrors (\a -> case f a of (x,y,_,dy) -> (x,y,dy)) xs + go fit = quadraticWithWeights (\a -> case f a of (x,y,dx,dy) -> (x,y,recip $ sq ((2 * p1 * x + p2) * dx) + sq dy)) xs where + V3 p1 p2 _ = fitParams fit + +-- | Calculate quadratic fit from 'QuadRegAcc'. +quadraticFit' :: QuadRegAcc -> Fit V3 +quadraticFit' QuadRegAcc {..} = Fit params errors ndf wssr where + matrix@(SM33 a11 + _ a22 + _ _ a33) = inv (SM33 x4 + x3 x2 + x2 x w) + + params@(V3 a b c) = mult matrix (V3 x2y xy y) + + errors = V3 sa sb sc + + wssr = max 0 (y2 - a * x2y - b * xy - c * y) + ndf = qra_n - 3 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + sb = sqrt (a22 * wssr / ndf') + sc = sqrt (a33 * wssr / ndf') + + w = getKBN qra_w + x = getKBN qra_x + x2 = getKBN qra_x2 + x3 = getKBN qra_x3 + x4 = getKBN qra_x4 + y = getKBN qra_y + xy = getKBN qra_xy + x2y = getKBN qra_x2y + y2 = getKBN qra_y2 + + -- is it useful? + -- r2 = (n * xy - x*y) / (sqrt (n * x2 - x*x) * sqrt (n * y2 - y*y)) + +quadRegAcc :: F.Foldable f => (a -> (Double, Double)) -> f a -> QuadRegAcc +quadRegAcc f = F.foldl' (\acc a -> case f a of (x,y) -> addQuadReg acc x y) zeroQuadRegAcc + +quadRegAccW :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> QuadRegAcc +quadRegAccW f = F.foldl' (\acc a -> case f a of (x,y,w) -> addQuadRegW acc x y w) zeroQuadRegAcc + +------------------------------------------------------------------------------- +-- Levenberg–Marquardt 1 +------------------------------------------------------------------------------- + +-- | Levenberg–Marquardt for functions with one parameter. +-- +-- See 'levenbergMarquardt2' for examples, this is very similar. +-- +-- For example we can fit \(f = x \mapsto \beta x + 1\), its derivative is \(\partial_\beta f = x \mapsto x\). +-- +-- >>> let scale a (x, y) = (y, a * x + 1, x) +-- >>> PP $ NE.last $ levenbergMarquardt1 scale 1 input2 +-- Fit 1.9685 0.04735 4 0.27914 +-- +-- Not bad, but worse then linear fit which fits the intercept point too. +-- +levenbergMarquardt1 + :: F.Foldable f + => (Double -> a -> (Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x)\) + -> Double -- ^ initial parameter, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit Double) -- ^ non-empty list of iteration results +levenbergMarquardt1 f b0 xs = loop lambda0 b0 acc0 where + acc0 = calcAcc b0 + + lambda0 = c11 + where + c11 = getKBN $ lm1_c11 acc0 + + calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d) -> addLM1Acc acc y g d) zeroLM1Acc xs + + loop lambda beta acc + | lmStop lambda wssr wssr' + = Fit beta errors ndf wssr NE.:| [] + + | wssr' >= wssr + = loop (lambda * 10) beta acc + + | otherwise + = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc' + + where + lambda1 z = (1 + lambda) * z + matrix = inv (lambda1 c11) + delta = mult matrix z1 + + beta' = add beta delta + acc' = calcAcc beta' + + a11 = inv c11 + errors = sa + + wssr = max 0 $ getKBN $ lm1_wssr acc + wssr' = getKBN $ lm1_wssr acc' + + ndf = lm1_n acc - 1 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + + c11 = getKBN $ lm1_c11 acc + z1 = getKBN $ lm1_z1 acc + +-- | 'levenbergMarquardt1' with weights. +levenbergMarquardt1WithWeights + :: F.Foldable f + => (Double -> a -> (Double, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), w_i\) + -> Double -- ^ initial parameter, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit Double) -- ^ non-empty list of iteration results +levenbergMarquardt1WithWeights f b0 xs = loop lambda0 b0 acc0 where + acc0 = calcAcc b0 + + lambda0 = c11 + where + c11 = getKBN $ lm1_c11 acc0 + + calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d, w) -> addLM1AccW acc y g d w) zeroLM1Acc xs + + loop lambda beta acc + | lmStop lambda wssr wssr' + = Fit beta errors ndf wssr NE.:| [] + + | wssr' >= wssr + = loop (lambda * 10) beta acc + + | otherwise + = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc' + + where + lambda1 z = (1 + lambda) * z + matrix = inv (lambda1 c11) + delta = mult matrix z1 + + beta' = add beta delta + acc' = calcAcc beta' + + a11 = inv c11 + errors = sa + + wssr = max 0 $ getKBN $ lm1_wssr acc + wssr' = getKBN $ lm1_wssr acc' + + ndf = lm1_n acc - 1 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + + c11 = getKBN $ lm1_c11 acc + z1 = getKBN $ lm1_z1 acc + +-- | 'levenbergMarquardt1' with Y-errors. +levenbergMarquardt1WithYerrors + :: F.Foldable f + => (Double -> a -> (Double, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), \delta y_i\) + -> Double -- ^ initial parameter, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit Double) -- ^ non-empty list of iteration results +levenbergMarquardt1WithYerrors f = levenbergMarquardt1WithWeights f' where + f' beta x = case f beta x of (y, fbetax, grad, dy) -> (y, fbetax, grad, recip $ sq dy) + +-- | 'levenbergMarquardt1' with XY-errors. +levenbergMarquardt1WithXYerrors + :: F.Foldable f + => (Double -> a -> (Double, Double, Double, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\) + -> Double -- ^ initial parameter, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit Double) -- ^ non-empty list of iteration results +levenbergMarquardt1WithXYerrors g = levenbergMarquardt1WithWeights g' where + g' beta x = case g beta x of (y, fbetax, grad, f', dx, dy) -> (y, fbetax, grad, recip $ sq (f' * dx) + sq dy) + +data LM1Acc = LM1Acc + { lm1_n :: !Int + , lm1_c11 :: !KBN + , lm1_z1 :: !KBN + , lm1_wssr :: !KBN + } + deriving Show + +zeroLM1Acc :: LM1Acc +zeroLM1Acc = LM1Acc 0 zeroKBN zeroKBN zeroKBN + +addLM1Acc :: LM1Acc -> Double -> Double -> Double -> LM1Acc +addLM1Acc LM1Acc {..} y f d1 = LM1Acc + { lm1_n = lm1_n + 1 + , lm1_c11 = addKBN lm1_c11 (d1 * d1) + , lm1_z1 = addKBN lm1_z1 (d1 * res) + , lm1_wssr = addKBN lm1_wssr (res * res) + } + where + res = y - f + +addLM1AccW :: LM1Acc -> Double -> Double -> Double -> Double -> LM1Acc +addLM1AccW LM1Acc {..} y f d1 w = LM1Acc + { lm1_n = lm1_n + 1 + , lm1_c11 = addKBN lm1_c11 (w * d1 * d1) + , lm1_z1 = addKBN lm1_z1 (w * d1 * res) + , lm1_wssr = addKBN lm1_wssr (w * res * res) + } + where + res = y - f + +------------------------------------------------------------------------------- +-- Levenberg–Marquardt 2 +------------------------------------------------------------------------------- + +-- | Levenberg–Marquardt for functions with two parameters. +-- +-- You can use this sledgehammer to do a a linear fit: +-- +-- >>> let lin (V2 a b) (x, y) = (y, a * x + b, V2 x 1) +-- +-- We can then use 'levenbergMarquardt2' to find a fit: +-- +-- >>> PP $ levenbergMarquardt2 lin (V2 1 1) input2 +-- Fit (V2 1.00000 1.00000) (V2 1.0175 2.5385) 3 29.470 +-- Fit (V2 1.0181 1.0368) (V2 0.98615 2.4602) 3 27.681 +-- Fit (V2 1.1557 1.2988) (V2 0.75758 1.8900) 3 16.336 +-- Fit (V2 1.5463 1.6577) (V2 0.29278 0.73043) 3 2.4400 +-- Fit (V2 1.9129 1.1096) (V2 0.11033 0.27524) 3 0.34645 +-- Fit (V2 2.0036 0.89372) (V2 0.09552 0.23830) 3 0.25970 +-- Fit (V2 2.0063 0.88687) (V2 0.09550 0.23826) 3 0.25962 +-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962 +-- +-- This is the same result what 'linearFit' returns: +-- +-- >>> PP $ linearFit id input2 +-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962 +-- +-- == Using AD +-- +-- You can use @ad@ to calculate derivatives for you. +-- +-- >>> import qualified Numeric.AD.Mode.Reverse.Double as AD +-- +-- We need a ('Traversable') homogenic triple to represent the two parameters and @x@: +-- +-- >>> data H3 a = H3 a a a deriving (Functor, Foldable, Traversable) +-- +-- Then we define a function @ad@ can operate with: +-- +-- >>> let linearF (H3 a b x) = a * x + b +-- +-- which we can use to fit the curve in generic way: +-- +-- >>> let lin' (V2 a b) (x, y) = case AD.grad' linearF (H3 a b x) of (f, H3 da db _f') -> (y, f, V2 da db) +-- >>> PP $ levenbergMarquardt2 lin' (V2 1 1) input2 +-- Fit (V2 1.00000 1.00000) (V2 1.0175 2.5385) 3 29.470 +-- Fit (V2 1.0181 1.0368) (V2 0.98615 2.4602) 3 27.681 +-- Fit (V2 1.1557 1.2988) (V2 0.75758 1.8900) 3 16.336 +-- Fit (V2 1.5463 1.6577) (V2 0.29278 0.73043) 3 2.4400 +-- Fit (V2 1.9129 1.1096) (V2 0.11033 0.27524) 3 0.34645 +-- Fit (V2 2.0036 0.89372) (V2 0.09552 0.23830) 3 0.25970 +-- Fit (V2 2.0063 0.88687) (V2 0.09550 0.23826) 3 0.25962 +-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962 +-- +-- == Non-polynomial example +-- +-- We can fit other curves too, for example an example from Wikipedia +-- https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm#Example +-- +-- >>> let rateF (H3 vmax km s) = (vmax * s) / (km + s) +-- >>> let rateF' (V2 vmax km) (x, y) = case AD.grad' rateF (H3 vmax km x) of (f, H3 vmax' km' _) -> (y, f, V2 vmax' km') +-- >>> let input = zip [0.038,0.194,0.425,0.626,1.253,2.500,3.740] [0.050,0.127,0.094,0.2122,0.2729,0.2665,0.3317] +-- >>> PP $ levenbergMarquardt2 rateF' (V2 0.9 0.2) input +-- Fit (V2 0.90000 0.20000) (V2 0.43304 0.43936) 5 1.4455 +-- Fit (V2 0.83306 0.25278) (V2 0.39164 0.49729) 5 1.0055 +-- Fit (V2 0.59437 0.43508) (V2 0.21158 0.53403) 5 0.18832 +-- Fit (V2 0.39687 0.56324) (V2 0.05723 0.25666) 5 0.01062 +-- Fit (V2 0.36289 0.56104) (V2 0.04908 0.24007) 5 0.00784 +-- Fit (V2 0.36190 0.55662) (V2 0.04887 0.23843) 5 0.00784 +-- Fit (V2 0.36184 0.55629) (V2 0.04885 0.23830) 5 0.00784 +-- Fit (V2 0.36184 0.55627) (V2 0.04885 0.23829) 5 0.00784 +-- +-- We get the same result as in the article: 0.362 and 0.556 +-- +-- The algorithm terminates when a scaling parameter \(\lambda\) becomes larger than 1e20 or smaller than 1e-20, or relative WSSR change is smaller than 1e-10, or sum-of-squared-residuals candidate becomes @NaN@ (i.e. when it would start to produce garbage). +-- You may want to terminate sooner, Numerical Recipes suggest to stop when WSSR decreases by a neglible amount absolutely or fractionally. +-- +levenbergMarquardt2 + :: F.Foldable f + => (V2 -> a -> (Double, Double, V2)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x)\) + -> V2 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V2) -- ^ non-empty list of iteration results +levenbergMarquardt2 f b0 xs = loop lambda0 b0 acc0 where + acc0 = calcAcc b0 + + calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d) -> addLM2Acc acc y g d) zeroLM2Acc xs + + lambda0 = max l1 l2 + where + V2 l1 l2 = eigenSM22 c11 c12 c22 + c11 = getKBN $ lm2_c11 acc0 + c12 = getKBN $ lm2_c11 acc0 + c22 = getKBN $ lm2_c22 acc0 + + loop lambda beta acc + | lmStop lambda wssr wssr' + = Fit beta errors ndf wssr NE.:| [] + + | wssr' >= wssr + = loop (lambda * 10) beta acc + + | otherwise + = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc' + + where + lambda1 z = (1 + lambda) * z + matrix = inv (SM22 (lambda1 c11) c12 (lambda1 c22)) + delta = mult matrix (V2 z1 z2) + + beta' = add beta delta + acc' = calcAcc beta' + + SM22 a11 _ a22 = inv (SM22 c11 c12 c22) + errors = V2 sa sb + + wssr = max 0 $ getKBN $ lm2_wssr acc + wssr' = getKBN $ lm2_wssr acc' + + ndf = lm2_n acc - 2 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + sb = sqrt (a22 * wssr / ndf') + + c11 = getKBN $ lm2_c11 acc + c12 = getKBN $ lm2_c12 acc + c22 = getKBN $ lm2_c22 acc + z1 = getKBN $ lm2_z1 acc + z2 = getKBN $ lm2_z2 acc + +-- | 'levenbergMarquardt2' with weights. +-- +-- Because 'levenbergMarquardt2' is an iterative algorithm, +-- not only we can use it to fit curves with known y-errors ('levenbergMarquardt2WithYerrors'), +-- but also with both x and y-errors ('levenbergMarquardt2WithXYerrors'). +-- +levenbergMarquardt2WithWeights + :: F.Foldable f + => (V2 -> a -> (Double, Double, V2, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), w_i\) + -> V2 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V2) -- ^ non-empty list of iteration results +levenbergMarquardt2WithWeights f b0 xs = loop lambda0 b0 acc0 where + acc0 = calcAcc b0 + + lambda0 = max l1 l2 + where + V2 l1 l2 = eigenSM22 c11 c12 c22 + c11 = getKBN $ lm2_c11 acc0 + c12 = getKBN $ lm2_c11 acc0 + c22 = getKBN $ lm2_c22 acc0 + + calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d, w) -> addLM2AccW acc y g d w) zeroLM2Acc xs + + loop lambda beta acc + | lmStop lambda wssr wssr' + = Fit beta errors ndf wssr NE.:| [] + + | wssr' >= wssr + = loop (lambda * 10) beta acc + + | otherwise + = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc' + + where + lambda1 z = (1 + lambda) * z + matrix = inv (SM22 (lambda1 c11) c12 (lambda1 c22)) + delta = mult matrix (V2 z1 z2) + + beta' = add beta delta + acc' = calcAcc beta' + + SM22 a11 _ a22 = inv (SM22 c11 c12 c22) + errors = V2 sa sb + + wssr = max 0 $ getKBN $ lm2_wssr acc + wssr' = getKBN $ lm2_wssr acc' + + ndf = lm2_n acc - 2 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + sb = sqrt (a22 * wssr / ndf') + + c11 = getKBN $ lm2_c11 acc + c12 = getKBN $ lm2_c12 acc + c22 = getKBN $ lm2_c22 acc + z1 = getKBN $ lm2_z1 acc + z2 = getKBN $ lm2_z2 acc + +-- | 'levenbergMarquardt2' with Y-errors. +levenbergMarquardt2WithYerrors + :: F.Foldable f + => (V2 -> a -> (Double, Double, V2, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \delta y_i\) + -> V2 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V2) -- ^ non-empty list of iteration results +levenbergMarquardt2WithYerrors f = levenbergMarquardt2WithWeights f' where + f' beta x = case f beta x of (y, fbetax, grad, dy) -> (y, fbetax, grad, recip $ sq dy) + +-- | 'levenbergMarquardt2' with XY-errors. +levenbergMarquardt2WithXYerrors + :: F.Foldable f + => (V2 -> a -> (Double, Double, V2, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\) + -> V2 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V2) -- ^ non-empty list of iteration results +levenbergMarquardt2WithXYerrors g = levenbergMarquardt2WithWeights g' where + g' beta x = case g beta x of (y, fbetax, grad, f', dx, dy) -> (y, fbetax, grad, recip $ sq (f' * dx) + sq dy) + +data LM2Acc = LM2Acc + { lm2_n :: !Int + , lm2_c11 :: !KBN + , lm2_c12 :: !KBN + , lm2_c22 :: !KBN + , lm2_z1 :: !KBN + , lm2_z2 :: !KBN + , lm2_wssr :: !KBN + } + deriving Show + +zeroLM2Acc :: LM2Acc +zeroLM2Acc = LM2Acc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN + +addLM2Acc :: LM2Acc -> Double -> Double -> V2 -> LM2Acc +addLM2Acc LM2Acc {..} y f (V2 d1 d2) = LM2Acc + { lm2_n = lm2_n + 1 + , lm2_c11 = addKBN lm2_c11 (d1 * d1) + , lm2_c12 = addKBN lm2_c12 (d1 * d2) + , lm2_c22 = addKBN lm2_c22 (d2 * d2) + , lm2_z1 = addKBN lm2_z1 (d1 * res) + , lm2_z2 = addKBN lm2_z2 (d2 * res) + , lm2_wssr = addKBN lm2_wssr (res * res) + } + where + res = y - f + +addLM2AccW :: LM2Acc -> Double -> Double -> V2 -> Double -> LM2Acc +addLM2AccW LM2Acc {..} y f (V2 d1 d2) w = LM2Acc + { lm2_n = lm2_n + 1 + , lm2_c11 = addKBN lm2_c11 (w * d1 * d1) + , lm2_c12 = addKBN lm2_c12 (w * d1 * d2) + , lm2_c22 = addKBN lm2_c22 (w * d2 * d2) + , lm2_z1 = addKBN lm2_z1 (w * d1 * res) + , lm2_z2 = addKBN lm2_z2 (w * d2 * res) + , lm2_wssr = addKBN lm2_wssr (w * res * res) + } + where + res = y - f + +------------------------------------------------------------------------------- +-- Levenberg–Marquardt 3 +------------------------------------------------------------------------------- + +-- | Levenberg–Marquardt for functions with three parameters. +-- +-- See 'levenbergMarquardt2' for examples, this is very similar. +-- +-- >>> let quad (V3 a b c) (x, y) = (y, a * x * x + b * x + c, V3 (x * x) x 1) +-- >>> PP $ NE.last $ levenbergMarquardt3 quad (V3 2 2 2) input3 +-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000 +-- +-- Same as quadratic fit, just less direct: +-- +-- >>> PP $ quadraticFit id input3 +-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000 +-- +levenbergMarquardt3 + :: F.Foldable f + => (V3 -> a -> (Double, Double, V3)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x)\) + -> V3 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V3) -- ^ non-empty list of iteration results +levenbergMarquardt3 f b0 xs = loop lambda0 b0 acc0 where + acc0 = calcAcc b0 + + calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d) -> addLM3Acc acc y g d) zeroLM3Acc xs + + -- frobenius norm is larger than largest eigen value. + -- calculating the eigen values for 3x3 (symmetric) matrix is becoming complicated. + lambda0 = sqrt $ sumKBN [ sq c11 + , 2 * sq c12, sq c22 + , 2 * sq c13, 2 * sq c23, sq c33] + where + c11 = getKBN $ lm3_c11 acc0 + c12 = getKBN $ lm3_c12 acc0 + c13 = getKBN $ lm3_c13 acc0 + c22 = getKBN $ lm3_c22 acc0 + c23 = getKBN $ lm3_c23 acc0 + c33 = getKBN $ lm3_c33 acc0 + + loop lambda beta acc + | lmStop lambda wssr wssr' + = Fit beta errors ndf wssr NE.:| [] + + | wssr' >= wssr + = loop (lambda * 10) beta acc + + | otherwise + = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc' + + where + lambda1 z = (1 + lambda) * z + matrix = inv (SM33 (lambda1 c11) + c12 (lambda1 c22) + c13 c23 (lambda1 c33)) + delta = mult matrix (V3 z1 z2 z3) + + beta' = add beta delta + acc' = calcAcc beta' + + SM33 a11 + _ a22 + _ _ a33 = inv (SM33 c11 + c12 c22 + c13 c23 c33) + errors = V3 sa sb sc + + wssr = max 0 $ getKBN $ lm3_wssr acc + wssr' = getKBN $ lm3_wssr acc' + + ndf = lm3_n acc - 3 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + sb = sqrt (a22 * wssr / ndf') + sc = sqrt (a33 * wssr / ndf') + + c11 = getKBN $ lm3_c11 acc + c12 = getKBN $ lm3_c12 acc + c13 = getKBN $ lm3_c13 acc + c22 = getKBN $ lm3_c22 acc + c23 = getKBN $ lm3_c23 acc + c33 = getKBN $ lm3_c33 acc + z1 = getKBN $ lm3_z1 acc + z2 = getKBN $ lm3_z2 acc + z3 = getKBN $ lm3_z3 acc + +-- | 'levenbergMarquardt3' with weights. +levenbergMarquardt3WithWeights + :: F.Foldable f + => (V3 -> a -> (Double, Double, V3, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), w_i\) + -> V3 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V3) -- ^ non-empty list of iteration results +levenbergMarquardt3WithWeights f b0 xs = loop lambda0 b0 acc0 where + acc0 = calcAcc b0 + + lambda0 = sqrt $ sumKBN [ sq c11 + , 2 * sq c12, sq c22 + , 2 * sq c13, 2 * sq c23, sq c33] + where + c11 = getKBN $ lm3_c11 acc0 + c12 = getKBN $ lm3_c12 acc0 + c13 = getKBN $ lm3_c13 acc0 + c22 = getKBN $ lm3_c22 acc0 + c23 = getKBN $ lm3_c23 acc0 + c33 = getKBN $ lm3_c33 acc0 + + calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d, w) -> addLM3AccW acc y g d w) zeroLM3Acc xs + + loop lambda beta acc + | lmStop lambda wssr wssr' + = Fit beta errors ndf wssr NE.:| [] + + | wssr' >= wssr + = loop (lambda * 10) beta acc + + | otherwise + = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc' + + where + lambda1 z = (1 + lambda) * z + matrix = inv (SM33 (lambda1 c11) + c12 (lambda1 c22) + c13 c23 (lambda1 c33)) + delta = mult matrix (V3 z1 z2 z3) + + beta' = add beta delta + acc' = calcAcc beta' + + SM33 a11 _ _ a22 _ a33 = inv (SM33 c11 c12 c13 c22 c23 c33) + errors = V3 sa sb sc + + wssr = max 0 $ getKBN $ lm3_wssr acc + wssr' = getKBN $ lm3_wssr acc' + + ndf = lm3_n acc - 3 + ndf' = fromIntegral ndf :: Double + + sa = sqrt (a11 * wssr / ndf') + sb = sqrt (a22 * wssr / ndf') + sc = sqrt (a33 * wssr / ndf') + + c11 = getKBN $ lm3_c11 acc + c12 = getKBN $ lm3_c12 acc + c13 = getKBN $ lm3_c13 acc + c22 = getKBN $ lm3_c22 acc + c23 = getKBN $ lm3_c23 acc + c33 = getKBN $ lm3_c33 acc + z1 = getKBN $ lm3_z1 acc + z2 = getKBN $ lm3_z2 acc + z3 = getKBN $ lm3_z3 acc + +-- | 'levenbergMarquardt3' with Y-errors. +levenbergMarquardt3WithYerrors + :: F.Foldable f + => (V3 -> a -> (Double, Double, V3, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \delta y_i\) + -> V3 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V3) -- ^ non-empty list of iteration results +levenbergMarquardt3WithYerrors f = levenbergMarquardt3WithWeights f' where + f' beta x = case f beta x of (y, fbetax, grad, dy) -> (y, fbetax, grad, recip $ sq dy) + +-- | 'levenbergMarquardt3' with XY-errors. +levenbergMarquardt3WithXYerrors + :: F.Foldable f + => (V3 -> a -> (Double, Double, V3, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\) + -> V3 -- ^ initial parameters, \(\beta_0\) + -> f a -- ^ data, \(d\) + -> NE.NonEmpty (Fit V3) -- ^ non-empty list of iteration results +levenbergMarquardt3WithXYerrors g = levenbergMarquardt3WithWeights g' where + g' beta x = case g beta x of (y, fbetax, grad, f', dx, dy) -> (y, fbetax, grad, recip $ sq (f' * dx) + sq dy) + +data LM3Acc = LM3Acc + { lm3_n :: !Int + , lm3_c11 :: !KBN + , lm3_c12 :: !KBN + , lm3_c13 :: !KBN + , lm3_c22 :: !KBN + , lm3_c23 :: !KBN + , lm3_c33 :: !KBN + , lm3_z1 :: !KBN + , lm3_z2 :: !KBN + , lm3_z3 :: !KBN + , lm3_wssr :: !KBN + } + deriving Show + +zeroLM3Acc :: LM3Acc +zeroLM3Acc = LM3Acc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN + +addLM3Acc :: LM3Acc -> Double -> Double -> V3 -> LM3Acc +addLM3Acc LM3Acc {..} y f (V3 d1 d2 d3) = LM3Acc + { lm3_n = lm3_n + 1 + , lm3_c11 = addKBN lm3_c11 (d1 * d1) + , lm3_c12 = addKBN lm3_c12 (d1 * d2) + , lm3_c13 = addKBN lm3_c12 (d1 * d3) + , lm3_c22 = addKBN lm3_c22 (d2 * d2) + , lm3_c23 = addKBN lm3_c22 (d2 * d3) + , lm3_c33 = addKBN lm3_c22 (d3 * d3) + , lm3_z1 = addKBN lm3_z1 (d1 * res) + , lm3_z2 = addKBN lm3_z2 (d2 * res) + , lm3_z3 = addKBN lm3_z3 (d3 * res) + , lm3_wssr = addKBN lm3_wssr (res * res) + } + where + res = y - f + +addLM3AccW :: LM3Acc -> Double -> Double -> V3 -> Double -> LM3Acc +addLM3AccW LM3Acc {..} y f (V3 d1 d2 d3) w = LM3Acc + { lm3_n = lm3_n + 1 + , lm3_c11 = addKBN lm3_c11 (w * d1 * d1) + , lm3_c12 = addKBN lm3_c12 (w * d1 * d2) + , lm3_c13 = addKBN lm3_c12 (w * d1 * d3) + , lm3_c22 = addKBN lm3_c22 (w * d2 * d2) + , lm3_c23 = addKBN lm3_c22 (w * d2 * d3) + , lm3_c33 = addKBN lm3_c22 (w * d3 * d3) + , lm3_z1 = addKBN lm3_z1 (w * d1 * res) + , lm3_z2 = addKBN lm3_z2 (w * d2 * res) + , lm3_z3 = addKBN lm3_z3 (w * d3 * res) + , lm3_wssr = addKBN lm3_wssr (w * res * res) + } + where + res = y - f + +------------------------------------------------------------------------------- +-- Output +------------------------------------------------------------------------------- + +-- | Result of a curve fit. +data Fit v = Fit + { fitParams :: !v -- ^ fit parameters + , fitErrors :: !v -- ^ asympotic standard errors, /assuming a good fit/ + , fitNDF :: !Int -- ^ number of degrees of freedom + , fitWSSR :: !Double -- ^ sum of squares of residuals + } + deriving Show + +instance (NFData v) => NFData (Fit v) where + rnf (Fit p e _ _) = rnf p `seq` rnf e + +------------------------------------------------------------------------------- +-- LinRegAcc +------------------------------------------------------------------------------- + +-- | Linear regression accumulator. +data LinRegAcc = LinRegAcc + { lra_n :: {-# UNPACK #-} !Int -- ^ \(n\) + , lra_w :: {-# UNPACK #-} !KBN -- ^ \(\sum w_i\) + , lra_x :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i \) + , lra_x2 :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i^2 \) + , lra_y :: {-# UNPACK #-} !KBN -- ^ \(\sum y_i \) + , lra_xy :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i y_i \) + , lra_y2 :: {-# UNPACK #-} !KBN -- ^ \(\sum y_i^2 \) + } + deriving Show + +instance NFData LinRegAcc where + rnf LinRegAcc {} = () + +-- | All-zeroes 'LinRegAcc'. +zeroLinRegAcc :: LinRegAcc +zeroLinRegAcc = LinRegAcc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN + +-- | Add a point to linreg accumulator. +addLinReg + :: LinRegAcc + -> Double -- ^ x + -> Double -- ^ y + -> LinRegAcc +addLinReg LinRegAcc {..} x y = LinRegAcc + { lra_n = lra_n + 1 + , lra_w = addKBN lra_w 1 + , lra_x = addKBN lra_x x + , lra_x2 = addKBN lra_x2 (x * x) + , lra_y = addKBN lra_y y + , lra_xy = addKBN lra_xy (x * y) + , lra_y2 = addKBN lra_y2 (y * y) + } + +-- | Add a weighted point to linreg accumulator. +addLinRegW + :: LinRegAcc + -> Double -- ^ x + -> Double -- ^ y + -> Double -- ^ w + -> LinRegAcc +addLinRegW LinRegAcc {..} x y w = LinRegAcc + { lra_n = lra_n + 1 + , lra_w = addKBN lra_w w + , lra_x = addKBN lra_x (w * x) + , lra_x2 = addKBN lra_x2 (w * x * x) + , lra_y = addKBN lra_y (w * y) + , lra_xy = addKBN lra_xy (w * x * y) + , lra_y2 = addKBN lra_y2 (w * y * y) + } + +------------------------------------------------------------------------------- +-- QuadRegAcc +------------------------------------------------------------------------------- + +-- | Quadratic regression accumulator. +data QuadRegAcc = QuadRegAcc + { qra_n :: {-# UNPACK #-} !Int -- ^ \(n\) + , qra_w :: {-# UNPACK #-} !KBN -- ^ \(\sum w_i\) + , qra_x :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i \) + , qra_x2 :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i^2 \) + , qra_x3 :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i^3 \) + , qra_x4 :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i^4 \) + , qra_y :: {-# UNPACK #-} !KBN -- ^ \(\sum y_i \) + , qra_xy :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i y_i \) + , qra_x2y :: {-# UNPACK #-} !KBN -- ^ \(\sum x_i^2 y_i \) + , qra_y2 :: {-# UNPACK #-} !KBN -- ^ \(\sum y_i^2 \) + } + deriving Show + +instance NFData QuadRegAcc where + rnf QuadRegAcc {} = () + +-- | All-zeroes 'QuadRegAcc'. +zeroQuadRegAcc :: QuadRegAcc +zeroQuadRegAcc = QuadRegAcc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN + +-- | Add a point to quadreg accumulator. +addQuadReg + :: QuadRegAcc + -> Double -- ^ x + -> Double -- ^ y + -> QuadRegAcc +addQuadReg QuadRegAcc {..} x y = QuadRegAcc + { qra_n = qra_n + 1 + , qra_w = addKBN qra_w 1 + , qra_x = addKBN qra_x x + , qra_x2 = addKBN qra_x2 x2 + , qra_x3 = addKBN qra_x3 (x * x2) + , qra_x4 = addKBN qra_x4 (x2 * x2) + , qra_y = addKBN qra_y y + , qra_xy = addKBN qra_xy (x * y) + , qra_x2y = addKBN qra_x2y (x2 * y) + , qra_y2 = addKBN qra_y2 (y * y) + } + where + x2 = x * x + +-- | Add a weighted point to quadreg accumulator. +addQuadRegW + :: QuadRegAcc + -> Double -- ^ x + -> Double -- ^ y + -> Double -- ^ w + -> QuadRegAcc +addQuadRegW QuadRegAcc {..} x y w = QuadRegAcc + { qra_n = qra_n + 1 + , qra_w = addKBN qra_w w + , qra_x = addKBN qra_x (w * x) + , qra_x2 = addKBN qra_x2 (w * x2) + , qra_x3 = addKBN qra_x3 (w * x * x2) + , qra_x4 = addKBN qra_x4 (w * x2 * x2) + , qra_y = addKBN qra_y (w * y) + , qra_xy = addKBN qra_xy (w * x * y) + , qra_x2y = addKBN qra_x2y (w * x2 * y) + , qra_y2 = addKBN qra_y2 (w * y * y) + } + where + x2 = x * x + +-- | Convert 'QuadRegAcc' to 'LinRegAcc'. +-- +-- Using this we can try quadratic and linear fits with a single data scan. +-- +quadRegAccToLin :: QuadRegAcc -> LinRegAcc +quadRegAccToLin QuadRegAcc {..} = LinRegAcc + { lra_n = qra_n + , lra_w = qra_w + , lra_x = qra_x + , lra_x2 = qra_x2 + , lra_y = qra_y + , lra_xy = qra_xy + , lra_y2 = qra_y2 + } + +------------------------------------------------------------------------------- +-- utils +------------------------------------------------------------------------------- + +sq :: Num a => a -> a +sq x = x * x +{-# INLINE sq #-} + +iterate1 :: (b -> b) -> b -> NE.NonEmpty b +iterate1 g x = NE.cons x (iterate1 g (g x)) + +eigenSM22 :: Double -> Double -> Double -> V2 +eigenSM22 a b c = V2 ((ac + discr) / 2) ((ac - discr) / 2) + where + ac = a + c + discr = sqrt (sq (a - c) + 4 * sq b) + +-- | Levenberg-Marquard stop condition +lmStop :: Double -> Double -> Double -> Bool +lmStop lambda wssr wssr' = + lambda < 1e-20 || lambda > 1e20 || isNaN wssr' || relDiff < 1e-10 + where + relDiff = abs (wssr' - wssr) / wssr diff --git a/src/Math/Regression/Simple/LinAlg.hs b/src/Math/Regression/Simple/LinAlg.hs new file mode 100644 --- /dev/null +++ b/src/Math/Regression/Simple/LinAlg.hs @@ -0,0 +1,402 @@ +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE MultiParamTypeClasses #-} +-- | Minimil linear algebra lib. +module Math.Regression.Simple.LinAlg ( + -- * Operations + Add (..), + Eye (..), + Mult (..), + Det (..), + Inv (..), + -- * Zeros + zerosLin, + zerosQuad, + optimaQuad, + -- * Two dimensions + V2 (..), + M22 (..), + SM22 (..), + -- * Three dimensions + V3 (..), + M33 (..), + SM33 (..), +) where + +import Control.DeepSeq (NFData (..)) +import Data.Complex (Complex (..)) + +------------------------------------------------------------------------------- +-- Classes +------------------------------------------------------------------------------- + +-- | Addition +class Add a where + zero :: a + add :: a -> a -> a + +-- | Identity +class Eye a where + eye :: a + +-- | Multiplication of different things. +class Eye a => Mult a b c | a b -> c where + mult :: a -> b -> c + +-- | Determinant +class Eye a => Det a where + det :: a -> Double + +-- | Inverse +class Det a => Inv a where + inv :: a -> a + +infixl 6 `add` +infixl 7 `mult` + +instance Eye Double where + eye = 1 + +instance Add Double where + zero = 0 + add = (+) + +instance Mult Double Double Double where + mult = (*) + +instance Det Double where + det = id + +instance Inv Double where + inv = recip + +------------------------------------------------------------------------------- +-- Zeros +------------------------------------------------------------------------------- + +-- | Solve linear equation. +-- +-- >>> zerosLin (V2 1 2) +-- -2.0 +-- +zerosLin :: V2 -> Double +zerosLin (V2 a b) = negate (b / a) + +-- | Solve quadratic equation. +-- +-- >>> zerosQuad (V3 2 0 (-1)) +-- Right (-0.7071067811865476,0.7071067811865476) +-- +-- >>> zerosQuad (V3 2 0 1) +-- Left ((-0.0) :+ (-0.7071067811865476),(-0.0) :+ 0.7071067811865476) +-- +-- Double root is not treated separately: +-- +-- >>> zerosQuad (V3 1 0 0) +-- Right (-0.0,0.0) +-- +-- >>> 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)) + | otherwise = Right ((- b - sqrtDelta) / da, (-b + sqrtDelta) / da) + where + delta = b*b - 4 * a * c + sqrtDelta = sqrt delta + sqrtNDelta = sqrt (- delta) + da = 2 * a + +-- | Find an optima point. +-- +-- >>> optimaQuad (V3 1 (-2) 0) +-- 1.0 +-- +-- compare to +-- +-- >>> zerosQuad (V3 1 (-2) 0) +-- Right (0.0,2.0) +-- +optimaQuad :: V3 -> Double +optimaQuad (V3 a b _) = zerosLin (V2 (2 * a) b) + +------------------------------------------------------------------------------- +-- 2 dimensions +------------------------------------------------------------------------------- + +-- | 2d vector. Strict pair of 'Double's. +-- +-- Also used to represent linear polynomial: @V2 a b@ \(= a x + b\). +-- +data V2 = V2 !Double !Double + deriving (Eq, Show) + +instance NFData V2 where + rnf V2 {} = () + +instance Add V2 where + zero = V2 0 0 + add (V2 x y) (V2 x' y') = V2 (x + x') (y + y') + {-# INLINE zero #-} + {-# INLINE add #-} + +instance Mult Double V2 V2 where + mult k (V2 x y) = V2 (k * x) (k * y) + {-# INLINE mult #-} + +-- | 2×2 matrix. +data M22 = M22 !Double !Double !Double !Double + deriving (Eq, Show) + +instance NFData M22 where + rnf M22 {} = () + +instance Add M22 where + zero = M22 0 0 0 0 + add (M22 a b c d) (M22 a' b' c' d') = M22 (a + a') (b + b') (c + c') (d + d') + {-# INLINE zero #-} + {-# INLINE add #-} + +instance Eye M22 where + eye = M22 1 0 0 1 + {-# INLINE eye #-} + +instance Det M22 where det = det2 +instance Inv M22 where inv = inv2 + +instance Mult Double M22 M22 where + mult k (M22 a b c d) = M22 (k * a) (k * b) (k * c) (k * d) + {-# INLINE mult #-} + +instance Mult M22 V2 V2 where + mult (M22 a b c d) (V2 u v) = V2 (a * u + b * v) (c * u + d * v) + {-# INLINE mult #-} + +-- | >>> M22 1 2 3 4 `mult` eye @M22 +-- M22 1.0 2.0 3.0 4.0 +instance Mult M22 M22 M22 where + mult (M22 a b c d) (M22 x y z w) = M22 + (a * x + b * z) (a * y + b * w) + (c * x + d * z) (c * y + d * w) + {-# INLINE mult #-} + +det2 :: M22 -> Double +det2 (M22 a b c d) = a * d - b * c +{-# INLINE det2 #-} + +inv2 :: M22 -> M22 +inv2 m@(M22 a b c d) = M22 + ( d / detm) (- b / detm) + (- c / detm) ( a / detm) + where + detm = det2 m +{-# INLINE inv2 #-} + +-- | Symmetric 2x2 matrix. +data SM22 = SM22 !Double !Double !Double + deriving (Eq, Show) + +instance NFData SM22 where + rnf SM22 {} = () + +instance Add SM22 where + zero = SM22 0 0 0 + add (SM22 a b d) (SM22 a' b' d') = SM22 (a + a') (b + b') (d + d') + {-# INLINE zero #-} + {-# INLINE add #-} + +instance Eye SM22 where + eye = SM22 1 0 1 + {-# INLINE eye #-} + +instance Det SM22 where det = detS2 +instance Inv SM22 where inv = invS2 + +instance Mult Double SM22 SM22 where + mult k (SM22 a b d) = SM22 (k * a) (k * b) (k * d) + {-# INLINE mult #-} + +instance Mult SM22 V2 V2 where + mult (SM22 a b d) (V2 u v) = V2 (a * u + b * v) (b * u + d * v) + {-# INLINE mult #-} + +detS2 :: SM22 -> Double +detS2 (SM22 a b d) = a * d - b * b +{-# INLINE detS2 #-} + +invS2 :: SM22 -> SM22 +invS2 m@(SM22 a b d) = SM22 + ( d / detm) + (- b / detm) ( a / detm) + where + detm = detS2 m +{-# INLINE invS2 #-} + +------------------------------------------------------------------------------- +-- 3 dimensions +------------------------------------------------------------------------------- + +-- | 3d vector. Strict triple of 'Double's. +-- +-- Also used to represent quadratic polynomial: @V3 a b c@ \(= a x^2 + b x + c\). +data V3 = V3 !Double !Double !Double + deriving (Eq, Show) + +instance NFData V3 where + rnf V3 {} = () + +instance Add V3 where + zero = V3 0 0 0 + add (V3 x y z) (V3 x' y' z') = V3 (x + x') (y + y') (z + z') + {-# INLINE zero #-} + {-# INLINE add #-} + +instance Mult Double V3 V3 where + mult k (V3 x y z) = V3 (k * x) (k * y) (k * z) + {-# INLINE mult #-} + +-- | 3×3 matrix. +data M33 = M33 + !Double !Double !Double + !Double !Double !Double + !Double !Double !Double + deriving (Eq, Show) + +instance NFData M33 where + rnf M33 {} = () + +instance Add M33 where + zero = M33 0 0 0 0 0 0 0 0 0 + + add (M33 a b c d e f g h i) (M33 a' b' c' d' e' f' g' h' i') = M33 + (a + a') (b + b') (c + c') + (d + d') (e + e') (f + f') + (g + g') (h + h') (i + i') + {-# INLINE zero #-} + {-# INLINE add #-} + +instance Eye M33 where + eye = M33 1 0 0 + 0 1 0 + 0 0 1 + {-# INLINE eye #-} + +instance Det M33 where det = det3 +instance Inv M33 where inv = inv3 + +instance Mult Double M33 M33 where + mult k (M33 a b c d e f g h i) = M33 + (k * a) (k * b) (k * c) + (k * d) (k * e) (k * f) + (k * g) (k * h) (k * i) + {-# INLINE mult #-} + +instance Mult M33 V3 V3 where + mult (M33 a b c + d e f + g h i) (V3 u v w) = V3 + (a * u + b * v + c * w) + (d * u + e * v + f * w) + (g * u + h * v + i * w) + {-# INLINE mult #-} + +-- TODO: instance Mult M33 M33 M33 where + +det3 :: M33 -> Double +det3 (M33 a b c + d e f + g h i) + = a * (e*i-f*h) - d * (b*i-c*h) + g * (b*f-c*e) +{-# INLINE det3 #-} + +inv3 :: M33 -> M33 +inv3 m@(M33 a b c + d e f + g h i) + = M33 a' b' c' + d' e' f' + g' h' i' + where + a' = cofactor e f h i / detm + b' = cofactor c b i h / detm + c' = cofactor b c e f / detm + d' = cofactor f d i g / detm + e' = cofactor a c g i / detm + f' = cofactor c a f d / detm + g' = cofactor d e g h / detm + h' = cofactor b a h g / detm + i' = cofactor a b d e / detm + cofactor q r s t = det2 (M22 q r s t) + detm = det3 m +{-# INLINE inv3 #-} + +-- | Symmetric 3×3 matrix. +data SM33 = SM33 + !Double + !Double !Double + !Double !Double !Double + deriving (Eq, Show) + +instance NFData SM33 where + rnf SM33 {} = () + +instance Add SM33 where + zero = SM33 0 0 0 0 0 0 + + add (SM33 a d e g h i) (SM33 a' d' e' g' h' i') = SM33 + (a + a') + (d + d') (e + e') + (g + g') (h + h') (i + i') + {-# INLINE zero #-} + {-# INLINE add #-} + +instance Eye SM33 where + eye = SM33 1 + 0 1 + 0 0 1 + {-# INLINE eye #-} + +instance Det SM33 where det = detS3 +instance Inv SM33 where inv = invS3 + +instance Mult Double SM33 SM33 where + mult k (SM33 a d e g h i) = SM33 + (k * a) + (k * d) (k * e) + (k * g) (k * h) (k * i) + {-# INLINE mult #-} + +instance Mult SM33 V3 V3 where + mult (SM33 a + d e + g h i) (V3 u v w) = V3 + (a * u + d * v + g * w) + (d * u + e * v + h * w) + (g * u + h * v + i * w) + {-# INLINE mult #-} + +detS3 :: SM33 -> Double +detS3 (SM33 a + d e + g h i) + = a * (e*i-h*h) - d * (d*i-g*h) + g * (d*h-g*e) +{-# INLINE detS3 #-} + +invS3 :: SM33 -> SM33 +invS3 m@(SM33 a + d e + g h i) + = SM33 a' + d' e' + g' h' i' + where + a' = cofactor e h h i / detm + d' = cofactor h d i g / detm + e' = cofactor a g g i / detm + g' = cofactor d e g h / detm + h' = cofactor d a h g / detm + i' = cofactor a d d e / detm + cofactor q r s t = det2 (M22 q r s t) + detm = detS3 m +{-# INLINE invS3 #-} diff --git a/src/Numeric/KBN.hs b/src/Numeric/KBN.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/KBN.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FunctionalDependencies #-} +-- | https://en.wikipedia.org/wiki/KBN_summation_algorithm +-- +-- @math-functions@ has KBN-Babuška-Neumaier summation algorithm as well +-- in @Numeric.Sum@ module. +module Numeric.KBN ( + KBN (..), + zeroKBN, + getKBN, + sumKBN, + addKBN, +) where + +import Control.DeepSeq (NFData (..)) + +import qualified Data.Foldable as F + +-- $setup +-- >>> import qualified Data.List + +-- | KBN summation accumulator. +data KBN = KBN !Double !Double + deriving Show + +instance NFData KBN where + rnf KBN {} = () + +getKBN :: KBN -> Double +getKBN (KBN x e) = x + e + +zeroKBN :: KBN +zeroKBN = KBN 0 0 + +-- | KBN summation algorithm. +-- +-- >>> sumKBN (replicate 10 0.1) +-- 1.0 +-- +-- >>> Data.List.foldl' (+) 0 (replicate 10 0.1) :: Double +-- 0.9999999999999999 +-- +-- >>> sumKBN [1, 1e100, 1, -1e100] +-- 2.0 +-- +-- >>> Data.List.foldl' (+) 0 [1, 1e100, 1, -1e100] +-- 0.0 +-- +sumKBN :: F.Foldable f => f Double -> Double +sumKBN = sumKBNWith id + +-- | Generalized version of 'sumKBN'. +sumKBNWith :: F.Foldable f => (a -> Double) -> f a -> Double +sumKBNWith f xs = getKBN (F.foldl' (\k a -> addKBN k (f a)) (KBN 0 0) xs) + +-- | Add a 'Double' to 'KBN' accumulator. +addKBN :: KBN -> Double -> KBN +addKBN (KBN acc c) x = KBN acc' c' where + acc' = acc + x + c' | abs acc >= abs x = c + ((acc - acc') + x) + | otherwise = c + ((x - acc') + acc) diff --git a/test/generate-test-data.hs b/test/generate-test-data.hs new file mode 100644 --- /dev/null +++ b/test/generate-test-data.hs @@ -0,0 +1,60 @@ +module Main where + +import Control.Monad (unless) +import Data.List (unfoldr, zip5) +import System.Environment (getArgs) +import Text.Printf (printf) + +import qualified System.Random.SplitMix as SM + +main :: IO () +main = do + args <- getArgs + unless (null args) $ do + let (g1 : g2 : g3 : g4 : _) = unfoldr (Just . SM.splitSMGen) (SM.mkSMGen 42) + + writeFile "gnuplot/linear.dat" $ unlines + [ printf "%9.05f %9.05f %9.05f %9.05f" (x + dx) (y + dy) devX devY + | (x, devX, devY, px, py) <- zip5 + [ 1 .. 20 :: Double ] + devsX + devsY + (normals (doubles g1)) + (normals (doubles g2)) + , let y = 3 * x + 5 + , let dx = devX * px + , let dy = devY * py + ] + + writeFile "gnuplot/quad.dat" $ unlines + [ printf "%9.05f %9.05f %9.05f %9.05f" (x + dx) (y + dy) devX devY + | (x, devX, devY, px, py) <- zip5 + [ 1 .. 20 :: Double ] + devsX + devsY + (normals (doubles g3)) + (normals (doubles g4)) + , let y = 0.1 * x * x - 3 * x + 5 + , let dx = devX * px + , let dy = devY * py + ] + +boxMuller :: Double -> Double -> (Double, Double) +boxMuller u v = (x, y) where + x = k * sin (2 * pi * v) + y = k * cos (2 * pi * u) + k = sqrt (negate 2 * log u) + +doubles :: SM.SMGen -> [Double] +doubles g = let (d, g') = SM.nextDouble g in d : doubles g' + +normals :: [Double] -> [Double] +normals [] = [] +normals [_] = [] +normals (u:v:ds) = let (x, y) = boxMuller u v in x : y : normals ds + +devsY :: [Double] +devsY = 1.2 : 1.4 : 1.6 : devsY + +devsX :: [Double] +devsX = 0.2 : 0.3 : devsX diff --git a/test/regression-simple-tests.hs b/test/regression-simple-tests.hs new file mode 100644 --- /dev/null +++ b/test/regression-simple-tests.hs @@ -0,0 +1,333 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveFoldable #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE DeriveTraversable #-} +module Main (main) where + +import Data.List (zip4) +import Test.Tasty (TestTree, defaultMain, testGroup, withResource) +import Test.Tasty.HUnit (assertEqual, testCase) + +import qualified Data.Foldable as F +import qualified Data.List.NonEmpty as NE +import qualified Data.Traversable as T + +#if __GLASGOW_HASKELL__ >= 704 +import qualified Numeric.AD.Mode.Reverse.Double as AD +import qualified Statistics.Distribution as S +import qualified Statistics.Distribution.ChiSquared as S +#endif + +import Math.Regression.Simple +import Numeric.KBN (sumKBN) + +------------------------------------------------------------------------------- +-- Main +------------------------------------------------------------------------------- + +main :: IO () +main = defaultMain $ testGroup "regression-simple" + [ linearTests + , quadraticTests + , lm1Tests + , lm2Tests + ] + +------------------------------------------------------------------------------- +-- data +------------------------------------------------------------------------------- + +withData :: FilePath -> (IO [(Double, Double, Double, Double)] -> TestTree) -> TestTree +withData fp = withResource acquire release where + acquire = do + contents <- readFile ("gnuplot/" ++ fp) + return $ map (quad . map read . words) $ lines contents + + release _ = return () + + quad :: [Double] -> (Double, Double, Double, Double) + quad (x:y:dx:dy:_) = (x,y,dx,dy) + quad _ = error "invalid data" + +------------------------------------------------------------------------------- +-- Linear +------------------------------------------------------------------------------- + +linearTests :: TestTree +linearTests = withData "linear.dat" $ \load -> testGroup "linear" + [ testCase "no-errors" $ do + linearData <- load + let fit = linearFit (\(x,y,_,_) -> (x,y)) linearData + assertEqual "params" (V2 2.95689 6.04617) (round' (fitParams fit)) + assertEqual "errors" (V2 7.9788e-2 0.95195) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 75.6356 (round' (fitWSSR fit)) + + , testCase "y-errors" $ do + linearData <- load + let fit = linearWithYerrors (\(x,y,_,dy) -> (x,y,dy)) linearData + assertEqual "params" (V2 2.97271 5.91878) (round' (fitParams fit)) + assertEqual "errors" (V2 7.722e-2 0.91882) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 38.8345 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 2.999e-3 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + + , testCase "xy-errors" $ do + linearData <- load + let fit = nth 5 $ linearWithXYerrors (\(x,y,dx,dy) -> (x,y,dx,dy)) linearData + assertEqual "params" (V2 2.97021 5.99061) (round' (fitParams fit)) + assertEqual "errors" (V2 7.6542e-2 0.90917) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 29.141 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 4.6683e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + + , testCase "yx-errors" $ do + linearData <- load + let fit = nth 5 $ linearWithXYerrors (\(x,y,dx,dy) -> (y,x,dy,dx)) linearData + assertEqual "params" (V2 0.33271 (-1.87107)) (round' (fitParams fit)) + assertEqual "errors" (V2 8.5724e-3 0.34855) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 29.3171 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 4.4639e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + ] + +nth :: Int -> NE.NonEmpty a -> a +nth n (x NE.:| xs) = go n x xs where + go _ z [] = z + go m z (y:ys) = if m <= 0 then z else go (m - 1) y ys + +------------------------------------------------------------------------------- +-- Quad +------------------------------------------------------------------------------- + +quadraticTests :: TestTree +quadraticTests = withData "quad.dat" $ \load -> testGroup "quad" + [ testCase "no-errors" $ do + quadraticData <- load + let fit = quadraticFit (\(x,y,_,_) -> (x,y)) quadraticData + assertEqual "params" (V3 0.11487 (-3.34246) 6.63601) (round' (fitParams fit)) + assertEqual "errors" (V3 1.0297e-2 0.22674 1.07032) (round' (fitErrors fit)) + assertEqual "ndf" 17 (round' (fitNDF fit)) + assertEqual "wssr" 33.9104 (round' (fitWSSR fit)) + + , testCase "y-errors" $ do + quadraticData <- load + let fit = quadraticWithYerrors (\(x,y,_,dy) -> (x,y,dy)) quadraticData + assertEqual "params" (V3 0.11156 (-3.27481) 6.25286) (round' (fitParams fit)) + assertEqual "errors" (V3 9.7603e-3 0.21331 0.99362) (round' (fitErrors fit)) + assertEqual "ndf" 17 (round' (fitNDF fit)) + assertEqual "wssr" 16.793 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + let q = S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit) + assertEqual "P" 0.46847 (round' (1 - q)) +#endif + + , testCase "xy-errors" $ do + quadraticData <- load + let fit = nth 5 $ quadraticWithXYerrors (\(x,y,dx,dy) -> (x,y,dx,dy)) quadraticData + assertEqual "params" (V3 0.11222 (-3.29575) 6.39876) (round' (fitParams fit)) + assertEqual "errors" (V3 9.9372e-3 0.22027 1.06196) (round' (fitErrors fit)) + assertEqual "ndf" 17 (round' (fitNDF fit)) + assertEqual "wssr" 15.6318 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + let q = S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit) + assertEqual "P" 0.55007 (round' (1 - q)) +#endif + ] + +------------------------------------------------------------------------------- +-- LM +------------------------------------------------------------------------------- + +lm1Tests :: TestTree +lm1Tests = withData "linear.dat" $ \load -> testGroup "lm1" + [ testCase "no-errors" $ do + linearData <- load + + let scale a (x, y, _, _) = case scaleGrad' (H2 a x) of + (f, H2 da _) -> (y, f, da) + + let fit = NE.last $ levenbergMarquardt1 scale 1 linearData + assertEqual "params" 3.03374 (round' (fitParams fit)) + assertEqual "errors" 3.8628e-2 (round' (fitErrors fit)) + assertEqual "ndf" 19 (round' (fitNDF fit)) + assertEqual "wssr" 80.7105 (round' (fitWSSR fit)) + + , testCase "y-errors" $ do + linearData <- load + let scaleY a (x, y, _, dy) = case scaleGrad' (H2 a x) of + (f, H2 da _) -> (y, f, da, dy) + let fit = NE.last $ levenbergMarquardt1WithYerrors scaleY 1 linearData + assertEqual "params" 3.04015 (round' (fitParams fit)) + assertEqual "errors" 3.7604e-2 (round' (fitErrors fit)) + assertEqual "ndf" 19 (round' (fitNDF fit)) + assertEqual "wssr" 40.9918 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 2.4195e-3 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + + , testCase "xy-errors" $ do + linearData <- load + let scaleXY a (x, y, dx, dy) = case scaleGrad' (H2 a x) of + (f, H2 da f') -> (y, f, da, f', dx, dy) + let fit = NE.last $ levenbergMarquardt1WithXYerrors scaleXY 1 linearData + assertEqual "params" 3.04315 (round' (fitParams fit)) + assertEqual "errors" 3.7477e-2 (round' (fitErrors fit)) + assertEqual "ndf" 19 (round' (fitNDF fit)) + assertEqual "wssr" 30.7021 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 4.3516e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + ] + +lm2Tests :: TestTree +lm2Tests = withData "linear.dat" $ \load -> testGroup "lm2" + [ testCase "no-errors" $ do + linearData <- load + let lin (V2 a b) (x, y, _, _) = case linearGrad' (H3 a b x) of + (f, H3 da db _) -> (y, f, V2 da db) + + let fit = NE.last $ levenbergMarquardt2 lin (V2 1 1) linearData + assertEqual "params" (V2 2.95689 6.04617) (round' (fitParams fit)) + assertEqual "errors" (V2 7.9788e-2 0.95195) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 75.6356 (round' (fitWSSR fit)) + + , testCase "y-errors" $ do + linearData <- load + let linY (V2 a b) (x, y, _, dy) = case linearGrad' (H3 a b x) of + (f, H3 da db _) -> (y, f, V2 da db, dy) + let fit = NE.last $ levenbergMarquardt2WithYerrors linY (V2 1 1) linearData + assertEqual "params" (V2 2.97271 5.91882) (round' (fitParams fit)) + assertEqual "errors" (V2 7.722e-2 0.91882) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 38.8345 (round' (fitWSSR fit)) + + , testCase "xy-errors" $ do + linearData <- load + let linXY (V2 a b) (x, y, dx, dy) = case linearGrad' (H3 a b x) of + (f, H3 da db f') -> (y, f, V2 da db, f', dx, dy) + let fit = NE.last $ levenbergMarquardt2WithXYerrors linXY (V2 1 1) linearData + assertEqual "params" (V2 2.97021 5.99061) (round' (fitParams fit)) + assertEqual "errors" (V2 7.6542e-2 0.90917) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 29.141 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 4.6683e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + + , testCase "yx-errors" $ do + -- with x and y flipped: + linearData <- load + let linYX (V2 a b) (y, x, dy, dx) = case linearGrad' (H3 a b x) of + (f, H3 da db f') -> (y, f, V2 da db, f', dx, dy) + let fit = NE.last $ levenbergMarquardt2WithXYerrors linYX (V2 1 1) linearData + assertEqual "params" (V2 0.33402 (-1.92156)) (round' (fitParams fit)) + assertEqual "errors" (V2 8.5785e-3 0.3488) (round' (fitErrors fit)) + assertEqual "ndf" 18 (round' (fitNDF fit)) + assertEqual "wssr" 29.1822 (round' (fitWSSR fit)) + +#if __GLASGOW_HASKELL__ >= 704 + assertEqual "P" 4.6197e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + +#if __GLASGOW_HASKELL__ >= 704 + , testCase "orear-example" $ do + let orearData :: [(Double, Double, Double, Double)] + orearData = zip4 + [22000, 22930,23880,25130,26390] + [-4.017,-2.742,-1.1478,1.491,6.873] + [440,470,500,530,540] + [0.50,0.25,0.08,0.09,1.90] + + let orearXY (V2 a b) (x, y, dx, dy) = case AD.grad' orearF (H3 a b x) of + (f, H3 da db f') -> (y, f, V2 da db, f', dx, dy) + + let wssr0 = sumKBN + [ sq (y - f) * w + | d <- orearData + , let a1 = 1e-3 + , let a2 = 6e5 + , let (y, f, _, f', dx, dy) = orearXY (V2 a1 a2) d + , let w = recip $ sq (f' * dx) + sq dy + ] + + assertEqual "wssr0" 3.82243 (round' wssr0) + + let fit = NE.last $ levenbergMarquardt2WithXYerrors orearXY (V2 1e-3 6e5) orearData + + assertEqual "params" (V2 1.0163e-3 593725.0) (round' (fitParams fit)) + assertEqual "errors" (V2 1.7025e-4 95284.8) (round' (fitErrors fit)) + assertEqual "ndf" 3 (round' (fitNDF fit)) + assertEqual "wssr" 2.18668 (round' (fitWSSR fit)) + + assertEqual "P" 0.53458 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) +#endif + + ] + +data H2 a = H2 a a deriving (Functor, F.Foldable, T.Traversable) +data H3 a = H3 a a a deriving (Functor, F.Foldable, T.Traversable) + +scaleF :: Num a => H2 a -> a +scaleF (H2 a x) = a * x + 5 + +scaleGrad' :: H2 Double -> (Double, H2 Double) +#if __GLASGOW_HASKELL__ >= 704 +scaleGrad' = AD.grad' scaleF +#else +scaleGrad' (H2 a x) = (a * x + 5, H2 x a) +#endif + +linearF :: Num a => H3 a -> a +linearF (H3 a b x) = a * x + b + +linearGrad' :: H3 Double -> (Double, H3 Double) +#if __GLASGOW_HASKELL__ >= 704 +linearGrad' = AD.grad' linearF +#else +linearGrad' (H3 a b x) = (a * x + b, H3 x 1 a) +#endif + +orearF :: Fractional a => H3 a -> a +orearF (H3 a b x) = a * x - b / x + +sq :: Num a => a -> a +sq x = x * x + +------------------------------------------------------------------------------- +-- Round +------------------------------------------------------------------------------- + +class Round a where + round' :: a -> a + +instance Round Double where + round' 0 = 0 + round' x = fromInteger (round (x * rat)) / rat + where + mag = truncate (logBase 10 (abs x)) :: Int + rat = 10 ^ (5 - mag) + +instance Round Int where + round' = id + +instance Round V2 where + round' (V2 x y) = V2 (round' x) (round' y) + +instance Round V3 where + round' (V3 x y z) = V3 (round' x) (round' y) (round' z)