1 i
;;; -*- mode: lisp -*-
3 ;;; Copyright (c) 2008--, by A.J. Rossini <blindglobe@gmail.com>
4 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
5 ;;; Since 1991, ANSI was finally finished. Modified to match ANSI
8 ;; matrix-matlisp -- extended functions for matrix and linear algebra
15 (defpackage :lisp-stat-matrix-matlisp
21 (in-package :lisp-stat-matrix-matlisp
)
23 ;; Might check with features that we've (pushnew :matlisp *features*)
24 ;; and then we can do #+matlisp as needed.
27 ;;; SWEEP Operator, MATLISP implementation.
30 (defun make-sweep-front (x y w n p mode has_w x_mean result
)
36 mode is: real, complex
41 (declare (fixnum n p mode has_w
))
56 (has-w (if (/= 0 has_w
) t nil
))
58 (declare (long-float val dxi dyi dv dw sum_w dxik dxjk dyj
59 dx_meani dx_meanj dy_mean
))
61 (setf x_data
(compound-data-seq x
))
62 (setf result_data
(compound-data-seq result
))
69 (setf dyi
(makedouble (aref y i
)))
71 (setf dw
(makedouble (aref w i
)))
73 (setf dyi
(* dyi dw
)))
75 (if (not has-w
) (setf sum_w
(float n
0.0)))
76 (if (<= sum_w
0.0) (error "non positive sum of weights"))
77 (setf dy_mean
(/ val sum_w
))
79 ;; find the column means
85 (setf dxi
(makedouble (aref x_data
(+ (* p i
) j
))))
87 (setf dw
(makedouble (aref w i
)))
88 (setf dxi
(* dxi dw
)))
90 (setf (aref x_mean j
) (/ val sum_w
)))
92 ;; put 1/sum_w in topleft, means on left, minus means on top
93 (setf (aref result_data
0) (/ 1.0 sum_w
))
96 (setf dxi
(makedouble (aref x_mean i
)))
97 (setf (aref result_data
(+ i
1)) (- dxi
))
98 (setf (aref result_data
(* (+ i
1) (+ p
2))) dxi
))
99 (setf (aref result_data
(+ p
1)) (- dy_mean
))
100 (setf (aref result_data
(* (+ p
1) (+ p
2))) dy_mean
)
102 ;; put sums of adjusted cross products in body
110 (setf dxik
(makedouble (aref x_data
(+ (* p k
) i
))))
111 (setf dxjk
(makedouble (aref x_data
(+ (* p k
) j
))))
112 (setf dx_meani
(makedouble (aref x_mean i
)))
113 (setf dx_meanj
(makedouble (aref x_mean j
)))
114 (setf dv
(* (- dxik dx_meani
) (- dxjk dx_meanj
)))
116 (setf dw
(makedouble (aref w k
)))
119 (setf (aref result_data
(+ (* (+ i
1) (+ p
2)) (+ j
1))) val
)
120 (setf (aref result_data
(+ (* (+ j
1) (+ p
2)) (+ i
1))) val
))
124 (setf dxik
(makedouble (aref x_data
(+ (* p j
) i
))))
125 (setf dyj
(makedouble (aref y j
)))
126 (setf dx_meani
(makedouble (aref x_mean i
)))
127 (setf dv
(* (- dxik dx_meani
) (- dyj dy_mean
)))
129 (setf dw
(makedouble (aref w j
)))
132 (setf (aref result_data
(+ (* (+ i
1) (+ p
2)) (+ p
1))) val
)
133 (setf (aref result_data
(+ (* (+ p
1) (+ p
2)) (+ i
1))) val
))
137 (setf dyj
(makedouble (aref y j
)))
138 (setf dv
(* (- dyj dy_mean
) (- dyj dy_mean
)))
140 (setf dw
(makedouble (aref w j
)))
143 (setf (aref result_data
(+ (* (+ p
1) (+ p
2)) (+ p
1))) val
)))
145 ;;; FIXME: use matlisp
146 (defun sweep-in-place-front (a rows cols mode k tol
)
147 "Sweep algorithm for linear regression."
148 (declare (long-float tol
))
149 (declare (fixnum rows cols mode k
))
157 (declare (long-float pivot aij aik akj akk
))
159 (if (> mode RE
) (error "not supported for complex data yet"))
160 (if (or (< k
0) (>= k rows
) (>= k cols
)) (error "index out of range"))
162 (setf tol
(max tol machine-epsilon
))
163 (setf data
(compound-data-seq a
))
165 (setf pivot
(makedouble (aref data
(+ (* cols k
) k
))))
168 ((or (> pivot tol
) (< pivot
(- tol
)))
173 (when (and (/= i k
) (/= j k
))
174 (setf aij
(makedouble (aref data
(+ (* cols i
) j
))))
175 (setf aik
(makedouble (aref data
(+ (* cols i
) k
))))
176 (setf akj
(makedouble (aref data
(+ (* cols k
) j
))))
177 (setf aij
(- aij
(/ (* aik akj
) pivot
)))
178 (setf (aref data
(+ (* cols i
) j
)) aij
))))
182 (setf aik
(makedouble (aref data
(+ (* cols i
) k
))))
184 (setf aik
(/ aik pivot
))
185 (setf (aref data
(+ (* cols i
) k
)) aik
)))
189 (setf akj
(makedouble (aref data
(+ (* cols k
) j
))))
191 (setf akj
(- (/ akj pivot
)))
192 (setf (aref data
(+ (* cols k
) j
)) akj
)))
194 (setf akk
(/ 1.0 pivot
))
195 (setf (aref data
(+ (* cols k
) k
)) akk
)
199 ;; FIXME: use matlisp
200 (defun make-sweep-matrix (x y
&optional w
)
201 "Args: (x y &optional weights)
202 X is matrix, Y and WEIGHTS are sequences. Returns the sweep matrix of the
203 (weighted) regression of Y on X"
206 (if w
(check-sequence w
))
207 (let ((n (num-rows x
))
209 (if (/= n
(length y
)) (error "dimensions do not match"))
210 (if (and w
(/= n
(length w
))) (error "dimensions do not match"))
211 (let ((mode (max (la-data-mode x
)
213 (if w
(la-data-mode w
) 0)))
214 (result (make-array (list (+ p
2) (+ p
2))))
215 (x-mean (make-array p
))
216 (y (coerce y
'vector
))
217 (w (if w
(coerce w
'vector
)))
219 (make-sweep-front x y w n p mode has-w x-mean result
)
222 (defun sweep-in-place (a k tol
)
226 (let ((rows (num-rows a
))
228 (mode (la-data-mode a
)))
229 (let ((swept (sweep-in-place-front a rows cols mode k tol
)))
230 (if (/= 0 swept
) t nil
))))
232 (defun sweep-operator (a columns
&optional tolerances
)
233 "Args: (a indices &optional tolerances)
235 A is a matrix, INDICES a sequence of the column indices to be
236 swept. Returns a list of the swept result and the list of the columns
237 actually swept. (See MULTREG documentation.) If supplied, TOLERANCES
238 should be a list of real numbers the same length as INDICES. An index
239 will only be swept if its pivot element is larger than the
240 corresponding element of TOLERANCES."
242 ;; Why we should use generics/methods....
243 (if (not (typep a
'matlisp
:matrix-real
))
244 (error "a is not a matrix"))
245 (if (and (not (typep columns
'sequence
))
246 (not (typep (setf columns
(list columns
)) 'sequence
)))
247 (error "columns not coerceable to sequence"))
249 (if (and (not (typep tolerances
'sequence
))
250 (not (setf tolerances
(list tolerances
))))
251 (error "tolerances not coercable to sequence.")))
254 (check-fixnum columns
)
255 (if tolerances
(check-real tolerances
))
258 (result (copy-array a
))
260 (columns (coerce columns
'list
) (cdr columns
))
261 (tolerances (if (consp tolerances
) (coerce tolerances
'list
))
262 (if (consp tolerances
) (cdr tolerances
))))
263 ((null columns
) (list result swept-columns
))
264 (let ((col (first columns
))
265 (tol (if (consp tolerances
) (first tolerances
) tol
)))
266 (if (sweep-in-place result col tol
)
267 (setf swept-columns
(cons col swept-columns
))))))