• incanter

# non-linear-model

## incanter.optimize

• (non-linear-model f y x start & options)

Determine the nonlinear least-squares estimates of the
parameters of a nonlinear model.
Based on R's nls (non-linear least squares) function.

Arguments:
f -- model function, takes two argumetns the first a list of parameters
that are to be estimated, and an x value.
y -- sequence of dependent data
x -- sequence of independent data
start -- start values for the parameters to be estimated

Options:
:method (default :gauss-newton) other option :newton-raphson
:tol (default 1E-5)
:max-iter (default 200)

Returns: a hash-map containing the following fields:
:method -- the method used
:coefs -- the parameter estimates
:hessian -- the estimated hessian, if available
:iterations -- the number of iterations performed
:fitted -- the fitted values of y (i.e. y-hat)
:x -- the independent data values
:y -- the dependent data values

Examples:

;; example 1
(use '(incanter core optimize datasets charts))
;; define the Michaelis-Menton model function
;; y = a + (b - a)*x/(c + x)
(defn f [theta x]
(let [[a b c] theta]
(plus a (div (mult x (minus b a)) (plus c x)))))

(def start [20 200 100])
(def data (to-matrix (get-dataset :thurstone)))
(def x (sel data :cols 1))
(def y (sel data :cols 0))
;; view the data
(def plot (scatter-plot x y))
(view plot)

(def nlm (non-linear-model f y x start))

;; example 2
(use '(incanter core optimize datasets charts))
;; Chwirut data set from NIST
(def data (to-matrix (get-dataset :chwirut)))
(def x (sel data :cols 1))
(def y (sel data :cols 0))

;; define model function: y = exp(-b1*x)/(b2+b3*x) + e
(defn f [theta x]
(let [[b1 b2 b3] theta]
(div (exp (mult (minus b1) x)) (plus b2 (mult b3 x)))))

(def plot (scatter-plot x y :legend true))
(view plot)

;; the newton-raphson algorithm fails to converge to the correct solution
;; using first set of start values from NIST, but the default gauss-newton
;; algorith converges to the correct solution.

(def start1 [0.1 0.01 0.02])
(add-lines plot x (f start1 x))
(def nlm1 (non-linear-model f y x start1))

;; both algorithms converges with second set of start values from NIST
(def start2 [0.15 0.008 0.010])
(add-lines plot x (f start2 x))
(def nlm2 (non-linear-model f y x start2))

### Source incanter/optimize.clj:590 top

(defn non-linear-model
"
Determine the nonlinear least-squares estimates of the
parameters of a nonlinear model.
Based on R's nls (non-linear least squares) function.

Arguments:
f -- model function, takes two argumetns the first a list of parameters
that are to be estimated, and an x value.
y -- sequence of dependent data
x -- sequence of independent data
start -- start values for the parameters to be estimated

Options:
:method (default :gauss-newton) other option :newton-raphson
:tol (default 1E-5)
:max-iter (default 200)

Returns: a hash-map containing the following fields:
:method -- the method used
:coefs  -- the parameter estimates
:hessian -- the estimated hessian, if available
:iterations -- the number of iterations performed
:fitted -- the fitted values of y (i.e. y-hat)
:x -- the independent data values
:y -- the dependent data values

Examples:

;; example 1
(use '(incanter core optimize datasets charts))
;; define the Michaelis-Menton model function
;; y = a + (b - a)*x/(c + x)
(defn f [theta x]
(let [[a b c] theta]
(plus a (div (mult x (minus b a)) (plus c x)))))

(def start [20 200 100])
(def data (to-matrix (get-dataset :thurstone)))
(def x (sel data :cols 1))
(def y (sel data :cols 0))
;; view the data
(def plot (scatter-plot x y))
(view plot)

(def nlm (non-linear-model f y x start))

;; example 2
(use '(incanter core optimize datasets charts))
;; Chwirut data set from NIST
(def data (to-matrix (get-dataset :chwirut)))
(def x (sel data :cols 1))
(def y (sel data :cols 0))

;; define model function: y = exp(-b1*x)/(b2+b3*x) + e
(defn f [theta x]
(let [[b1 b2 b3] theta]
(div (exp (mult (minus b1) x)) (plus b2 (mult b3 x)))))

(def plot (scatter-plot x y :legend true))
(view plot)

;; the newton-raphson algorithm fails to converge to the correct solution
;; using first set of start values from NIST, but the default gauss-newton
;; algorith converges to the correct solution.

(def start1 [0.1 0.01 0.02])
(add-lines plot x (f start1 x))
(def nlm1 (non-linear-model f y x start1))

;; both algorithms converges with second set of start values from NIST
(def start2 [0.15 0.008 0.010])
(add-lines plot x (f start2 x))
(def nlm2 (non-linear-model f y x start2))

"
([f y x start & options]
(let [opts (when options (apply assoc {} options))
method (or (:method opts) :guass-newton) ;; other option is :newton-raphson
tol (or (:tol opts) 1E-5)
max-iter (or (:max-iter opts) 200)
nls (if (= method :newton-raphson)
(nls-newton-raphson f (gradient f start) (hessian f start) start x y :tol tol :max-iter max-iter)
(nls-gauss-newton f start x y :tol tol :max-iter max-iter))
fitted (map #(f (:theta nls) %) x)]
{:method method
:coefs (:theta nls)
:hessian (:hessian nls)
:iterations (:iterations nls)