diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,8 @@ +# 0.2.1 + +* Better lambda0 guesses, hopefully +* Support GHC-9.6 + # 0.2 Large refactoring: 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.2 +version: 0.2.1 synopsis: Simple linear and quadratic regression category: Math description: @@ -38,8 +38,9 @@ || ==8.8.4 || ==8.10.7 || ==9.0.2 - || ==9.2.5 - || ==9.4.4 + || ==9.2.7 + || ==9.4.5 + || ==9.6.1 source-repository head type: git @@ -55,7 +56,7 @@ Numeric.KBN build-depends: - , base >=4.3 && <4.18 + , base >=4.3 && <4.19 , deepseq if !impl(ghc >=8.0) diff --git a/src/Math/Regression/Simple.hs b/src/Math/Regression/Simple.hs --- a/src/Math/Regression/Simple.hs +++ b/src/Math/Regression/Simple.hs @@ -405,15 +405,16 @@ -- 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 -> a -> (Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x_i)\) -> 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 + lambda0 = sqrt (c11 / fromIntegral n) where + n = lm1_n acc0 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 @@ -453,15 +454,16 @@ -- | '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 -> a -> (Double, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x_i), 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 + lambda0 = sqrt (c11 / fromIntegral n) where + n = lm1_n acc0 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 @@ -501,7 +503,7 @@ -- | '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 -> a -> (Double, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x_i), \delta y_i\) -> Double -- ^ initial parameter, \(\beta_0\) -> f a -- ^ data, \(d\) -> NE.NonEmpty (Fit Double) -- ^ non-empty list of iteration results @@ -511,7 +513,7 @@ -- | '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 -> a -> (Double, Double, Double, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x_i), \partial_x f(\beta, x_i), \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 @@ -563,12 +565,10 @@ -- -- >>> 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 1.2782 1.4831) (V2 0.57784 1.4416) 3 9.5041 +-- Fit (V2 1.7254 1.4730) (V2 0.18820 0.46952) 3 1.0082 +-- Fit (V2 1.9796 0.95226) (V2 0.09683 0.24157) 3 0.26687 +-- Fit (V2 2.0060 0.88759) (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: @@ -595,12 +595,10 @@ -- >>> 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 1.2782 1.4831) (V2 0.57784 1.4416) 3 9.5041 +-- Fit (V2 1.7254 1.4730) (V2 0.18820 0.46952) 3 1.0082 +-- Fit (V2 1.9796 0.95226) (V2 0.09683 0.24157) 3 0.26687 +-- Fit (V2 2.0060 0.88759) (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 @@ -613,13 +611,12 @@ -- >>> 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 +-- Fit (V2 0.61786 0.36360) (V2 0.23270 0.50259) 5 0.26730 +-- Fit (V2 0.39270 0.49787) (V2 0.05789 0.24170) 5 0.01237 +-- Fit (V2 0.36121 0.54525) (V2 0.04835 0.23315) 5 0.00785 +-- Fit (V2 0.36168 0.55530) (V2 0.04880 0.23790) 5 0.00784 +-- Fit (V2 0.36182 0.55620) (V2 0.04885 0.23826) 5 0.00784 +-- Fit (V2 0.36184 0.55626) (V2 0.04885 0.23829) 5 0.00784 -- -- We get the same result as in the article: 0.362 and 0.556 -- @@ -628,7 +625,7 @@ -- 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 -> a -> (Double, Double, V2)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i)\) -> V2 -- ^ initial parameters, \(\beta_0\) -> f a -- ^ data, \(d\) -> NE.NonEmpty (Fit V2) -- ^ non-empty list of iteration results @@ -637,11 +634,10 @@ 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 + lambda0 = sqrt $ (c11 + c22) / fromIntegral n / 2 where - V2 l1 l2 = eigenSM22 c11 c12 c22 + n = lm2_n acc0 c11 = getKBN $ lm2_c11 acc0 - c12 = getKBN $ lm2_c11 acc0 c22 = getKBN $ lm2_c22 acc0 loop lambda beta acc @@ -688,18 +684,17 @@ -- 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 -> a -> (Double, Double, V2, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i), 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 + lambda0 = sqrt $ (c11 + c22) / fromIntegral n / 2 where - V2 l1 l2 = eigenSM22 c11 c12 c22 + n = lm2_n acc0 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 @@ -743,7 +738,7 @@ -- | '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 -> a -> (Double, Double, V2, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i), \delta y_i\) -> V2 -- ^ initial parameters, \(\beta_0\) -> f a -- ^ data, \(d\) -> NE.NonEmpty (Fit V2) -- ^ non-empty list of iteration results @@ -753,7 +748,7 @@ -- | '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 -> a -> (Double, Double, V2, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i), \partial_x f(\beta, x_i), \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 @@ -810,7 +805,7 @@ -- -- >>> 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 +-- 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: -- @@ -819,7 +814,7 @@ -- 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 -> a -> (Double, Double, V3)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i)\) -> V3 -- ^ initial parameters, \(\beta_0\) -> f a -- ^ data, \(d\) -> NE.NonEmpty (Fit V3) -- ^ non-empty list of iteration results @@ -828,17 +823,11 @@ 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] + lambda0 = sqrt $ (c11 + c22 + c33) / fromIntegral n / 3 where + n = lm3_n acc0 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 @@ -891,22 +880,18 @@ -- | '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 -> a -> (Double, Double, V3, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i), 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] + lambda0 = sqrt $ (c11 + c22 + c33) / fromIntegral n / 3 where + n = lm3_n acc0 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 @@ -957,7 +942,7 @@ -- | '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 -> a -> (Double, Double, V3, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i), \delta y_i\) -> V3 -- ^ initial parameters, \(\beta_0\) -> f a -- ^ data, \(d\) -> NE.NonEmpty (Fit V3) -- ^ non-empty list of iteration results @@ -967,7 +952,7 @@ -- | '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 -> a -> (Double, Double, V3, Double, Double, Double)) -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x_i), \partial_x f(\beta, x_i), \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 @@ -1192,12 +1177,6 @@ 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 diff --git a/test/regression-simple-tests.hs b/test/regression-simple-tests.hs --- a/test/regression-simple-tests.hs +++ b/test/regression-simple-tests.hs @@ -159,7 +159,7 @@ (f, H2 da _) -> (y, f, da) let fit = NE.last $ levenbergMarquardt1 scale 1 linearData - assertEqual "params" 3.03374 (round' (fitParams fit)) + assertEqual "params" 3.03373 (round' (fitParams fit)) assertEqual "errors" 3.8628e-2 (round' (fitErrors fit)) assertEqual "ndf" 19 (round' (fitNDF fit)) assertEqual "wssr" 80.7105 (round' (fitWSSR fit)) @@ -191,6 +191,16 @@ #if __GLASGOW_HASKELL__ >= 704 assertEqual "P" 4.3516e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) #endif + + , testCase "issue-8" $ do + let dat = [(1e5,1e6),(1e6,1e7)] :: [(Double, Double)] + let func a (x, y) = (y, a * x * log x, x * log x) + + let fit = NE.last $ levenbergMarquardt1 func 1.0 dat + assertEqual "params" 0.72482 (round' (fitParams fit)) + assertEqual "errors" 1.1981e-2 (round' (fitErrors fit)) + assertEqual "ndf" 1 (round' (fitNDF fit)) + assertEqual "wssr" 2.75862e10 (round' (fitWSSR fit)) ] lm2Tests :: TestTree @@ -211,7 +221,7 @@ 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 "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)) @@ -236,13 +246,13 @@ 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 "params" (V2 0.3334 (-1.8971)) (round' (fitParams fit)) + assertEqual "errors" (V2 8.5742e-3 0.34862) (round' (fitErrors fit)) assertEqual "ndf" 18 (round' (fitNDF fit)) - assertEqual "wssr" 29.1822 (round' (fitWSSR fit)) + assertEqual "wssr" 29.2361 (round' (fitWSSR fit)) #if __GLASGOW_HASKELL__ >= 704 - assertEqual "P" 4.6197e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) + assertEqual "P" 4.5568e-2 (round' (1 - S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit))) #endif #if __GLASGOW_HASKELL__ >= 704 @@ -318,10 +328,11 @@ instance Round Double where round' 0 = 0 - round' x = fromInteger (round (x * rat)) / rat + round' x + | mag > 5 = let rat = 10 ^ (mag - 5) in fromInteger (round (x / rat)) * rat + | otherwise = let rat = 10 ^ (5 - mag) in fromInteger (round (x * rat)) / rat where mag = truncate (logBase 10 (abs x)) :: Int - rat = 10 ^ (5 - mag) instance Round Int where round' = id