regression-simple 0.1.1 → 0.2

raw patch · 9 files changed

+2166/−579 lines, 9 filesdep +addep +deepseqdep +math-functionsdep −vectordep ~base

Dependencies added: ad, deepseq, math-functions, regression-simple, semigroups, splitmix, statistics, tasty, tasty-hunit

Dependencies removed: vector

Dependency ranges changed: base

Files

• changelog.md +8/−1
• gnuplot/linear.dat +20/−0
• gnuplot/quad.dat +20/−0
• regression-simple.cabal +58/−8
• src/Math/Regression/Simple.hs +1203/−570
• src/Math/Regression/Simple/LinAlg.hs +402/−0
• src/Numeric/KBN.hs +62/−0
• test/generate-test-data.hs +60/−0
• test/regression-simple-tests.hs +333/−0

@@ -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.-

@@ -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

@@ -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

@@ -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 <oleg.grenrus@iki.fi> 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

@@ -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

@@ -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 #-}

@@ -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)  

@@ -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

@@ -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)