more tests to track down multiple regression failures
[CommonLispStat.git] / ls-demo.lisp
blob486f01f8d2f9de953989db3c00a85fd63ae477c2
1 ;;; -*- mode: lisp -*-
2 ;;; Copyright (c) 2007, by A.J. Rossini <blindglobe@gmail.com>
3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
4 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
6 ;;; Time-stamp: <2008-03-09 09:21:50 user>
7 ;;; Creation: <2007-01-01 09:21:50 user> WRONG
8 ;;; File: ls-demo.lisp
9 ;;; Author: AJ Rossini <blindglobe@gmail.com>
10 ;;; Copyright: (c) 2007, AJ Rossini. BSD.
11 ;;; Purpose: demonstrations of how one might use CLS.
13 ;;; What is this talk of 'release'? Klingons do not make software
14 ;;; 'releases'. Our software 'escapes', leaving a bloody trail of
15 ;;; designers and quality assurance people in its wake.
17 (load "init.lisp")
18 ;; init needs to be more like the asdf-loader for lisp-stat, though it
19 ;; is pretty close.
21 ;;; non-rigorous check for exports.
22 ;;; This is generally not how I expect it to be used.
24 (in-package :cl-user)
25 (lisp-stat:binomial-quant 0.95 3 0.4) ;;; 3
26 (lisp-stat:binomial-quant 0 3 0.4) ;;; 0
27 (lisp-stat:normal-rand 20) ;;; 20 numbers :-)
29 ;;;; THIS is how I expect it to be used, either with work in ls-user,
30 ;;;; or a cloned package similar to ls-user.
32 (in-package :ls-user)
34 ;;;; Matrix algebra.
36 (chol-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
37 ;; (#2A((1.7888543819998317 0.0 0.0)
38 ;; (1.6770509831248424 0.11180339887498929 0.0)
39 ;; (2.23606797749979 2.23606797749979 3.332000937312528e-8))
40 ;; 5.000000000000003)
42 (defvar my-chol-decomp-test (chol-decomp #2A((2 3 4) (1 2 4) (2 4 5))))
43 my-chol-decomp-test
44 (nth 0 my-chol-decomp-test)
45 (nth 1 my-chol-decomp-test)
48 (lu-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
49 ;; (#2A((2.0 3.0 4.0) (1.0 1.0 1.0) (0.5 0.5 1.5)) #(0 2 2) -1.0 NIL)
51 (lu-solve
52 (lu-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
53 #(2 3 4))
54 ;; #(-2.333333333333333 1.3333333333333335 0.6666666666666666)
56 (inverse #2A((2 3 4) (1 2 4) (2 4 5)))
57 ;; #2A((2.0 -0.33333333333333326 -1.3333333333333335)
58 ;; (-1.0 -0.6666666666666666 1.3333333333333333)
59 ;; (0.0 0.6666666666666666 -0.3333333333333333))
61 (sv-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
62 ;; (#2A((-0.5536537653489974 0.34181191712789266 -0.7593629708013371)
63 ;; (-0.4653437312661058 -0.8832095891230851 -0.05827549615722014)
64 ;; (-0.6905959164998124 0.3211003503429828 0.6480523475178517))
65 ;; #(9.699290438141343 0.8971681569301373 0.3447525123483081)
66 ;; #2A((-0.30454218417339873 0.49334669582252344 -0.8147779426198863)
67 ;; (-0.5520024849987308 0.6057035911404464 0.5730762743603965)
68 ;; (-0.7762392122368734 -0.6242853493399995 -0.08786630745236332))
69 ;; T)
71 (qr-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
72 ;; (#2A((-0.6666666666666665 0.7453559924999298 5.551115123125783e-17)
73 ;; (-0.3333333333333333 -0.2981423969999719 -0.894427190999916)
74 ;; (-0.6666666666666666 -0.5962847939999439 0.44721359549995787))
75 ;; #2A((-3.0 -5.333333333333334 -7.333333333333332)
76 ;; (0.0 -0.7453559924999292 -1.1925695879998877)
77 ;; (0.0 0.0 -1.3416407864998738)))
79 (rcondest #2A((2 3 4) (1 2 4) (2 4 5)))
80 ;; 6.8157451e7
81 ;;; CURRENTLY FAILS!!
83 (eigen #2A((2 3 4) (1 2 4) (2 4 5)))
84 ;; (#(10.656854249492381 -0.6568542494923802 -0.9999999999999996)
85 ;; (#(0.4999999999999998 0.4999999999999997 0.7071067811865475)
86 ;; #(-0.49999999999999856 -0.5000000000000011 0.7071067811865474)
87 ;; #(0.7071067811865483 -0.7071067811865466 -1.2560739669470215e-15))
88 ;; NIL)
90 (spline #(1.0 1.2 1.3 1.8 2.1 2.5)
91 #(1.2 2.0 2.1 2.0 1.1 2.8) :xvals 6)
92 ;; ((1.0 1.3 1.6 1.9 2.2 2.5)
93 ;; (1.2 2.1 2.2750696543866313 1.6465231041904045 1.2186576148879609 2.8))
95 ;;; using KERNEL-SMOOTH-FRONT, not KERNEL-SMOOTH-CPORT
96 (kernel-smooth #(1.0 1.2 1.3 1.8 2.1 2.5)
97 #(1.2 2.0 2.1 2.0 1.1 2.8) :xvals 5)
98 ;; ((1.0 1.375 1.75 2.125 2.5)
99 ;; (1.6603277642110226 1.9471748095239771 1.7938127405752287
100 ;; 1.5871511322219498 2.518194783156392))
102 (kernel-dens #(1.0 1.2 2.5 2.1 1.8 1.2) :xvals 5)
103 ;; ((1.0 1.375 1.75 2.125 2.5)
104 ;; (0.7224150453621405 0.5820045548233707 0.38216411702854214
105 ;; 0.4829822708587095 0.3485939156929503))
107 (fft #(1.0 1.2 2.5 2.1 1.8))
108 ;; #(#C(1.0 0.0) #C(1.2 0.0) #C(2.5 0.0) #C(2.1 0.0) #C(1.8 0.0))
110 (lowess #(1.0 1.2 2.5 2.1 1.8 1.2) #(1.2 2.0 2.1 2.0 1.1 2.8))
111 ;; (#(1.0 1.2 1.2 1.8 2.1 2.5))
115 ;;;; Special functions
117 ;; Log-gamma function
119 (log-gamma 3.4) ;;1.0923280596789584
123 ;;;; Probability functions
125 ;; looking at these a bit more, perhaps a more CLOSy style is needed, i.e.
126 ;; (quantile :list-or-cons loc :type type (one of 'empirical 'normal 'cauchy, etc...))
127 ;; similar for the cdf, density, and rand.
128 ;; Probably worth figuring out how to add a new distribution
129 ;; efficiently, i.e. by keeping some kind of list.
131 ;; Normal distribution
133 (normal-quant 0.95) ;;1.6448536279366268
134 (normal-cdf 1.3) ;;0.9031995154143897
135 (normal-dens 1.3) ;;0.17136859204780736
136 (normal-rand 2) ;;(-0.40502015f0 -0.8091404f0)
138 (bivnorm-cdf 0.2 0.4 0.6) ;;0.4736873734160288
140 ;; Cauchy distribution
142 (cauchy-quant 0.95) ;;6.313751514675031
143 (cauchy-cdf 1.3) ;;0.7912855998398473
144 (cauchy-dens 1.3) ;;0.1183308127104695
145 (cauchy-rand 2) ;;(-1.06224644160405 -0.4524695943939537)
147 ;; Gamma distribution
149 (gamma-quant 0.95 4.3) ;;8.178692439291645
150 (gamma-cdf 1.3 4.3) ;;0.028895150986674906
151 (gamma-dens 1.3 4.3) ;;0.0731517686447374
152 (gamma-rand 2 4.3) ;;(2.454918912880936 4.081365384357454)
154 ;; Chi-square distribution
156 (chisq-quant 0.95 3) ;;7.814727903379012
157 (chisq-cdf 1 5) ;;0.03743422675631789
158 (chisq-dens 1 5) ;;0.08065690818083521
159 (chisq-rand 2 4) ;;(1.968535826180572 2.9988646156942997)
161 ;; Beta distribution
163 (beta-quant 0.95 3 2) ;;0.9023885371149876
164 (beta-cdf 0.4 2 2.4) ;;0.4247997418541529
165 (beta-dens 0.4 2 2.4) ;;1.5964741858913518
166 (beta-rand 2 2 2.4) ;;(0.8014897077282279 0.6516371997922659)
168 ;; t distribution
170 (t-quant 0.95 3) ;;2.35336343484194
171 (t-cdf 1 2.3) ;;0.794733624298342
172 (t-dens 1 2.3) ;;0.1978163816318102
173 (t-rand 2 2.3) ;;(-0.34303672776089306 -1.142505872436518)
175 ;; F distribution
177 (f-quant 0.95 3 5) ;;5.409451318117459
178 (f-cdf 1 3.2 5.4) ;;0.5347130905510765
179 (f-dens 1 3.2 5.4) ;;0.37551128864591415
180 (f-rand 2 3 2) ;;(0.7939093442091963 0.07442694152491144)
182 ;; Poisson distribution
184 (poisson-quant 0.95 3.2) ;;6
185 (poisson-cdf 1 3.2) ;;0.17120125672252395
186 (poisson-pmf 1 3.2) ;;0.13043905274097067
187 (poisson-rand 5 3.2) ;;(2 1 2 0 3)
189 ;; Binomial distribution
191 (binomial-quant 0.95 3 0.4) ;;; DOESN'T RETURN
192 (binomial-quant 0 3 0.4) ;;; -2147483648
193 (binomial-cdf 1 3 0.4) ;;0.6479999999965776
194 (binomial-pmf 1 3 0.4) ;;0.4320000000226171
195 (binomial-rand 5 3 0.4) ;;(2 2 0 1 2)
197 ;;;; OBJECT SYSTEM
199 (in-package :ls-user)
200 (defproto *test-proto*)
201 *test-proto*
202 (defmeth *test-proto* :make-data (&rest args) nil)
204 (defvar my-proto-instance nil)
205 (setf my-proto-instance (send *test-proto* :new))
206 (send *test-proto* :own-slots)
207 (send *test-proto* :has-slot 'proto-name)
208 (send *test-proto* :has-slot 'PROTO-NAME)
209 (send *test-proto* :has-slot 'make-data)
210 (send *test-proto* :has-slot 'MAKE-DATA)
211 (send *test-proto* :has-method 'make-data)
212 (send *test-proto* :has-method 'MAKE-DATA)
215 (defproto2 *test-proto3* (list) (list) (list) "test doc" t)
216 (defproto2 *test-proto4*)
217 *test-proto2*
218 (defmeth *test-proto* :make-data (&rest args) nil)
220 (defvar my-proto-instance nil)
221 (setf my-proto-instance (send *test-proto* :new))
222 (send *test-proto* :own-slots)
223 (send *test-proto* :has-slot 'proto-name)
224 (send *test-proto* :has-slot 'PROTO-NAME)
227 ;;;; Testing
229 (in-package :lisp-stat-unittests)
230 (testsuites)
231 (print-tests)
232 (run-tests)
233 (last-test-status)
234 ;;(failures)
236 (describe (run-tests :suite 'lisp-stat-testsupport))
237 (describe (run-tests :suite 'lisp-stat-testsupport2))
239 (testsuite-tests 'lisp-stat)
240 (run-tests :suite 'lisp-stat)
241 (describe (run-tests :suite 'lisp-stat))
243 (run-tests :suite 'lisp-stat-probdistn)
244 (describe (run-tests :suite 'lisp-stat-probdistn))
245 (run-tests :suite 'lisp-stat-spec-fns)
246 (describe (run-tests :suite 'lisp-stat-spec-fns))
248 (find-testsuite 'lisp-stat-lin-alg)
249 (testsuite-tests 'lisp-stat-lin-alg)
250 (run-tests :suite 'lisp-stat-lin-alg)
251 (describe (run-tests :suite 'lisp-stat-lin-alg))
253 ;;;; Data Analysis test
255 (in-package :ls-user)
257 ;; LispStat 1 approach to variables
259 (progn
260 (def iron (list 61 175 111 124 130 173 169 169 160 224 257 333 199))
261 iron
262 (def aluminum (list 13 21 24 23 64 38 33 61 39 71 112 88 54))
263 aluminum
264 (def absorbtion (list 4 18 14 18 26 26 21 30 28 36 65 62 40))
265 absorbtion
267 ;; LispStat 1 approach to data frames... (list of lists).
269 (DEF DIABETES
270 (QUOTE ((80 97 105 90 90 86 100 85 97 97 91 87 78 90 86 80 90 99 85 90 90 88 95 90 92 74 98 100 86 98 70 99 75 90 85 99 100 78 106 98 102 90 94 80 93 86 85 96 88 87 94 93 86 86 96 86 89 83 98 100 110 88 100 80 89 91 96 95 82 84 90 100 86 93 107 112 94 93 93 90 99 93 85 89 96 111 107 114 101 108 112 105 103 99 102 110 102 96 95 112 110 92 104 75 92 92 92 93 112 88 114 103 300 303 125 280 216 190 151 303 173 203 195 140 151 275 260 149 233 146 124 213 330 123 130 120 138 188 339 265 353 180 213 328 346)
271 (356 289 319 356 323 381 350 301 379 296 353 306 290 371 312 393 364 359 296 345 378 304 347 327 386 365 365 352 325 321 360 336 352 353 373 376 367 335 396 277 378 360 291 269 318 328 334 356 291 360 313 306 319 349 332 323 323 351 478 398 426 439 429 333 472 436 418 391 390 416 413 385 393 376 403 414 426 364 391 356 398 393 425 318 465 558 503 540 469 486 568 527 537 466 599 477 472 456 517 503 522 476 472 455 442 541 580 472 562 423 643 533 1468 1487 714 1470 1113 972 854 1364 832 967 920 613 857 1373 1133 849 1183 847 538 1001 1520 557 670 636 741 958 1354 1263 1428 923 1025 1246 1568)
272 (124 117 143 199 240 157 221 186 142 131 221 178 136 200 208 202 152 185 116 123 136 134 184 192 279 228 145 172 179 222 134 143 169 263 174 134 182 241 128 222 165 282 94 121 73 106 118 112 157 292 200 220 144 109 151 158 73 81 151 122 117 208 201 131 162 148 130 137 375 146 344 192 115 195 267 281 213 156 221 199 76 490 143 73 237 748 320 188 607 297 232 480 622 287 266 124 297 326 564 408 325 433 180 392 109 313 132 285 139 212 155 120 28 23 232 54 81 87 76 42 102 138 160 131 145 45 118 159 73 103 460 42 13 130 44 314 219 100 10 83 41 77 29 124 15)
273 (3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 2 2 3 2 2 3 3 3 3 2 3 3 3 3 3 2 3 3 3 3 3 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1))))
276 (DEF DLABS (QUOTE ("GLUFAST" "GLUTEST" "INSTEST" "CCLASS")))
277 (format t "loaded data.~%")
280 ;; Simple univariate variable-specific descriptions.
281 (fivnum absorbtion)
282 (median absorbtion)
283 (sort-data absorbtion)
284 (rank absorbtion)
285 (standard-deviation absorbtion)
286 (interquartile-range absorbtion)
288 (lisp-stat-matrix::bind-columns aluminum iron)
289 (bind-columns aluminum iron)
290 (apply #'bind-columns (list aluminum iron))
291 (lisp-stat-matrix::bind-columns #2a((1 2)(3 4)) #(5 6))
292 (bind-columns #2a((1 2)(3 4)) #(5 6))
295 (defvar fit1 nil)
296 (setf fit1 (regression-model absorbtion iron))
297 (send fit1 :display)
298 (send fit1 :residuals)
300 iron
301 (defvar fit1a nil)
302 (setf fit1a (regression-model absorbtion iron :print nil))
303 (send fit1a :doc)
304 (setf (send fit1a :doc) "this") ;; FIXME: this is a more naturualo
305 (send fit1a :x)
306 (send fit1a :y)
307 (send fit1a :compute)
308 (send fit1a :sweep-matrix)
309 (send fit1a :basis)
310 (send fit1a :residuals)
311 (send fit1a :display)
313 #+nil(progn
314 (array-dimension #2A ((1)) 0)
315 ;; more tests
318 ;;; FIXME: need to get multiple-linear regression working -- clearly
319 ;;; simple linear is working above!
320 (defvar m nil "holding variable.")
321 (def m (regression-model (list iron aluminum) absorbtion :print nil))
322 (send m :compute)
323 (send m :sweep-matrix)
324 (format t "~%~A~%" (send m :sweep-matrix))
325 (send m :display) ;; ERROR...
326 (def m (regression-model (bind-columns iron aluminum) absorbtion))
327 (send m :help)
328 (send m :help 'display)
329 (send m :plot-residuals)
332 (typep aluminum 'sequence)
333 (typep iron 'sequence)
334 (matrixp iron)
336 *variables*
338 (variables)
339 (undef 'iron)
340 (variables)
343 ;;; Example array calcs
345 #+nil(progn
346 (functionp #'and)
347 (= (array-dimensions #2A((2 3 3 ) (2 4 4)))
348 (array-dimensions #2A((2 3 3 ) (2 5 4))))
349 (reduce #'and (= (array-dimensions #2A((2 3) (2 4)))
350 (array-dimensions #2A((2 3 3 ) (2 5 4)))))
352 (defvar my-t-ar nil)
353 (setf my-t-ar #3A(((2 3) (2 2) (2 1))
354 ((2 3) (2 2) (2 1))))
355 (defvar my-t-ar2 nil)
356 (setf my-t-ar2 #2A((1 2 3 4)
357 (5 6 7 8)))
359 (array-dimensions my-t-ar)
360 (array-dimensions my-t-ar2)
362 (aref my-t-ar2 1 2) ;; GOOD
363 (aref my-t-ar2 (list 1 2)) ;; BAD
364 (apply #'aref my-t-ar2 (list 1 2)) ;; GOOD
365 ;; For extracting multiple elements
366 (mapcar #'(lambda (x) (apply #'aref my-t-ar2 x))
367 (list (list 0 0) (list 0 1)))
370 (aref my-t-ar 1 2 1)
371 (aref my-t-ar 1 2 1)
372 (aref my-t-ar 1 1 0)
374 (array-dimensions #3A(((2 3) (2 2) (2 1))
375 ((2 3) (2 2) (2 1))))
377 (reduce #'and (= #(2 3) #(2 4)))
378 (= #(2 3) #(2 3))
380 (let ((a #2A((2 3 3 ) (2 5 4)))
381 (b #2A((2 3 3 ) (2 5 4))))
382 (let ((a-rank (array-rank a))
383 (a-dim (array-dimensions a))
384 (a-b-elt-eq (loop for i in 0 to (aref a-dim 0)
385 for j in 0 to (aref a-dim 1)
386 collect (numerical= (apply #'aref a (list i j))
387 (apply #'aref b (list i j))))))
388 (every #'(lambda (x) x) a-b-elt-eq))))
390 (every #'(lambda (x) x) (list T T T))
391 (every #'(lambda (x) x) (list T T nil))
393 (and T T)
394 (mapcar #'(lambda (&rest args) (and args))
395 (list (= #(2 3) #(2 4))))
396 (= #(2 3) #(2 3))
398 ;;; examples of using CLEM
400 (in-package :clem-user)
402 (defvar m1 (make-instance 'double-float-matrix :rows 10 :cols 5))
405 (defvar m2 (make-instance 'number-matrix :rows 10 :cols 5))
409 ;;; Not defined but documented? Actually somewhere in clem/print.lisp
410 (setf *matrix-print-row-limit* 2)
411 (setf *matrix-print-col-limit* 2)
413 (print m2)
416 (mat-log m2)
417 (mat-abs m2)
418 (min m2)
419 (max m2)
421 (setf (mref m2 1 1) 5)
423 (setf (mref m2 0 0) 5)
427 ;;; prop list demo
430 (defvar tplist (list :this 'though :that 'there :thee 'tony))
431 (setf (symbol-plist tplist) (list :this 'though :that 'there :thee 'tony))
432 (get 'tplist :THIS)
433 tplist
434 (defvar tlist (list :this 'though :that 'there :thee 'tony))
435 (setf tlist (list :this 'though :that 'there :thee 'tony))
436 (listp tlist)
437 (getf tlist :THIS)
439 ;;; CL-SDL demos
442 (clc:clc-require :sdl-demos)
443 ;;(sdl-test:start) ; locks up SBCL.?
444 ;;; where <n> is 2-11, 16, :solar-system, :vertex-arrays,
445 (nehe:run-tutorial 2)