1.2.3-SNAPSHOT Arrow_down_16x16

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
:gradient -- the estimated gradient
:hessian -- the estimated hessian, if available
:iterations -- the number of iterations performed
:fitted -- the fitted values of y (i.e. y-hat)
:rss -- the residual sum-of-squares
: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))
(add-lines plot x (:fitted nlm))


;; example 2
(use '(incanter core optimize datasets charts))
;; Chwirut data set from NIST
;; http://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Chwirut1.dat
(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))
(add-lines plot x (:fitted nlm1))

;; 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))
(add-lines plot x (:fitted nlm2))


0 Examples top

Log in to add / edit an example.

See Also top

Log in to add a see also.

Plus_12x12 Minus_12x12 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
    :gradient  -- the estimated gradient
    :hessian -- the estimated hessian, if available
    :iterations -- the number of iterations performed
    :fitted -- the fitted values of y (i.e. y-hat)
    :rss -- the residual sum-of-squares
    :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))
    (add-lines plot x (:fitted nlm))


    ;; example 2
    (use '(incanter core optimize datasets charts))
    ;; Chwirut data set from NIST
    ;; http://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Chwirut1.dat
    (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))
    (add-lines plot x (:fitted nlm1))

    ;; 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))
    (add-lines plot x (:fitted nlm2))


"
  ([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)
       :gradient (:gradient nls)
       :hessian (:hessian nls)
       :iterations (:iterations nls)
       :rss (:rss nls)
       :fitted fitted
       :x x
       :y y})))
Vars in incanter.optimize/non-linear-model: defn let
Used in 0 other vars

Comments top

No comments for non-linear-model. Log in to add a comment.