From 2e48ffd1a395817288b8271299eebd0e483407af Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 6 Apr 2010 18:24:16 +0000 Subject: some changes in GSL.Fitting --- examples/fitting.hs | 6 ++-- lib/Numeric/GSL/Fitting.hs | 72 ++++++++++++++++++++++++++------------ lib/Numeric/GSL/gsl-aux.c | 1 + lib/Numeric/LinearAlgebra/Tests.hs | 12 ++++--- 4 files changed, 62 insertions(+), 29 deletions(-) diff --git a/examples/fitting.hs b/examples/fitting.hs index 8298c52..a8f6b1c 100644 --- a/examples/fitting.hs +++ b/examples/fitting.hs @@ -8,15 +8,15 @@ sigma = 0.1 ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs) + scalar sigma * (randomVector 0 Gaussian 40) -dat :: [([Double],[Double],Double)] +dat :: [([Double],([Double],Double))] -dat = zipWith3 (,,) xs ys (repeat sigma) +dat = zip xs (zip ys (repeat sigma)) expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] -(sol,path) = fitModel 1E-4 1E-4 20 (resM expModel, resD expModelDer) dat [1,0,0] +(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] main = do print dat diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs index 6e91b2d..b1f3a32 100644 --- a/lib/Numeric/GSL/Fitting.hs +++ b/lib/Numeric/GSL/Fitting.hs @@ -14,17 +14,17 @@ Nonlinear Least-Squares Fitting The example program in the GSL manual (see examples/fitting.hs): @dat = [ - ([0.0],[6.0133918608118675],0.1), - ([1.0],[5.5153769909966535],0.1), - ([2.0],[5.261094606015287],0.1), + ([0.0],([6.0133918608118675],0.1)), + ([1.0],([5.5153769909966535],0.1)), + ([2.0],([5.261094606015287],0.1)), ... - ([39.0],[1.0619821710802808],0.1)] + ([39.0],([1.0619821710802808],0.1))] expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] -(sol,path) = fitModel 1E-4 1E-4 20 (resM expModel, resD expModelDer) dat [1,0,0] +(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] \> path (6><5) @@ -46,7 +46,7 @@ module Numeric.GSL.Fitting ( -- * Levenberg-Marquardt nlFitting, FittingMethod(..), -- * Utilities - fitModel, resM, resD + fitModelScaled, fitModel ) where import Data.Packed.Internal @@ -57,7 +57,8 @@ import Numeric.GSL.Internal ------------------------------------------------------------------------- -data FittingMethod = LevenbergMarquardt -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom. +data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom. + | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately. deriving (Enum,Eq,Show,Bounded) @@ -104,32 +105,52 @@ checkdim2 n p m ------------------------------------------------------------ -err model dat vsol = zip sol errs where +err (model,deriv) dat vsol = zip sol errs where sol = toList vsol c = max 1 (chi/sqrt (fromIntegral dof)) dof = length dat - (rows cov) - chi = pnorm PNorm2 (fromList $ cost (fst model) dat sol) - js = fromLists $ jacobian (snd model) dat sol + chi = pnorm PNorm2 (fromList $ cost (resMs model) dat sol) + js = fromLists $ jacobian (resDs deriv) dat sol cov = inv $ trans js <> js errs = toList $ scalar c * sqrt (takeDiag cov) +-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and +-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of +-- instances (x, (y,sigma)) to be fitted. + +fitModelScaled + :: Double -- ^ absolute tolerance + -> Double -- ^ relative tolerance + -> Int -- ^ maximum number of iterations allowed + -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives) + -> [(x, ([Double], Double))] -- ^ instances + -> [Double] -- ^ starting point + -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path +fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where + (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit + (fromList . cost (resMs model) dt . toList) + (fromLists . jacobian (resDs deriv) dt . toList) + (fromList xin) + + + -- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and --- Jacobian are automatically built from a model f vs x = 0 and its derivatives, and a list of --- instances x to be fitted. +-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of +-- instances (x,y) to be fitted. fitModel :: Double -- ^ absolute tolerance -> Double -- ^ relative tolerance -> Int -- ^ maximum number of iterations allowed -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives) - -> [x] -- ^ instances + -> [(x, [Double])] -- ^ instances -> [Double] -- ^ starting point - -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path -fitModel epsabs epsrel maxit model dt xin = (err model dt sol, path) where + -> ([Double], Matrix Double) -- ^ solution and optimization path +fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit - (fromList . cost (fst model) dt . toList) - (fromLists . jacobian (snd model) dt . toList) + (fromList . cost (resM model) dt . toList) + (fromLists . jacobian (resD deriv) dt . toList) (fromList xin) cost model ds vs = concatMap (model vs) ds @@ -137,11 +158,18 @@ cost model ds vs = concatMap (model vs) ds jacobian modelDer ds vs = concatMap (modelDer vs) ds -- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'. -resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double], Double) -> [Double] -resM m v = \(x,ys,s) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s +resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double] +resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s --- | Associated derivative for 'resM'. -resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double], Double) -> [[Double]] -resD m v = \(x,_,s) -> map (map (/s)) (m v x) +-- | Associated derivative for 'resMs'. +resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]] +resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x) +-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent +-- to 'resMs' with all sigmas = 1. +resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double] +resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b +-- | Associated derivative for 'resM'. +resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]] +resD m v = \(x,_) -> m v x diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index f258fb1..f5a95f4 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -812,6 +812,7 @@ int nlfit(int method, int f(int, double*, int, double*), switch(method) { case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; } + case 1 : { T = gsl_multifit_fdfsolver_lmder; break; } default: ERROR(BAD_CODE); } diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index c9408a6..60c60e6 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs @@ -143,19 +143,23 @@ odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.36434890871 --------------------------------------------------------------------- -fittingTest = utest "levmar" ok +fittingTest = utest "levmar" (ok1 && ok2) where xs = map return [0 .. 39] sigma = 0.1 ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs) + scalar sigma * (randomVector 0 Gaussian 40) - dat = zipWith3 (,,) xs ys (repeat sigma) + dats = zip xs (zip ys (repeat sigma)) + dat = zip xs ys expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] expModelDer [a,lambda,_b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] - sol = fst $ fitModel 1E-4 1E-4 20 (resM expModel, resD expModelDer) dat [1,0,0] - ok = and (zipWith f sol [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d + sols = fst $ fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dats [1,0,0] + sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] + + ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d + ok2 = pnorm PNorm2 (fromList (map fst sols) - fromList sol) < 1E-5 --------------------------------------------------------------------- -- cgit v1.2.3