home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / mac / science / xlspstr1.sit / nonlin.lsp < prev    next >
Lisp/Scheme  |  1990-08-29  |  11KB  |  280 lines
  1. ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
  2. ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
  3. ;;;; You may give out copies of this software; for conditions see the file
  4. ;;;; COPYING included with this distribution.
  5.  
  6. (provide "nonlin")
  7.  
  8. (require "regression")
  9.  
  10. ;;;;
  11. ;;;;
  12. ;;;; Nonlinear Regression Model Prototype
  13. ;;;;
  14. ;;;;
  15.  
  16. (defproto nreg-model-proto 
  17.           '(mean-function theta-hat epsilon count-limit verbose)
  18.           '()
  19.           regression-model-proto)
  20.  
  21. (defun nreg-model (mean-function y theta 
  22.                                  &key 
  23.                                  (epsilon .0001) 
  24.                                  (print t) 
  25.                                  (count-limit 20)
  26.                                  parameter-names
  27.                                  response-name
  28.                                  case-labels
  29.                                  weights
  30.                                  included 
  31.                                  (verbose print))
  32. "Args: (mean-function y theta &key (epsilon .0001) (count-limit 20) 
  33.                       (print t) parameter-names response-name case-labels
  34.                       weights included (vetbose print))
  35. Fits nonlinear regression model with MEAN-FUNCTION and response Y using initial
  36. parameter guess THETA. Returns model object."
  37.   (let ((m (send nreg-model-proto :new)))
  38.     (send m :mean-function mean-function)
  39.     (send m :y y)
  40.     (send m :new-initial-guess theta)
  41.     (send m :epsilon epsilon)
  42.     (send m :count-limit count-limit)
  43.     (send m :parameter-names parameter-names)
  44.     (send m :response-name response-name)
  45.     (send m :case-labels case-labels)
  46.     (send m :weights weights)
  47.     (send m :included included)
  48.     (send m :verbose verbose)
  49.     (if print (send m :display))
  50.     m))
  51.  
  52. (defmeth nreg-model-proto :save ()
  53. "Message args: ()
  54. Returns an expression that will reconstruct the regression model."
  55.   `(nreg-model ',(send self :mean-function)
  56.                ',(send self :y)
  57.                ',(send self :coef-estimates)
  58.                :epsilon ',(send self :epsilon)
  59.                :count-limit ',(send self :count-limit)
  60.                :predictor-names ',(send self :predictor-names)
  61.                :response-name ',(send self :response-name)
  62.                :case-labels ',(send self :case-labels)
  63.                :weights ',(send self :weights)
  64.                :included ',(send self :included)
  65.                :verbose ',(send self :verbose)))
  66.  
  67. ;;;
  68. ;;; Computing Method
  69. ;;;
  70.               
  71. (defmeth nreg-model-proto :compute ()
  72. "Message args: ()
  73. Recomputes the estimates. For internal use by other messages"
  74.   (let* ((y (send self :y))
  75.          (weights (send self :weights))
  76.          (inc (if-else (send self :included) 1 0))
  77.          (w (if weights (* inc weights) inc)))
  78.     (setf (slot-value 'theta-hat)
  79.           (nlreg (send self :mean-function) 
  80.                  y 
  81.                  (slot-value 'theta-hat) 
  82.                  (send self :epsilon)
  83.                  (send self :count-limit)
  84.                  w
  85.                  (send self :verbose)))
  86.     (setf (slot-value 'x) 
  87.           (funcall (make-jacobian (slot-value 'mean-function) 
  88.                                   (length (slot-value 'theta-hat))) 
  89.                    (slot-value 'theta-hat)))
  90.     (setf (slot-value 'intercept) nil)
  91.     (call-next-method)
  92.     (let ((r (send self :residuals)))
  93.       (setf (slot-value 'residual-sum-of-squares)
  94.         (sum (* inc r r))))))
  95.  
  96. ;;;
  97. ;;; Slot Accessors and Mutators
  98. ;;;
  99.  
  100. (defmeth nreg-model-proto :new-initial-guess (guess)
  101. "Message args: (guess)
  102. Sets a new initial uess for parmeters."
  103.   (setf (slot-value 'theta-hat) guess)
  104.   (send self :needs-computing t))
  105.  
  106. (defmeth nreg-model-proto :theta-hat ()
  107. "Message args: ()
  108. Returns current parameter estimate."
  109.   (if (send self :needs-computing) (send self :compute))
  110.   (coerce (slot-value 'theta-hat) 'list))
  111.  
  112. (defmeth nreg-model-proto :mean-function (&optional f)
  113. "Message args: (&optional f)
  114. With no argument returns the mean function as supplied to m. With an 
  115. argument F sets the mean function of m to F and recomputes the
  116. estimates."
  117.   (when (and f (functionp f))
  118.         (setf (slot-value 'mean-function) f)
  119.         (send self :needs-computing t))
  120.   (slot-value 'mean-function))
  121.  
  122. (defmeth nreg-model-proto :epsilon (&optional eps)
  123. "Message args: (&optional eps)
  124. With no argument returns the tolerance as supplied to m. With an argument
  125. EPS sets the tolerance of m to EPS and recomputes the estimates."
  126.   (when (and eps (numberp eps))
  127.         (setf (slot-value 'epsilon) eps)
  128.         (send self :needs-computing t))
  129.    (slot-value 'epsilon))
  130.  
  131. (defmeth nreg-model-proto :count-limit (&optional count)
  132. "Message args: (&optional new-count)
  133. With no argument returns the iteration count limit as supplied to m. With
  134. an argument COUNT sets the limit to COUNT and recomputes the
  135. estimates."
  136.   (when (and count (numberp count))
  137.         (setf (slot-value 'count-limit) count)
  138.         (send self :needs-computing t))
  139.   (slot-value 'count-limit))
  140.  
  141. (defmeth nreg-model-proto :parameter-names (&optional (names nil set))
  142. "Method args: (&optional names)
  143. Sets or returns parameter names."
  144.   (if set (setf (slot-value 'predictor-names) names))
  145.   (let ((p-names (slot-value 'predictor-names))
  146.         (p (length (slot-value 'theta-hat))))
  147.     (if (not (and p-names (= p (length p-names))))
  148.         (setf (slot-value 'predictor-names)
  149.               (mapcar #'(lambda (a) (format nil "Parameter ~a" a)) 
  150.                       (iseq 0 (- p 1))))))
  151.   (slot-value 'predictor-names))
  152.  
  153. (defmeth nreg-model-proto :verbose (&optional (val nil set))
  154. "Method args: (&optional val)
  155. Sets or retrieves verbose setting. If T iteration info is printed during
  156. optimization."
  157.   (if set (setf (slot-value 'verbose) val))
  158.   (slot-value 'verbose))
  159.  
  160. ;;;
  161. ;;; Overrides for Linear Regression Methods
  162. ;;;
  163.  
  164. (defmeth nreg-model-proto :x ()
  165. "Message args: ()
  166. Returns the Jacobian matrix at theta-hat."
  167.    (call-next-method))
  168.  
  169. (defmeth nreg-model-proto :intercept (&rest args)
  170. "Message args: ()
  171. Always returns nil. (For compatibility with linear regression.)"
  172.   nil)
  173.  
  174. (defmeth nreg-model-proto :fit-values () 
  175. "Message args: ()
  176. Returns the fitted values for the model."
  177.   (coerce (funcall (send self :mean-function) (send self :theta-hat)) 
  178.           'list))
  179.  
  180. (defmeth nreg-model-proto :coef-estimates (&optional guess)
  181. "Message args: (&optional guess)
  182. With no argument returns the current parameter estimate. With an 
  183. argument GUESS takes it as a new initial guess and recomputes
  184. the estimates."
  185.   (if guess (send self :new-initial-guess guess)) 
  186.   (send self :theta-hat))
  187.  
  188. (defmeth nreg-model-proto :predictor-names () (send self :parameter-names))
  189.  
  190. ;;;;
  191. ;;;;
  192. ;;;; Linear Regression Coefficients
  193. ;;;;
  194. ;;;;
  195.  
  196. (defun regression-coefficients (x y &key (intercept T) weights)
  197. "Args: (x y &key (intercept T) weights)
  198. Returns the coefficients of the regression of the sequence Y on the columns of
  199. the matrix X."
  200.   (let* ((m (if weights
  201.                 (make-sweep-matrix x y weights)
  202.                 (make-sweep-matrix x y)))
  203.          (n (array-dimension x 1)))
  204.     (coerce (compound-data-seq
  205.          (if intercept
  206.          (select (car (sweep-operator m (iseq 1 n)))
  207.              (1+ n)
  208.              (iseq 0 n))
  209.               (select (car (sweep-operator m (iseq 0 n)))
  210.              (1+ n)
  211.              (iseq 1 n))))
  212.             'vector)))
  213.  
  214. ;;;;
  215. ;;;;
  216. ;;;; Nonlinear Regression Functions
  217. ;;;;
  218. ;;;;
  219. (defun nlreg1 (f j y initial-beta epsilon count-limit weights verbose)
  220. "Args: (mean-function jacobian y initial-beta 
  221.         epsilon count-limit weights verbose)
  222. MEAN-FUNCTION returns the mean response vector for a given parameter vector.
  223. JACOBIAN returns the jacobian of MEAN-FUNCTION at a given parameter vector.
  224. Y is the observed response vector. Returns the estimated parameter vector 
  225. obtained by a Gauss-Newton algorithm with backtracking that continues until
  226. the COUNT-LIMIT is reached or no component of the parameter vector changes
  227. by more than EPSILON."
  228.   (labels ((rss (beta)                   ; residual sum of squares
  229.                 (let ((res (- y (funcall f beta)))) 
  230.                   (sum (if weights (* res res weights) (* res res)))))
  231.            (next-beta (beta delta rss)   ; next beta by backtracking
  232.                       (do* ((lambda 1 (/ lambda 2))
  233.                             (new-rss (rss (+ beta delta)) 
  234.                                      (rss (+ beta (* lambda delta)))))
  235.                            ((or (< new-rss rss) (< lambda .0001))
  236.                             (if (>= lambda .0001)
  237.                                 (+ beta (* lambda delta))
  238.                                 beta)))))
  239.     (do* ((delbeta (regression-coefficients
  240.                     (funcall j initial-beta)
  241.                     (- y (funcall f initial-beta))
  242.                     :intercept nil
  243.                     :weights weights)
  244.                    (regression-coefficients
  245.                     (funcall j beta)
  246.                     (- y (funcall f beta))
  247.                     :intercept nil
  248.                     :weights weights))
  249.           (beta initial-beta (next-beta beta delbeta rss))
  250.           (rss (rss beta) (rss beta))
  251.           (count 0 (1+ count)))
  252.          ((or (> count count-limit) (> epsilon (max (abs delbeta))))
  253.           (if (and verbose (> count count-limit))
  254.               (format t "Iteration limit exceeded.~%"))
  255.           beta)
  256.          (if verbose (format t "Residual sum of squares: ~10g~%" rss)))))
  257.  
  258. (defun make-jacobian (f n)
  259. "Args: (f n)
  260. F is a function of an N-vector. Returns a function that approximates the
  261. jacobian function iof F by a symmetric difference."
  262.   (let* ((h .0001)
  263.          (del (* h (column-list (identity-matrix n)))))
  264.     #'(lambda (b) 
  265.       (let ((b+ (mapcar #'(lambda (x) (+ b x)) del))
  266.             (b- (mapcar #'(lambda (x) (- b x)) del)))
  267.         (apply #'bind-columns (/ (- (mapcar f b+) (mapcar f b-)) (* 2 h)))))))
  268.  
  269. (defun nlreg (f y guess &optional 
  270.                 (epsilon .0001) (count-limit 20) weights verbose)
  271. "Args: (mean-function y guess &optional 
  272.         (epsilon .0001) (count-limit 20) weights verbose)
  273. MEAN-FUNCTION returns the mean response vector for a given parameter vector.
  274. Y is the observed response vector. Returns the estimated parameter vector 
  275. obtained by a Gauss-Newton algorithm that continues until the ITERATION-LIMIT
  276. is reached or no component of the parameter vector changes by more than
  277. EPSILON. The jacobian of MEAN-FUNCTION is approximated by a symmetric difference."
  278.   (nlreg1 f (make-jacobian f (length guess)) y guess 
  279.           epsilon count-limit weights verbose))
  280.