1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
(setf eps 1e-10)
(setf x '#2a((1 2)(3 4)(5 6)(7 8)))
(setf y '#2a((9 10 11) (12 13 14)))
(setf xv '#(15 16))
(setf yv '#(18 19))
(defun mm (x y)
(let* ((m (array-dimension x 0))
(k (array-dimension x 1))
(n (array-dimension y 1))
(z (make-array (list m n) :initial-element 0)))
(dotimes (i m)
(dotimes (j n)
(dotimes (l k)
(incf (aref z i j) (* (aref x i l) (aref y l j))))))
z))
(defun mv (x y)
(let* ((m (array-dimension x 0))
(n (array-dimension x 1))
(z (make-array m :initial-element 0)))
(dotimes (i m)
(dotimes (j n)
(incf (aref z i) (* (aref x i j) (aref y j)))))
z))
(defun vm (x y)
(let* ((m (array-dimension y 0))
(n (array-dimension y 1))
(z (make-array n :initial-element 0)))
(dotimes (i m)
(dotimes (j n)
(incf (aref z j) (* (aref x i) (aref y i j)))))
z))
(defun ip (x y &optional (conj t))
(sum (* x (if conj (conjugate y) y))))
(defun test-report (name val)
;(format t "~&~a: ~9t~f~%" name val)
(check #'< val eps)
)
(test-report "MATMULT"
(let ((cx (complex x (* 2 x)))
(cy (complex y (* 3 y)))
(cxv (complex xv (* 4 xv)))
(cyv (complex yv (* 5 yv))))
(max (abs (- (matmult x y) (mm x y)))
(abs (- (matmult cx cy) (mm cx cy)))
(abs (- (matmult x yv) (mv x yv)))
(abs (- (matmult cx cyv) (mv cx cyv)))
(abs (- (matmult xv y) (vm xv y)))
(abs (- (matmult cxv cy) (vm cxv cy)))
(abs (- (matmult xv yv) (ip xv yv)))
(abs (- (matmult cxv cyv) (ip cxv cyv nil))))))
(test-report "CROSS-PRODUCT"
(max (abs (- (cross-product x) (mm (transpose x) x)))
(let ((cx (complex x (* 2 x))))
(abs (- (cross-product cx)
(mm (transpose cx) (conjugate cx)))))
(let ((cx (complex x (* 2 x))))
(abs (- (cross-product cx nil)
(mm (transpose cx) cx))))))
(test-report "INNER-PRODUCT"
(let ((cxv (complex xv (* 4 xv)))
(cyv (complex yv (* 5 yv))))
(max (abs (- (inner-product xv yv) (ip xv yv)))
(abs (- (inner-product cxv cyv) (ip cxv cyv)))
(abs (- (inner-product cxv cyv nil) (ip cxv cyv nil))))))
(let* ((x (list (iseq 1 4) (^ (iseq 1 4) 2)))
(y (^ (iseq 1 4) 3))
(w (iseq 4 1))
(m1 (regression-model x y :print nil))
(mw (regression-model x y :print nil :weights w))
(xm (apply #'bind-columns (repeat 1 4) x))
(xmt (transpose xm))
(wm (diagonal w))
(b1 (matmult (inverse (matmult xmt xm)) xmt y))
(bw (matmult (inverse (matmult xmt wm xm)) xmt wm y)))
(test-report "REGRESS"
(max (abs (- (send m1 :coef-estimates) b1))
(abs (- (send mw :coef-estimates) bw)))))
|
