summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL/Fitting.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL/Fitting.hs')
-rw-r--r--lib/Numeric/GSL/Fitting.hs179
1 files changed, 0 insertions, 179 deletions
diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs
deleted file mode 100644
index c4f3a91..0000000
--- a/lib/Numeric/GSL/Fitting.hs
+++ /dev/null
@@ -1,179 +0,0 @@
1{- |
2Module : Numeric.GSL.Fitting
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5
6Maintainer : Alberto Ruiz (aruiz at um dot es)
7Stability : provisional
8Portability : uses ffi
9
10Nonlinear Least-Squares Fitting
11
12<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>
13
14The example program in the GSL manual (see examples/fitting.hs):
15
16@
17dat = [
18 ([0.0],([6.0133918608118675],0.1)),
19 ([1.0],([5.5153769909966535],0.1)),
20 ([2.0],([5.261094606015287],0.1)),
21 ...
22 ([39.0],([1.0619821710802808],0.1))]
23
24expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
25
26expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
27
28(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
29@
30
31>>> path
32(6><5)
33 [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
34 , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662
35 , 3.0, 9.5807893736187, 4.948995119561291, 0.11942927999921617, 1.0945766509238248
36 , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608
37 , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375
38 , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ]
39>>> sol
40[(5.045357818922331,6.027976702418132e-2),
41(0.10404905846029407,3.157045047172834e-3),
42(1.0192487112786812,3.782067731353722e-2)]
43
44-}
45-----------------------------------------------------------------------------
46
47module Numeric.GSL.Fitting (
48 -- * Levenberg-Marquardt
49 nlFitting, FittingMethod(..),
50 -- * Utilities
51 fitModelScaled, fitModel
52) where
53
54import Data.Packed.Internal
55import Numeric.LinearAlgebra
56import Numeric.GSL.Internal
57
58import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
59import Foreign.C.Types
60import System.IO.Unsafe(unsafePerformIO)
61
62-------------------------------------------------------------------------
63
64data 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.
65 | 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.
66 deriving (Enum,Eq,Show,Bounded)
67
68
69-- | Nonlinear multidimensional least-squares fitting.
70nlFitting :: FittingMethod
71 -> Double -- ^ absolute tolerance
72 -> Double -- ^ relative tolerance
73 -> Int -- ^ maximum number of iterations allowed
74 -> (Vector Double -> Vector Double) -- ^ function to be minimized
75 -> (Vector Double -> Matrix Double) -- ^ Jacobian
76 -> Vector Double -- ^ starting point
77 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
78
79nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
80
81nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
82 let p = dim xiv
83 n = dim (f xiv)
84 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
85 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
86 rawpath <- createMatrix RowMajor maxit (2+p)
87 app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
88 let it = round (rawpath @@> (maxit-1,0))
89 path = takeRows it rawpath
90 [sol] = toRows $ dropRows (it-1) path
91 freeHaskellFunPtr fp
92 freeHaskellFunPtr jp
93 return (subVector 2 p sol, path)
94
95foreign import ccall safe "nlfit"
96 c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
97
98-------------------------------------------------------
99
100checkdim1 n _p v
101 | dim v == n = v
102 | otherwise = error $ "Error: "++ show n
103 ++ " components expected in the result of the function supplied to nlFitting"
104
105checkdim2 n p m
106 | rows m == n && cols m == p = m
107 | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
108 ++ " Jacobian expected in nlFitting"
109
110------------------------------------------------------------
111
112err (model,deriv) dat vsol = zip sol errs where
113 sol = toList vsol
114 c = max 1 (chi/sqrt (fromIntegral dof))
115 dof = length dat - (rows cov)
116 chi = norm2 (fromList $ cost (resMs model) dat sol)
117 js = fromLists $ jacobian (resDs deriv) dat sol
118 cov = inv $ trans js <> js
119 errs = toList $ scalar c * sqrt (takeDiag cov)
120
121
122
123-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
124-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
125-- instances (x, (y,sigma)) to be fitted.
126
127fitModelScaled
128 :: Double -- ^ absolute tolerance
129 -> Double -- ^ relative tolerance
130 -> Int -- ^ maximum number of iterations allowed
131 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
132 -> [(x, ([Double], Double))] -- ^ instances
133 -> [Double] -- ^ starting point
134 -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
135fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
136 (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
137 (fromList . cost (resMs model) dt . toList)
138 (fromLists . jacobian (resDs deriv) dt . toList)
139 (fromList xin)
140
141
142
143-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
144-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
145-- instances (x,y) to be fitted.
146
147fitModel :: Double -- ^ absolute tolerance
148 -> Double -- ^ relative tolerance
149 -> Int -- ^ maximum number of iterations allowed
150 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
151 -> [(x, [Double])] -- ^ instances
152 -> [Double] -- ^ starting point
153 -> ([Double], Matrix Double) -- ^ solution and optimization path
154fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
155 (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
156 (fromList . cost (resM model) dt . toList)
157 (fromLists . jacobian (resD deriv) dt . toList)
158 (fromList xin)
159
160cost model ds vs = concatMap (model vs) ds
161
162jacobian modelDer ds vs = concatMap (modelDer vs) ds
163
164-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
165resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
166resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
167
168-- | Associated derivative for 'resMs'.
169resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
170resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
171
172-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
173-- to 'resMs' with all sigmas = 1.
174resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
175resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
176
177-- | Associated derivative for 'resM'.
178resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
179resD m v = \(x,_) -> m v x
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-15 09:19:16 +0000