-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathlinfit.lisp
More file actions
executable file
·188 lines (170 loc) · 7.29 KB
/
linfit.lisp
File metadata and controls
executable file
·188 lines (170 loc) · 7.29 KB
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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
;; linfit.lisp -- Linear regression with Peter Stetson's iteration
;;
;; DM/RAL 11/10
;; ---------------------------------------------------------------
(defpackage #:com.ral.linfit
(:use #:common-lisp)
(:export
#:wmnsd
#:regression
#:regress-fixed-slope
#:regress-fixed-intercept))
(in-package #:com.ral.linfit)
(defun stetson-refinement (wts devs sigma &key (alpha 2) (beta 2))
;; take the resids and wts and produce new wts
(map 'vector (lambda (dev wt)
(/ wt (+ 1 (expt (/ (abs dev) alpha sigma) beta))))
devs wts))
(defun regression (xs ys wts &key (alpha 2) (beta 2) (eps 1d-3))
;; return y-weighted mean and slope though data centroid along with stdev of fit.
;; iterate until the relative improvement in overall chisq
;; is less than eps. Be careful not to set eps below the sqrt machine precision.
(let* ((wts (if (numberp wts)
(make-array (length ys)
:initial-element wts)
wts))
;; (nel (length xs))
;; (dof (- nel 2))
)
(um:nlet iter ((swts wts)
(sigprev nil)
(niter 0))
(let* ((twt (matrix:sum swts))
(neff (/ (* twt twt) (matrix:<*> swts swts)))
(xmn (/ (matrix:<*> xs swts) twt))
(xsmrm (matrix:- xs xmn))
(ywmn (/ (matrix:<*> ys swts) twt))
(ysmrm (matrix:- ys ywmn))
(slope (/ (matrix:<*> ysmrm xsmrm swts)
(matrix:<*> xsmrm xsmrm swts)))
(devs (matrix:- ysmrm (matrix:* slope xsmrm)))
(wsigma (sqrt (/ (matrix:<*> swts devs devs)
twt (/ (- neff 2) neff)))))
(if (and (< niter 100)
(plusp wsigma)
(or (null sigprev)
(> (abs (- wsigma sigprev)) (* eps sigprev))))
(go-iter (stetson-refinement wts devs wsigma
:alpha alpha
:beta beta)
wsigma
(1+ niter))
;; else
(progn
(format t "~%linfit:regression niter = ~A" niter)
(values xmn ywmn slope wsigma niter neff (/ wsigma (sqrt neff)) ) ))))))
(defun regress-fixed-slope (xs ys wts slope &key (alpha 2) (beta 2) (eps 1d-3))
;; return weighted best intercept b in y = slope * x + b, along with stdev of fit
;; iterate until the relative improvement in overall chisq
;; is less than eps. Be careful not to set eps below the sqrt machine precision.
(let* ((wts (if (numberp wts)
(make-array (length ys)
:initial-element wts)
wts))
(diffs (map 'vector (lambda (x y)
(- y (* slope x)))
xs ys))
(nel (length xs))
(dof (1- nel)))
(um:nlet iter ((swts wts)
(sigprev nil)
(niter 0))
(let* ((twt (matrix:sum swts))
(b (/ (matrix:<*> diffs swts) twt))
(devs (matrix:- diffs b))
(wsigma (sqrt (/ (matrix:<*> swts devs devs)
twt (/ dof nel)))))
(if (and (< niter 100)
(plusp wsigma)
(or (null sigprev)
(> (abs (- wsigma sigprev)) (* eps sigprev))))
(go-iter (stetson-refinement wts devs wsigma
:alpha alpha
:beta beta)
wsigma
(1+ niter))
;; else
(progn
(format t "~%linfit:regress-fixed-slope niter = ~A" niter)
(values b wsigma niter) ))))))
(defun regress-fixed-intercept (xs ys wts intercept &key (alpha 2) (beta 2) (eps 1d-3))
;; return weighted best intercept b in y = slope * x + b, along with stdev of fit
;; iterate until the relative improvement in overall chisq
;; is less than eps. Be careful not to set eps below the sqrt machine precision.
(let* ((wts (if (numberp wts)
(make-array (length ys)
:initial-element wts)
wts))
(nel (length xs))
(dof (1- nel)))
(um:nlet iter ((swts wts)
(sigprev nil)
(niter 0))
(let* ((twt (matrix:sum swts))
(slope (/ (matrix:<*> (matrix:- ys intercept) xs swts)
(matrix:<*> xs xs swts)))
(devs (matrix:- ys (matrix:+ (matrix:* slope xs) intercept)))
(wsigma (sqrt (/ (matrix:<*> swts devs devs)
twt (/ dof nel)))))
(if (and (< niter 100)
(plusp wsigma)
(or (null sigprev)
(> (abs (- wsigma sigprev)) (* eps sigprev))))
(go-iter (stetson-refinement wts devs wsigma
:alpha alpha
:beta beta)
wsigma
(1+ niter))
;; else
(progn
(format t "~%linfit:regress-fixed-intercept niter = ~A" niter)
(values slope wsigma niter) ))))))
(defun wmnsd (ys wts &key (alpha 2) (beta 2) (eps 1d-3))
;; return weighted mean and estimated standard deviation of data
;; iterate until the relative improvement in overall chisq
;; is less than eps. Be careful not to set eps below the sqrt machine precision.
(let* ((wts (if (numberp wts)
(make-array (length ys)
:initial-element wts)
wts))
(nel (length ys))
(dof (1- nel)))
(um:nlet iter ((swts wts)
(sigprev 0)
(niter 0))
(let* ((twt (matrix:sum swts))
(ywmn (/ (matrix:<*> ys swts) twt))
(devs (matrix:- ys ywmn))
(wsigma (sqrt (/ (matrix:<*> swts devs devs)
twt (/ dof nel)))))
(cond ((or (<= (abs (- wsigma sigprev)) (* eps sigprev))
(> niter 100))
(format t "~%linfit:wmnsd niter = ~A" niter)
(values ywmn wsigma
(/ wsigma (sqrt twt)) ;; est S.D. of mean
twt
niter) )
(t
(go-iter (stetson-refinement wts devs wsigma
:alpha alpha
:beta beta)
wsigma
(1+ niter)))
)))
))
#|
(let* ((xs #(1 2 3)))
(matrix:+ 3 xs))
(let* ((ys (vm:unoise 1000 1.0)))
(wmnsd ys 1))
(let* ((ys (vm:gnoise 1000 :sd 1.0)))
(wmnsd ys 1 :alpha 2 :beta 2))
(let* ((ys (vm:gnoise 1000 :sd 1.0)))
(let ((mn (vm:mean ys)))
(list mn (vm:stdev ys mn))))
(let* ((ys (matrix:+ (vm:gnoise 1000) (vm:unoise 1000 2))))
(wmnsd ys 1))
(let* ((ys (matrix:+ (vm:gnoise 1000) (vm:unoise 1000 2))))
(let ((mn (vm:mean ys)))
(list mn (vm:stdev ys mn))))
|#