So: you want to perform ordinary least squares linear regression.
I’ve needed this code a bunch of times, and the matrix library actually
suggests that something like this might be useful, but I my implementation
technique is “transliterate from wikipedia,” and I have a strong suspicion that
a numerical methods person would cry if they saw this code. On the other hand,
I’m pretty sure that it’s ballpark correct; my test cases suggest that I’m not
doing anything wildly incorrect.
Is there a better version of this floating around somewhere? If not, will
seeing this crappy code spur someone else to do a better job?
John
#lang typed/racket
(require math/matrix)
(provide XY->Beta
XYBeta-rms-error)
;; given a matrix X of inputs and a column-vector Y of responses,
;; compute Beta, the column vector mapping X to y-hat (minimizing
;; the sum of the squares of the residuals), which can be used
;; to estimate the response based on an input
(: XY->Beta ((Matrix Number) (Matrix Number) -> (Matrix Number)))
(define (XY->Beta X Y)
(unless (= (matrix-num-rows X) (matrix-num-rows Y))
(raise-argument-error
'XY->B "Y matrix with same # of rows as X matrix"
1 X Y
))
(unless (= (matrix-num-cols Y) 1)
(raise-argument-error
'XY->B "Y matrix with 1 column"
1 X Y))
;; this might be a really slow or inaccurate way to compute this....
(matrix*
(matrix* (matrix-inverse (matrix* (matrix-transpose X) X))
(matrix-transpose X))
Y))
;; given a matrix X of inputs, a column-vector Y of responses,
;; and a column vector B intended to map inputs to outputs,
;; return the rms error.
(: XYBeta-rms-error
((Matrix Number) (Matrix Number) (Matrix Number) -> Real))
(define (XYBeta-rms-error X y Beta)
(define epsilon
(matrix- y (matrix* X Beta)))
(match (matrix->list (matrix* (matrix-transpose epsilon) epsilon))
[(list s)
(cond [(= s 0.0) 0]
[else (sqrt (/ (cast s Positive-Real)
(cast (matrix-num-rows y) Positive-Real)))])]))
(module+ test
(require typed/rackunit
math/array)
;; example:
;; suppose y = 2x_1 - 8x_2, with small perturbations
;; 23 -4 | 78.9 (error +0.9)
;; 4 -2 | 23.2 (error -0.8)
;; 0 2 | -16.3 (error -0.3)
(define X
(matrix [[23 -4]
[4 -2]
[0 2]]))
(define Y
(matrix [[78.9]
[23.2]
[-16.3]]))
(define estimated-b (array->list (XY->Beta X Y)))
(check-= (cast (car estimated-b) Real) 2 0.2)
(check-= (cast (cadr estimated-b) Real) -8 0.2)
(check-= (XYBeta-rms-error X Y (col-matrix [2 -8]))
(sqrt (* 1/3 (+ 0.81 0.64 0.09)))
1e-5)
;; example from Wikipedia page on ordinary least squares:
(define heights
'(1.47 1.50 1.52 1.55 1.57 1.60 1.63 1.65 1.68 1.70 1.73 1.75 1.78 1.80
1.83))
(define sqr-hts
(map (λ ([x : Real]) (* x x)) heights))
(define units
(for/list : (Listof Real) ([h (in-list heights)]) 1.0))
(define X2
(matrix-transpose (list*->array (list units heights sqr-hts) real?)))
(define weights
(col-matrix [52.21 53.12 54.48 55.84 57.20 58.57 59.93 61.29 63.11 64.47
66.28 68.10 69.92 72.19 74.46]))
;; wikipedia lists expected result:
;; β 1 {\displaystyle \beta _{1}} \beta _{1} 128.8128 16.3083
7.8986 0.0000
;; β 2 {\displaystyle \beta _{2}} \beta _{2} –143.1620 19.8332
–7.2183 0.0000
;; β 3 {\displaystyle \beta _{3}} \beta _{3} 61.9603 6.0084 10.3122
0.0000
(check-equal?
(matrix->list (XY->Beta X2 weights))
(list 128.8128035850009
-143.16202287195165
61.960325444600926))
)
--
You received this message because you are subscribed to the Google Groups
"Racket Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to racket-users+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
