A bit more philosophy and task for visualization, in response to Peter S. Still...
[CommonLispStat.git] / examples / ls-demo.lisp
blobd757c5442170d81128222d03a4047d1a406857e5
1 ;;; -*- mode: lisp -*-
2 ;;; Copyright (c) 2006-2008, 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: <2011-12-06 16:04:17 tony>
7 ;;; Creation: <2009-09-17 22:19:31 tony> (sometime earlier, but serious now)
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.
18 ;; Start in the usual user-package for loading
19 (in-package :cl-user)
21 ;; You've configured CLS, right?
22 (asdf:oos 'asdf:load-op 'cls)
24 ;; Go somewhere so that you have the functions available.
25 (in-package :ls-user)
28 ;; we'll be loading from directories in the CLS homedir, so we want to
29 ;; make it easier to reach.
30 (defparameter *my-cls-homedir*
31 "/home/tony/sandbox/CLS.git/" ; <- value with trailing
32 ; directory separator
33 "documentation: change this to localize") ; <- doc
34 ;; so
35 (concatenate 'string *my-cls-homedir* "Data/example.csv")
36 ;; implies
37 (defun localized-pathto (x)
38 "return a string denoting the complete path.
39 FIXME: UNIX-centric (though might work on Mac OSX). Might want to
40 return a pathspec, not a string/namespec"
41 (check-type x string)
42 (concatenate 'string *my-cls-homedir* x))
45 ;;; DataFrames
47 (defparameter *df-test*
48 (make-instance 'dataframe-array
49 :storage #2A (('a "test0" 0 0d0)
50 ('b "test1" 1 1d0)
51 ('c "test2" 2 2d0)
52 ('d "test3" 3 3d0)
53 ('e "test4" 4 4d0))
54 :doc "test reality"
55 :case-labels (list "0" "1" "2" "3" "4")
56 :var-labels (list "symb-var" "strng-var" "int-var" "dfloat-var")
57 :var-types (list 'symbol 'string 'integer 'double-float)))
59 *df-test* ; but with SBCL, ints become floats?
60 (caselabels *df-test*)
61 (varlabels *df-test*)
62 (vartypes *df-test*) ;; currently errors
64 (setf (xref *df-test* 0 0) -1d0) ; for dataframes, we might want to do
65 ; some type checking to prevent what
66 ; I just did!
68 (setf (xref *df-test* 0 0) (quote 'a)) ; so that we can restore the
69 ; quoted value.
70 *df-test*
72 ;; Making from arrays and matrix-likes
73 (make-dataframe #2A((1 2 3 4 5)
74 (10 20 30 40 50))) ;; another error...
75 (make-dataframe (rand 4 3)) ;; another error
77 ;;; HERE#1
78 ;;; == READ DATA
80 ;;; read in a CSV dataframe...
82 ;; a better approach is:
83 (asdf:oos 'asdf:load-op 'rsm-string)
84 (listoflist->array
85 (rsm.string:file->string-table
86 (localized-pathto "Data/example-mixed.csv")
87 :delims ","))
89 (rsm.string:file->number-table
90 (localized-pathto "Data/example-numeric.csv")
91 :delims ",")
95 (rsm.string:file->number-table
96 (localized-pathto "Data/R-chickwts.csv")
97 :delims ",")
98 (rsm.string:file->string-table
99 (localized-pathto "Data/R-chickwts.csv")
100 :delims ",")
102 (defparameter *my-df-2*
103 (make-instance 'dataframe-array
104 :storage
105 (listoflist->array
106 (rsm.string:file->string-table
107 (localized-pathto "Data/example-mixed.csv")))
108 :doc "This is an interesting dataframe-array"))
109 ;; *my-df-2*
111 (defparameter *my-df-3*
112 (make-instance 'dataframe-array
113 :storage
114 (listoflist->array
115 (transpose-listoflist
116 (rsm.string:file->number-table
117 (localized-pathto "Data/example-numeric.csv"))))
118 :doc "This is an interesting dataframe-array"))
119 ;; *my-df-3*
122 (defparameter *my-df-4*
123 (make-instance 'dataframe-array
124 :storage
125 (listoflist->array
126 (rsm.string:file->number-table
127 (localized-pathto "Data/R-chickwts.csv")
128 :delims ","))
129 :doc "This is an interesting dataframe-array that currently fails"))
130 ;; *my-df-4*
132 (aref (dataset *my-df-4*) 0 1)
134 (xref (dataset *my-df-4*) 0 1) ;; preferred to use the x* tools.
137 (defparameter *my-df-5*
138 (make-instance 'dataframe-array
139 :storage
140 (listoflist->array
141 (transpose-listoflist
142 (rsm.string:file->number-table
143 (localized-pathto "Data/R-swiss.csv"))))
144 :doc "This is an interesting dataframe-array that currently fails"))
145 ;; *my-df-5*
149 (defparameter *mat-1*
150 (make-matrix 3 3
151 :initial-contents #2A((2d0 3d0 4d0) (3d0 2d0 4d0) (4d0 4d0 5d0))))
153 (defparameter *mat-1*
154 (make-matrix 3 3
155 :initial-contents #2A((2d0 3d0 -4d0)
156 (3d0 2d0 -4d0)
157 (4d0 4d0 -5d0))))
158 (mref *mat-1* 2 0)
160 (xref *mat-1* 2 0) ;; fails, but should work.
164 (defparameter *mat-2*
165 (let ((m (rand 3 3)))
166 (m* m (transpose m)))) ;; fails, but should work.
168 (defparameter *mat-2*
169 (let ((m (rand 3 3)))
170 (m* m (transpose-matrix m)))) ;; works, it's transpose-matrix
172 (axpy 100.0d0 *mat-2* (eye 3 3))
174 (potrf (copy *mat-2*)) ;; factor
175 (potri (copy *mat-2*)) ;; invert
176 (minv-cholesky (copy *mat-2*))
177 (m* (minv-cholesky (copy *mat-2*)) *mat-2*)
179 (defparameter *mat-3*
180 (make-matrix
182 :initial-contents '((16d0 13d0 12d0)
183 (13d0 22d0 7d0)
184 (12d0 7d0 17d0))))
186 (potrf (copy *mat-3*)) ;; factor
189 *mat-3* =>
190 #<LA-SIMPLE-MATRIX-DOUBLE 3 x 3
191 16.0 13.0 12.0
192 13.0 22.0 7.0
193 12.0 7.0 17.0>
195 (potrf (copy *mat-3*)) =>
196 (#<LA-SIMPLE-MATRIX-DOUBLE 3 x 3
197 4.0 3.25 3.0
198 13.0 3.3819373146171707 -0.8131433980500301
199 12.0 7.0 2.7090215603069034>
200 "U" NIL)
202 ;; and compare with...
204 > testm <- matrix(data=c(16,13,12,13,22,7,12,7,17),nrow=3)
205 > chol(testm)
206 [,1] [,2] [,3]
207 [1,] 4 3.250000 3.0000000
208 [2,] 0 3.381937 -0.8131434
209 [3,] 0 0.000000 2.7090216
212 ;; which suggests that the major difference is that R zero's out the
213 ;; appropriate terms, and that CLS does not.
217 (potri (copy *mat-2*)) ;; invert
218 (minv-cholesky (copy *mat-2*))
219 (m* (minv-cholesky (copy *mat-2*)) *mat-2*)
223 (lu-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
224 ;; => (#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)
225 (lu-solve
226 (lu-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
227 #(2 3 4))
228 ;; => #(-2.333333333333333 1.3333333333333335 0.6666666666666666)
230 (getrf
231 (make-matrix 3 3
232 :initial-contents #2A((2d0 3d0 4d0) (1d0 2d0 4d0) (2d0 4d0 5d0))))
234 #| => ; so not so good for the vector, but matrix matches.
235 (#<LA-SIMPLE-MATRIX-DOUBLE 3 x 3
236 2.0 3.0 4.0
237 1.0 1.0 1.0
238 0.5 0.5 1.5>
239 #<FNV-INT32 (3) 1 3 3> NIL)
242 (msolve-lu
243 (make-matrix 3 3
244 :initial-contents #2A((2d0 3d0 4d0)
245 (1d0 2d0 4d0)
246 (2d0 4d0 5d0)))
247 (make-vector 3 :type :column
248 :initial-contents '((2d0)
249 (3d0)
250 (4d0))))
252 #| =>
253 #<LA-SIMPLE-VECTOR-DOUBLE (3 x 1)
254 -2.3333333333333335
255 1.3333333333333335
256 0.6666666666666666>
261 ;;; LU common applications
263 (defun minv-lu (a)
264 "invert A using LU Factorization"
265 (let ((a-fac (getrf (copy a))))
266 (first (getri (first a-fac) (second a-fac)))))
268 #+nil (progn
269 (let ((m1 (rand 3 3)))
270 (m* m1 (minv-lu m1))))
272 (defun msolve-lu (a b)
273 "Compute `x1' solving `A x = b', with LU factorization."
274 (let ((a-fac (getrf (copy a))))
275 (first (getrs (first a-fac) b (second a-fac)))))
279 ;; (inverse #2A((2 3 4) (1 2 4) (2 4 5)))
280 ;; #2A((2.0 -0.33333333333333326 -1.3333333333333335)
281 ;; (-1.0 -0.6666666666666666 1.3333333333333333)
282 ;; (0.0 0.6666666666666666 -0.3333333333333333))
284 (minv-lu
285 (make-matrix
287 :initial-contents #2A((2d0 3d0 4d0)
288 (1d0 2d0 4d0)
289 (2d0 4d0 5d0))))
293 #<LA-SIMPLE-MATRIX-DOUBLE 3 x 3
294 2.0 -0.3333333333333333 -1.3333333333333333
295 -1.0 -0.6666666666666666 1.3333333333333333
296 0.0 0.6666666666666666 -0.3333333333333333>
298 ;; so is correct.
302 ;;;;;HERE#2
304 (factorize
305 (make-matrix 3 3
306 :initial-contents #2A((2d0 3d0 4d0)
307 (1d0 2d0 4d0)
308 (2d0 4d0 5d0)))
309 :by :svd)
311 ;; (sv-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
312 ;; (#2A((-0.5536537653489974 0.34181191712789266 -0.7593629708013371)
313 ;; (-0.4653437312661058 -0.8832095891230851 -0.05827549615722014)
314 ;; (-0.6905959164998124 0.3211003503429828 0.6480523475178517))
315 ;; #(9.699290438141343 0.8971681569301373 0.3447525123483081)
316 ;; #2A((-0.30454218417339873 0.49334669582252344 -0.8147779426198863)
317 ;; (-0.5520024849987308 0.6057035911404464 0.5730762743603965)
318 ;; (-0.7762392122368734 -0.6242853493399995 -0.08786630745236332))
319 ;; T)
323 (qr-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
324 ;; (#2A((-0.6666666666666665 0.7453559924999298 5.551115123125783e-17)
325 ;; (-0.3333333333333333 -0.2981423969999719 -0.894427190999916)
326 ;; (-0.6666666666666666 -0.5962847939999439 0.44721359549995787))
327 ;; #2A((-3.0 -5.333333333333334 -7.333333333333332)
328 ;; (0.0 -0.7453559924999292 -1.1925695879998877)
329 ;; (0.0 0.0 -1.3416407864998738)))
331 (rcondest #2A((2 3 4) (1 2 4) (2 4 5)))
332 ;; 6.8157451e7
333 ;;; CURRENTLY FAILS!!
335 (eigen #2A((2 3 4) (1 2 4) (2 4 5)))
336 ;; (#(10.656854249492381 -0.6568542494923802 -0.9999999999999996)
337 ;; (#(0.4999999999999998 0.4999999999999997 0.7071067811865475)
338 ;; #(-0.49999999999999856 -0.5000000000000011 0.7071067811865474)
339 ;; #(0.7071067811865483 -0.7071067811865466 -1.2560739669470215e-15))
340 ;; NIL)
342 (spline #(1.0 1.2 1.3 1.8 2.1 2.5)
343 #(1.2 2.0 2.1 2.0 1.1 2.8) :xvals 6)
344 ;; ((1.0 1.3 1.6 1.9 2.2 2.5)
345 ;; (1.2 2.1 2.2750696543866313 1.6465231041904045 1.2186576148879609 2.8))
347 ;;; using KERNEL-SMOOTH-FRONT, not KERNEL-SMOOTH-CPORT
348 (kernel-smooth #(1.0 1.2 1.3 1.8 2.1 2.5)
349 #(1.2 2.0 2.1 2.0 1.1 2.8) :xvals 5)
350 ;; ((1.0 1.375 1.75 2.125 2.5)
351 ;; (1.6603277642110226 1.9471748095239771 1.7938127405752287
352 ;; 1.5871511322219498 2.518194783156392))
354 (kernel-dens #(1.0 1.2 2.5 2.1 1.8 1.2) :xvals 5)
355 ;; ((1.0 1.375 1.75 2.125 2.5)
356 ;; (0.7224150453621405 0.5820045548233707 0.38216411702854214
357 ;; 0.4829822708587095 0.3485939156929503))
359 (fft #(1.0 1.2 2.5 2.1 1.8))
360 ;; #(#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))
362 (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))
363 ;; (#(1.0 1.2 1.2 1.8 2.1 2.5))
367 ;;;; Special functions
369 ;; Log-gamma function
371 (log-gamma 3.4) ;;1.0923280596789584
375 ;;;; Probability functions
377 ;; looking at these a bit more, perhaps a more CLOSy style is needed, i.e.
378 ;; (quantile :list-or-cons loc :type type (one of 'empirical 'normal 'cauchy, etc...))
379 ;; similar for the cdf, density, and rand.
380 ;; Probably worth figuring out how to add a new distribution
381 ;; efficiently, i.e. by keeping some kind of list.
383 ;; Normal distribution
385 (normal-quant 0.95) ;;1.6448536279366268
386 (normal-cdf 1.3) ;;0.9031995154143897
387 (normal-dens 1.3) ;;0.17136859204780736
388 (normal-rand 2) ;;(-0.40502015f0 -0.8091404f0)
390 (bivnorm-cdf 0.2 0.4 0.6) ;;0.4736873734160288
392 ;; Cauchy distribution
394 (cauchy-quant 0.95) ;;6.313751514675031
395 (cauchy-cdf 1.3) ;;0.7912855998398473
396 (cauchy-dens 1.3) ;;0.1183308127104695
397 (cauchy-rand 2) ;;(-1.06224644160405 -0.4524695943939537)
399 ;; Gamma distribution
401 (gamma-quant 0.95 4.3) ;;8.178692439291645
402 (gamma-cdf 1.3 4.3) ;;0.028895150986674906
403 (gamma-dens 1.3 4.3) ;;0.0731517686447374
404 (gamma-rand 2 4.3) ;;(2.454918912880936 4.081365384357454)
406 ;; Chi-square distribution
408 (chisq-quant 0.95 3) ;;7.814727903379012
409 (chisq-cdf 1 5) ;;0.03743422675631789
410 (chisq-dens 1 5) ;;0.08065690818083521
411 (chisq-rand 2 4) ;;(1.968535826180572 2.9988646156942997)
413 ;; Beta distribution
415 (beta-quant 0.95 3 2) ;;0.9023885371149876
416 (beta-cdf 0.4 2 2.4) ;;0.4247997418541529
417 (beta-dens 0.4 2 2.4) ;;1.5964741858913518
418 (beta-rand 2 2 2.4) ;;(0.8014897077282279 0.6516371997922659)
420 ;; t distribution
422 (t-quant 0.95 3) ;;2.35336343484194
423 (t-cdf 1 2.3) ;;0.794733624298342
424 (t-dens 1 2.3) ;;0.1978163816318102
425 (t-rand 2 2.3) ;;(-0.34303672776089306 -1.142505872436518)
427 ;; F distribution
429 (f-quant 0.95 3 5) ;;5.409451318117459
430 (f-cdf 1 3.2 5.4) ;;0.5347130905510765
431 (f-dens 1 3.2 5.4) ;;0.37551128864591415
432 (f-rand 2 3 2) ;;(0.7939093442091963 0.07442694152491144)
434 ;; Poisson distribution
436 (poisson-quant 0.95 3.2) ;;6
437 (poisson-cdf 1 3.2) ;;0.17120125672252395
438 (poisson-pmf 1 3.2) ;;0.13043905274097067
439 (poisson-rand 5 3.2) ;;(2 1 2 0 3)
441 ;; Binomial distribution
443 (binomial-quant 0.95 3 0.4) ;;; DOESN'T RETURN
444 (binomial-quant 0 3 0.4) ;;; -2147483648
445 (binomial-cdf 1 3 0.4) ;;0.6479999999965776
446 (binomial-pmf 1 3 0.4) ;;0.4320000000226171
447 (binomial-rand 5 3 0.4) ;;(2 2 0 1 2)
449 ;;;; OBJECT SYSTEM
451 (in-package :ls-user)
452 (defproto *test-proto*)
453 *test-proto*
454 (defmeth *test-proto* :make-data (&rest args) nil)
456 (defparameter my-proto-instance nil)
457 (setf my-proto-instance (send *test-proto* :new))
458 (send *test-proto* :own-slots)
459 (lsos::ls-object-slots *test-proto*)
460 (lsos::ls-object-methods *test-proto*)
461 (lsos::ls-object-parents *test-proto*)
462 (lsos::ls-object-preclist *test-proto*)
463 ;;; The following fail and I do not know why?
464 (send *test-proto* :has-slot 'proto-name)
465 (send *test-proto* :has-slot 'PROTO-NAME)
466 (send *test-proto* :has-slot 'make-data)
467 (send *test-proto* :has-slot 'MAKE-DATA)
468 (send *test-proto* :has-method 'make-data)
469 (send *test-proto* :has-method 'MAKE-DATA)
472 (defproto2 *test-proto3* (list) (list) (list) "test doc" t)
473 (defproto2 *test-proto4*)
474 *test-proto2*
475 (defmeth *test-proto* :make-data (&rest args) nil)
477 (defparameter my-proto-instance nil)
478 (setf my-proto-instance (send *test-proto* :new))
479 (send *test-proto* :own-slots)
480 (send *test-proto* :has-slot 'proto-name)
481 (send *test-proto* :has-slot 'PROTO-NAME)
484 ;;;; Testing
486 (in-package :lisp-stat-unittests)
487 (testsuites)
488 (print-tests)
489 (run-tests)
490 (last-test-status)
491 ;;(failures)
493 (describe (run-tests :suite 'lisp-stat-ut-testsupport))
494 (describe (run-tests :suite 'lisp-stat-ut-testsupport2))
496 (testsuite-tests 'lisp-stat-ut)
497 (run-tests :suite 'lisp-stat-ut)
498 (describe (run-tests :suite 'lisp-stat-ut))
500 (run-tests :suite 'lisp-stat-ut-probdistn)
501 (describe (run-tests :suite 'lisp-stat-ut-probdistn))
502 (run-tests :suite 'lisp-stat-ut-spec-fns)
503 (describe (run-tests :suite 'lisp-stat-ut-spec-fns))
505 (find-testsuite 'lisp-stat-ut-lin-alg)
506 (testsuite-tests 'lisp-stat-ut-lin-alg)
507 (run-tests :suite 'lisp-stat-ut-lin-alg)
508 (describe (run-tests :suite 'lisp-stat-ut-lin-alg))
510 ;;;; Data Analysis test
512 (in-package :ls-user)
514 ;; LispStat 1 approach to variables
516 (progn
517 (def iron (list 61 175 111 124 130 173 169 169 160 224 257 333 199))
518 iron
519 (def aluminum (list 13 21 24 23 64 38 33 61 39 71 112 88 54))
520 aluminum
521 (def absorbtion (list 4 18 14 18 26 26 21 30 28 36 65 62 40))
522 absorbtion
524 ;; LispStat 1 approach to data frames... (list of lists).
526 (DEF DIABETES
527 (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)
528 (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)
529 (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)
530 (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))))
533 (DEF DLABS (QUOTE ("GLUFAST" "GLUTEST" "INSTEST" "CCLASS")))
534 (format t "loaded data.~%")
535 ) ;; eval at this point.
537 ;; Simple univariate variable-specific descriptions.
538 (fivnum absorbtion)
539 (median absorbtion)
540 (sort-data absorbtion)
541 (rank absorbtion)
542 (standard-deviation absorbtion)
543 (interquartile-range absorbtion)
545 (lisp-stat-matrix::bind-columns aluminum iron)
546 (bind-columns aluminum iron)
547 (apply #'bind-columns (list aluminum iron))
548 (lisp-stat-matrix::bind-columns #2a((1 2)(3 4)) #(5 6))
549 (bind-columns #2a((1 2)(3 4)) #(5 6))
552 (defparameter fit1 nil)
553 (setf fit1 (regression-model absorbtion iron))
554 (send fit1 :display)
555 (send fit1 :residuals)
557 iron
558 (defparameter fit1a nil)
559 (setf fit1a (regression-model absorbtion iron :print nil))
560 (send fit1a :doc)
561 ;; (setf (send fit1a :doc) "this") ;; FIXME: this error...
562 (send fit1a :doc "this") ;; FIXME: this is a more natural
563 (send fit1a :doc)
564 (send fit1a :x)
565 (send fit1a :y)
566 (send fit1a :compute)
567 (send fit1a :sweep-matrix)
568 (send fit1a :basis)
569 (send fit1a :residuals)
570 (send fit1a :display)
572 #+nil(progn
573 ;; syntax example
574 (array-dimension #2A ((1)) 0)
577 ;;; FIXME: need to get multiple-linear regression working -- clearly
578 ;;; simple linear is working above!
579 (defvar m nil "holding variable.")
580 (def m (regression-model (list iron aluminum) absorbtion :print nil))
581 (send m :compute)
582 (send m :sweep-matrix)
583 (format t "~%~A~%" (send m :sweep-matrix))
585 ;; ERROR... FIX-ME!!
586 (send m :basis) ;; this should be positive?
587 (send m :coef-estimates)
589 (send m :display)
590 (def m (regression-model (bind-columns iron aluminum) absorbtion))
591 (send m :help)
592 (send m :help :display)
593 (send m :help :basis)
594 ;; No graphics! But handle the error gracefully...
595 (send m :plot-residuals)
598 (typep aluminum 'sequence)
599 (typep iron 'sequence)
600 (matrixp iron)
602 *variables*
604 (variables)
605 (undef 'iron)
606 (variables)
608 ;;; Plotting!
610 (asdf:oos 'asdf:compile-op 'cl-cairo2 :force t)
611 (asdf:oos 'asdf:load-op 'cl-cairo2)
613 ;; The above can be used to generate PDF, PS, PNG, and X11/Microsoft
614 ;; displays (the latter being a proof of concept, of limited use for
615 ;; "real work".
617 ;; and this below, as well.
618 (asdf:oos 'asdf:load-op 'cl-plplot)
620 ;;; Using R!
622 (asdf:oos 'asdf:compile-op 'rclg :force t)
623 (asdf:oos 'asdf:load-op 'rclg)
626 (in-package :rclg-user)
628 ;; rclg-init::*r-started*
630 ;;;#3 Start R within Lisp
632 (start-rclg)
633 ;; rclg-init::*r-started*
634 (rclg-init::check-stack)
635 (r "Cstack_info")
636 (defparameter *x* (r seq 1 10))
637 (defparameter *y* (r rnorm 10))
639 (r plot *x* *y*)
642 (defparameter *r-version* (r "version"))
644 ;; This is for illustrative purposes only. It is not a "good" use of rnbi.
645 ;; Really, you'll want rnbi to hold anonymous intermeditae results, like:
646 (r plot *x* (rnbi rnorm 10))
648 (r "Sys.getenv" "LD_LIBRARY_PATH")
649 (r "Sys.getenv" "LD_PRELOAD")
651 (r "ls")
652 (r ls)
653 (r "search")
655 (r "geterrmessage")
657 (r "library" "stats")
658 (r library "MASS")
659 (r "library" "Biobase")
661 (setf my.lib "Biobase")
662 my.lib
663 (r library my.lib)
665 (r "ls")
667 (r "print.default" 3)
668 (r "rnorm" 10)
670 ;; Working in the R space
672 (r assign "x" 5)
673 (r assign "x2" (list 1 2 3 5))
675 (r assign "x2" #(1 2 3 5 3 4 5))
676 (r assign "z" "y") ;; unlike the above, this assigns character data
677 (r "ls")
678 (r ls)
680 (setf my.r.x2 (r get "x2")) ;; moving data from R to CL
681 (r assign "x2" my.r.x2) ;; moving data from CL to R
683 ;; The following is not the smartest thing to do!
684 ;;(r q)
688 ;;; How might we do statistics with Common Lisp?
689 ;;; How might we work with a data.frame?
690 ;;; What could the structures be?
691 ;;; How much hinting, and of what type, should drive the data
692 ;;; analysis?
694 (defpackage :my-data-analysis-example
695 (:documentation "Example work-package for a data analysis")
696 (:use :common-lisp :lisp-stat)
697 (:export results figures report))
699 (in-package :my-data-analysis-example)
701 (defvar my-dataset1 (read-file "data/test1.lisp"))
702 ;; or
703 (defvar my-dataset2 (read-file "data/test1.csv" :type 'csv))
705 ;;; manipulate
707 (setf my-dataset2 (set-description my-datasets2
708 :dependent-variables (list of symbols)))
709 (setf my-dataset2 (set-description my-datasets2
710 :independent-variables (list of symbols)))
712 ;; the following could be true in many cases.
713 (assert
714 (list-intersection (get-description my-datasets2 :independent-variables)
715 (get-description my-datasets2 :dependent-variables)))
717 ;; but we could phrase better,i.e.
719 (get-description
720 my-datasets2
721 :predicate-list-on-variable-metadata (list (and 'independent-variables
722 'dependent-variables)))
725 ;; statistical relations re: input/output, as done above, is one
726 ;; issue, another one is getting the right approach for statistical
727 ;; typing, i.e.
728 (get-description
729 my-datasets2
730 :predicate-list-on-variable-metadata (list 'ordinal-variables))
733 ;; so we could use a set of logical ops to selection from variable
734 ;; metadata, i.e.
735 ;; and, or, not
736 ;; do we really need the simplifying extensions?
739 ;;; output to REPL
741 (report my-dataset1 :style 'five-num)
742 (report my-dataset1 :style 'univariate)
743 (report my-dataset1 :style 'bivariate)
744 (report my-dataset1 :style 'metadata)
746 ;;; to file?
748 (report my-dataset1
749 :style 'five-num
750 :format 'pdf
751 :stream (filename-as-stream "my-dataset1-5num.pdf"))
752 (report my-dataset1 :style 'univariate)
753 (report my-dataset1 :style 'bivariate)
754 (report my-dataset1 :style 'metadata)
756 ;;; so report could handle datasets... and models?
758 (report my-model :style 'formula)
759 (report my-model :style 'simulate
760 (list :parameters (:eta 5 :mu 4 :sigma (list 2 1 0.5))
761 :number-of-reps 10))
762 ;; should return a list of parameters along with range information,
763 ;; useful for auto-building the above. Note that there are 3 types
764 ;; of parameters that can be considered -- we can have values which
765 ;; define ddata, we can have values which define fixed values and some
766 ;; could be things tht we estimate.
769 (defgeneric report (object &optional style format stream)
770 (:documentation "method for reporting on data"))
772 (defmethod report ((object dataset)
773 (style report-dataset-style-type)
774 (format output-format-type)
775 ((stream *repl*) output-stream-type))
776 "dataset reporting")
779 (defmethod report ((object model)
780 (style report-model-style-type)
781 (format output-format-type)
782 ((stream *repl*) output-stream-type))
783 "model reporting")
785 (defmethod report ((object analysis-instance)
786 (style report-analysis-style-type)
787 (format output-format-type)
788 ((stream *repl*) output-stream-type))
789 "model + dataset reporting")
792 ;; parameters are just things which get filled with values, repeatedly
793 ;; with data, or by considering to need estimation.
794 (parameters my-model)
795 (parameters my-model :type 'data)
796 (parameters my-model :type 'fixed)
797 (parameters my-model :type 'estimate)
798 (parameters my-model :type '(estimate fixed))
799 (parameters my-model :list-types) ;; useful for list-based extraction
800 ;; of particular types
802 (setf my-model-data-instance
803 (compute model data :specification (list :spec 'linear-model
804 :depvar y
805 :indepvar (list x1 x2))))
806 (report my-model-data-instance)
809 ;;; So how might we use this? Probably need to consider the
810 ;;; serialization of any lisp objects generated, perhaps via some form
811 ;;; of memoization...?
812 (in-package :cl-user)
814 (my-data-analysis-example:report :type 'full)
815 (my-data-analysis-example:report :type 'summary)
816 (my-data-analysis-example:figures :type 'pdf :file "results-figs.pdf")
818 (my-data-analysis-example:report)
820 ;;; more stuff
822 (send m :display)
823 (def m (regression-model (bind-columns iron aluminum) absorbtion))
824 (send m :help)
825 (send m :help :display)
826 (send m :help :basis)
828 (send m :plot-residuals)
830 (progn
831 ;; General Lisp, there is also a need to add, remove symbols from the
832 ;; workspace/namespace. This is a fundamental skill, similar to
833 ;; stopping, which is critical.
835 ;; boundp, fboundp
836 ;; makunbound, fmakunbound
840 (progn
841 ;;; A study in array vs list access
842 (defparameter *x* (list 1 2 3))
843 (defparameter *y* #(1 2 3))
844 (defparameter *z* (list 1 (list 2 3) (list 4 5 (list 6 7)) ))
845 (length *x*)
846 (length *y*)
847 (length *z*) ; => need a means to make this 7.
848 (length (reduce #'cons *z*)) ; => not quite -- missing iterative
850 (nelts *x*)
851 (nth 1 *x*)
852 (aref *y* 1)
853 (setf (nth 1 *x*) 6)
855 (setf (aref *y* 1) 6)
859 (in-package :ls-user)
861 (progn
862 (defparameter *x* (make-vector 5 :initial-contents '((1d0 2d0 3d0 4d0 5d0))))
863 ;; estimating a mean, simple way.
864 (/ (loop for i from 0 to (- (nelts *x*) 1)
865 summing (vref *x* i))
866 (nelts *x*))
868 (defun mean (x)
869 (checktype x 'vector-like)
870 (/ (loop for i from 0 to (- (nelts *x*) 1)
871 summing (vref *x* i))
872 (nelts *x*)))
874 ;; estimating variance, Moments
875 (let ((meanx (mean *x*))
876 (n (nelts *x*)))
877 (/ (loop for i from 0 to (1- n)
878 summing (* (- (vref *x* i) meanx)
879 (- (vref *x* i) meanx)))
882 ;; estimating variance, Moments
883 (let ((meanx (mean *x*))
884 (nm1 (1- (nelts *x*))))
885 (/ (loop for i from 0 to nm1
886 summing (* (- (vref *x* i) meanx)
887 (- (vref *x* i) meanx) ))
888 nm1))
892 ;;;;;;;;;;;;;;; Data stuff
894 (progn ;; Data setup
896 ;; Making data-frames (i.e. cases (rows) by variables (columns))
897 ;; takes a bit of getting used to. For this, it is important to
898 ;; realize that we can do the following:
899 ;; #1 - consider the possibility of having a row, and transposing
900 ;; it, so the list-of-lists is: ((1 2 3 4 5)) (1 row, 5 columns)
901 ;; #2 - naturally list-of-lists: ((1)(2)(3)(4)(5)) (5 rows, 1 column)
902 ;; see src/data/listoflist.lisp for code to process this particular
903 ;; data structure.
904 (defparameter *indep-vars-1-matrix*
905 (transpose (make-matrix 1 (length iron)
906 :initial-contents
907 (list (mapcar #'(lambda (x) (coerce x 'double-float))
908 iron))))
909 "creating iron into double float, straightforward")
911 (documentation '*indep-vars-1-matrix* 'variable)
912 ;; *indep-vars-1-matrix*
914 ;; or directly:
915 (defparameter *indep-vars-1a-matrix*
916 (make-matrix (length iron) 1
917 :initial-contents
918 (mapcar #'(lambda (x) (list (coerce x 'double-float)))
919 iron)))
920 ;; *indep-vars-1a-matrix*
922 ;; and mathematically, they seem equal:
923 (m= *indep-vars-1-matrix* *indep-vars-1a-matrix*) ; => T
924 ;; but of course not completely...
925 (eql *indep-vars-1-matrix* *indep-vars-1a-matrix*) ; => NIL
926 (eq *indep-vars-1-matrix* *indep-vars-1a-matrix*) ; => NIL
928 ;; and verify...
929 (print *indep-vars-1-matrix*)
930 (print *indep-vars-1a-matrix*)
932 (documentation 'lisp-matrix:bind2 'function) ; by which we mean:
933 (documentation 'bind2 'function)
934 (bind2 *indep-vars-1-matrix* *indep-vars-1a-matrix* :by :column) ; 2 col
935 (bind2 *indep-vars-1-matrix* *indep-vars-1a-matrix* :by :row) ; 1 long col
937 ;; the weird way
938 (defparameter *indep-vars-2-matrix*
939 (transpose (make-matrix 2 (length iron)
940 :initial-contents
941 (list
942 (mapcar #'(lambda (x) (coerce x 'double-float))
943 iron)
944 (mapcar #'(lambda (x) (coerce x 'double-float))
945 aluminum)))))
946 ;; *indep-vars-2-matrix*
948 ;; the "right"? way
949 (defparameter *indep-vars-2-matrix*
950 (make-matrix (length iron) 2
951 :initial-contents
952 (mapcar #'(lambda (x y)
953 (list (coerce x 'double-float)
954 (coerce y 'double-float)))
955 iron aluminum)))
956 ;; *indep-vars-2-matrix*
959 ;; The below FAILS due to coercion issues; it just isn't lispy, it's R'y.
961 (defparameter *dep-var* (make-vector (length absorbtion)
962 :initial-contents (list absorbtion)))
964 ;; BUT below, this should be the right type.
965 (defparameter *dep-var*
966 (make-vector (length absorbtion)
967 :type :row
968 :initial-contents
969 (list
970 (mapcar #'(lambda (x) (coerce x 'double-float))
971 absorbtion))))
972 ;; *dep-var*
975 (defparameter *dep-var-int*
976 (make-vector (length absorbtion)
977 :type :row
978 :element-type 'integer
979 :initial-contents (list absorbtion)))
981 (typep *dep-var* 'matrix-like) ; => T
982 (typep *dep-var* 'vector-like) ; => T
984 (typep *indep-vars-1-matrix* 'matrix-like) ; => T
985 (typep *indep-vars-1-matrix* 'vector-like) ; => T
986 (typep *indep-vars-2-matrix* 'matrix-like) ; => T
987 (typep *indep-vars-2-matrix* 'vector-like) ; => F
989 iron
990 ;; following fails, need to ensure that we work on list elts, not just
991 ;; elts within a list:
993 ;; (coerce iron 'real)
995 ;; the following is a general list-conversion coercion approach -- is
996 ;; there a more efficient way?
997 ;; (coerce 1 'real)
998 ;; (mapcar #'(lambda (x) (coerce x 'double-float)) iron)
1000 (princ "Data Set up"))
1005 (progn ;; Data setup
1007 (describe 'make-matrix)
1009 (defparameter *indep-vars-2-matrix*
1010 (make-matrix (length iron) 2
1011 :initial-contents
1012 (mapcar #'(lambda (x y)
1013 (list (coerce x 'double-float)
1014 (coerce y 'double-float)))
1015 iron aluminum)))
1018 (defparameter *dep-var*
1019 (make-vector (length absorbtion)
1020 :type :row
1021 :initial-contents
1022 (list
1023 (mapcar #'(lambda (x) (coerce x 'double-float))
1024 absorbtion))))
1026 (make-dataframe *dep-var*)
1027 (make-dataframe (transpose *dep-var*))
1029 (defparameter *dep-var-int*
1030 (make-vector (length absorbtion)
1031 :type :row
1032 :element-type 'integer
1033 :initial-contents (list absorbtion)))
1036 (defparameter *xv+1a*
1037 (make-matrix
1039 :initial-contents #2A((1d0 1d0)
1040 (1d0 3d0)
1041 (1d0 2d0)
1042 (1d0 4d0)
1043 (1d0 3d0)
1044 (1d0 5d0)
1045 (1d0 4d0)
1046 (1d0 6d0))))
1048 (defparameter *xv+1b*
1049 (bind2
1050 (ones 8 1)
1051 (make-matrix
1053 :initial-contents '((1d0)
1054 (3d0)
1055 (2d0)
1056 (4d0)
1057 (3d0)
1058 (5d0)
1059 (4d0)
1060 (6d0)))
1061 :by :column))
1063 (m= *xv+1a* *xv+1b*) ; => T
1065 (princ "Data Set up"))
1069 ;;;; LM
1071 (progn
1073 (defparameter *y*
1074 (make-vector
1076 :type :row
1077 :initial-contents '((1d0 2d0 3d0 4d0 5d0 6d0 7d0 8d0))))
1080 (defparameter *xv+1*
1081 (make-matrix
1083 :initial-contents '((1d0 1d0)
1084 (1d0 3d0)
1085 (1d0 2d0)
1086 (1d0 4d0)
1087 (1d0 3d0)
1088 (1d0 5d0)
1089 (1d0 4d0)
1090 (1d0 6d0))))
1093 ;; so something like (NOTE: matrices are transposed to begin with, hence the incongruety)
1094 (defparameter *xtx-2* (m* (transpose *xv+1*) *xv+1*))
1095 ;; #<LA-SIMPLE-MATRIX-DOUBLE 2 x 2
1096 ;; 8.0d0 28.0d0
1097 ;; 28.0d0 116.0d0>
1099 (defparameter *xty-2* (m* (transpose *xv+1*) (transpose *y*)))
1100 ;; #<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
1101 ;; 36.0d0
1102 ;; 150.0d0>
1104 (defparameter *rcond-2* 0.000001)
1105 (defparameter *betahat-2* (gelsy *xtx-2* *xty-2* *rcond-2*))
1106 ;; *xtx-2* => "details of complete orthogonal factorization"
1107 ;; according to man page:
1108 ;; #<LA-SIMPLE-MATRIX-DOUBLE 2 x 2
1109 ;; -119.33147112141039d0 -29.095426104883202d0
1110 ;; 0.7873402682880205d0 -1.20672274167718d0>
1112 ;; *xty-2* => output becomes solution:
1113 ;; #<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
1114 ;; -0.16666666666668312d0
1115 ;; 1.333333333333337d0>
1117 *betahat-2* ; which matches R, see below
1119 (documentation 'gelsy 'function)
1122 ;; (#<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
1123 ;; -0.16666666666668312 1.333333333333337>
1124 ;; 2)
1126 ;; ## Test case in R:
1127 ;; x <- c( 1.0, 3.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0)
1128 ;; y <- c( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
1129 ;; lm(y~x)
1130 ;; ## => Call: lm(formula = y ~ x)
1132 ;; Coefficients: (Intercept) x
1133 ;; -0.1667 1.3333
1135 ;; summary(lm(y~x))
1136 ;; ## =>
1138 ;; Call:
1139 ;; lm(formula = y ~ x)
1141 ;; Residuals:
1142 ;; Min 1Q Median 3Q Max
1143 ;; -1.833e+00 -6.667e-01 -3.886e-16 6.667e-01 1.833e+00
1145 ;; Coefficients:
1146 ;; Estimate Std. Error t value Pr(>|t|)
1147 ;; (Intercept) -0.1667 1.1587 -0.144 0.89034
1148 ;; x 1.3333 0.3043 4.382 0.00466 **
1149 ;; ---
1150 ;; Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
1152 ;; Residual standard error: 1.291 on 6 degrees of freedom
1153 ;; Multiple R-squared: 0.7619, Adjusted R-squared: 0.7222
1154 ;; F-statistic: 19.2 on 1 and 6 DF, p-value: 0.004659
1158 ;; which suggests one might do (modulo ensuring correct
1159 ;; orientations). When this is finalized, it should migrate to
1160 ;; CLS.
1164 (defparameter *n* 20) ; # rows = # obsns
1165 (defparameter *p* 10) ; # cols = # vars
1166 (defparameter *x-temp* (rand *n* *p*))
1167 (defparameter *b-temp* (rand *p* 1))
1168 (defparameter *y-temp* (m* *x-temp* *b-temp*))
1169 ;; so Y=Xb + \eps
1170 (defparameter *rcond* (* (coerce (expt 2 -52) 'double-float)
1171 (max (nrows *x-temp*) (ncols *y-temp*))))
1172 (defparameter *orig-x* (copy *x-temp*))
1173 (defparameter *orig-b* (copy *b-temp*))
1174 (defparameter *orig-y* (copy *y-temp*))
1176 (defparameter *lm-result* (lm *x-temp* *y-temp*))
1177 (princ (first *lm-result*))
1178 (princ (second *lm-result*))
1179 (princ (third *lm-result*))
1180 (v= (third *lm-result*)
1181 (v- (first (first *lm-result*))
1182 (first (second *lm-result*))))
1187 ;; Some issues exist in the LAPACK vs. LINPACK variants, hence R
1188 ;; uses LINPACK primarily, rather than LAPACK. See comments in R
1189 ;; source for issues.
1192 ;; Goal is to start from X, Y and then realize that if
1193 ;; Y = X \beta, then, i.e. 8x1 = 8xp px1 + 8x1
1194 ;; XtX \hat\beta = Xt Y
1195 ;; so that we can solve the equation W \beta = Z where W and Z
1196 ;; are known, to estimate \beta.
1198 ;; the above is known to be numerically instable -- some processing
1199 ;; of X is preferred and should be done prior. And most of the
1200 ;; transformation-based work does precisely that.
1202 ;; recall: Var[Y] = E[(Y - E[Y])(Y-E[Y])t]
1203 ;; = E[Y Yt] - 2 \mu \mut + \mu \mut
1204 ;; = E[Y Yt] - \mu \mut
1206 ;; Var Y = E[Y^2] - \mu^2
1209 ;; For initial estimates of covariance of \hat\beta:
1211 ;; \hat\beta = (Xt X)^-1 Xt Y
1212 ;; with E[ \hat\beta ]
1213 ;; = E[ (Xt X)^-1 Xt Y ]
1214 ;; = E[(Xt X)^-1 Xt (X\beta)]
1215 ;; = \beta
1217 ;; So Var[\hat\beta] = ...
1218 ;; (Xt X)
1219 ;; and this gives SE(\beta_i) = (* (sqrt (mref Var i i)) adjustment)
1222 ;; from docs:
1224 (setf *temp-result*
1225 (let ((*default-implementation* :foreign-array))
1226 (let* ((m 10)
1227 (n 10)
1228 (a (rand m n))
1229 (x (rand n 1))
1230 (b (m* a x))
1231 (rcond (* (coerce (expt 2 -52) 'double-float)
1232 (max (nrows a) (ncols a))))
1233 (orig-a (copy a))
1234 (orig-b (copy b))
1235 (orig-x (copy x)))
1236 (list x (gelsy a b rcond))
1237 ;; no applicable conversion?
1238 ;; (m- (#<FA-SIMPLE-VECTOR-DOUBLE (10 x 1))
1239 ;; (#<FA-SIMPLE-VECTOR-DOUBLE (10 x 1)) )
1240 (v- x (first (gelsy a b rcond))))))
1243 (princ *temp-result*)
1245 (setf *temp-result*
1246 (let ((*default-implementation* :lisp-array))
1247 (let* ((m 10)
1248 (n 10)
1249 (a (rand m n))
1250 (x (rand n 1))
1251 (b (m* a x))
1252 (rcond (* (coerce (expt 2 -52) 'double-float)
1253 (max (nrows a) (ncols a))))
1254 (orig-a (copy a))
1255 (orig-b (copy b))
1256 (orig-x (copy x)))
1257 (list x (gelsy a b rcond))
1258 (m- x (first (gelsy a b rcond)))
1260 (princ *temp-result*)
1263 (defparameter *xv*
1264 (make-vector
1266 :type :row ;; default, not usually needed!
1267 :initial-contents '((1d0 3d0 2d0 4d0 3d0 5d0 4d0 6d0))))
1269 (defparameter *y*
1270 (make-vector
1272 :type :row
1273 :initial-contents '((1d0 2d0 3d0 4d0 5d0 6d0 7d0 8d0))))
1275 ;; so something like (NOTE: matrices are transposed to begin with, hence the incongruety)
1276 (defparameter *xtx-1* (m* *xv* (transpose *xv*)))
1277 (defparameter *xty-1* (m* *xv* (transpose *y*)))
1278 (defparameter *rcond-in* (* (coerce (expt 2 -52) 'double-float)
1279 (max (nrows *xtx-1*)
1280 (ncols *xty-1*))))
1282 (defparameter *betahat* (gelsy *xtx-1* *xty-1* *rcond-in*))
1284 ;; (#<LA-SIMPLE-VECTOR-DOUBLE (1 x 1)
1285 ;; 1.293103448275862>
1286 ;; 1)
1288 ;; ## Test case in R:
1289 ;; x <- c( 1.0, 3.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0)
1290 ;; y <- c( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
1291 ;; lm(y~x-1)
1292 ;; ## =>
1293 ;; Call:
1294 ;; lm(formula = y ~ x - 1)
1296 ;; Coefficients:
1297 ;; x
1298 ;; 1.293
1300 (first *betahat*))
1305 (type-of #2A((1 2 3 4 5)
1306 (10 20 30 40 50)))
1308 (type-of (rand 10 20))
1310 (typep #2A((1 2 3 4 5)
1311 (10 20 30 40 50))
1312 'matrix-like)
1314 (typep (rand 10 20) 'matrix-like)
1316 (typep #2A((1 2 3 4 5)
1317 (10 20 30 40 50))
1318 'array)
1320 (typep (rand 10 20) 'array)
1323 ;;;;;;;;;;;;;;;;; ===========
1325 (defparameter *my-df-trees*
1326 (make-instance 'dataframe-array
1327 :storage
1328 (listoflist->array
1329 (transpose-listoflist
1330 (rsm.string:file->number-table
1331 (localized-pathto "Data/trees.csv"))))
1332 :doc "This is an interesting dataframe-array that currently fails"))
1334 ;; *my-df-trees*
1336 (defparameter *my-df-trees2*
1337 (make-instance 'dataframe-array
1338 :storage
1339 (listoflist->array
1340 (transpose-listoflist
1341 (rsm.string:file->number-table
1342 (localized-pathto "Data/trees.csv")
1343 :delims ",")))
1344 :doc "This is an interesting dataframe-array that currently fails"))
1345 ;; *my-df-trees2*
1346 ;; (dataset *my-df-trees2*)
1348 (defparameter *my-df-trees2a*
1349 (make-instance 'dataframe-array
1350 :storage
1351 (listoflist->array
1352 (rsm.string:file->number-table
1353 (localized-pathto "Data/trees.csv")
1354 :delims ","))
1355 :doc "This is an interesting dataframe-array that currently fails"))
1357 ;; *my-df-trees2a*
1358 ;; (dataset *my-df-trees2a*)