3 ;;; Time-stamp: <2009-03-30 08:17:46 tony>
4 ;;; Creation: <2008-09-08 08:06:30 tony>
6 ;;; Author: AJ Rossini <blindglobe@gmail.com>
7 ;;; Copyright: (c) 2007-2008, AJ Rossini <blindglobe@gmail.com>. BSD.
8 ;;; Purpose: Stuff that needs to be made working sits inside the progns...
10 ;;; What is this talk of 'release'? Klingons do not make software
11 ;;; 'releases'. Our software 'escapes', leaving a bloody trail of
12 ;;; designers and quality assurance people in its wake.
14 ;;; This file contains the current challenges to solve, including a
15 ;;; description of the setup and the work to solve....
20 ;;(asdf:oos 'asdf:compile-op 'lispstat)
21 ;;(asdf:oos 'asdf:load-op 'lispstat)
23 (in-package :lisp-stat-unittests
)
25 ;; tests = 87, failures = 9, errors = 22
27 (describe (run-tests :suite
'lisp-stat-ut
))
28 (run-tests :suite
'lisp-stat-ut
)
31 ;; FIXME: Example: currently not relevant, yet
34 :test-case
'lisp-stat-unittests
::create-proto
35 :suite
'lisp-stat-unittests
::lisp-stat-ut-proto
))
38 (describe 'lisp-stat-ut
)
44 (describe (lift::run-tests
:suite
'lisp-stat-ut-dataclos
))
45 (lift::run-tests
:suite
'lisp-stat-ut-dataclos
)
49 :test-case
'lisp-stat-unittests
::create-proto
50 :suite
'lisp-stat-unittests
::lisp-stat-ut-proto
))
52 (defparameter *my-df-1
*
53 (make-instance 'dataframe-array
54 :storage
#2A
((1 2 3 4 5)
56 :doc
"This is an interesting dataframe-array"
57 :case-labels
(list "x" "y")
58 :var-labels
(list "a" "b" "c" "d" "e")))
60 (make-dataframe #2A
((1 2 3 4 5)
63 (make-dataframe (rand 4 3))
65 (defparameter *my-df-1
*
66 (make-dataframe #2A
((1 2 3 4 5)
68 :caselabels
(list "x" "y")
69 :varlabels
(list "a" "b" "c" "d" "e")
70 :doc
"This is an interesting dataframe-array"))
72 (caselabels *my-df-1
*)
76 (defparameter *my-df-2
*
77 (make-instance 'dataframe-array
79 (make-array-from-listoflists
80 (cybertiggyr-dsv::load-escaped
81 "/media/disk/Desktop/sandbox/CLS.git/Data/example-mixed.csv"))
82 :doc
"This is an interesting dataframe-array"))
83 #|
:case-labels
(list "x" "y")
84 :var-labels
(list "a" "b" "c" "d" "e")
92 (describe 'make-matrix
)
94 (defparameter *indep-vars-2-matrix
*
95 (make-matrix (length iron
) 2
97 (mapcar #'(lambda (x y
)
98 (list (coerce x
'double-float
)
99 (coerce y
'double-float
)))
103 (defparameter *dep-var
*
104 (make-vector (length absorbtion
)
108 (mapcar #'(lambda (x) (coerce x
'double-float
))
111 (make-dataframe *dep-var
*)
112 (make-dataframe (transpose *dep-var
*))
114 (defparameter *dep-var-int
*
115 (make-vector (length absorbtion
)
117 :element-type
'integer
118 :initial-contents
(list absorbtion
)))
121 (defparameter *xv
+1a
*
124 :initial-contents
#2A
((1d0 1d0
)
133 (defparameter *xv
+1b
*
138 :initial-contents
'((1d0)
148 (m= *xv
+1a
* *xv
+1b
*) ; => T
150 (princ "Data Set up"))
156 ;; REVIEW: general Lisp use guidance
158 (fdefinition 'make-matrix
)
159 (documentation 'make-matrix
'function
)
161 #| Examples from CLHS
, a bit of guidance.
163 ;; This function assumes its callers have checked the types of the
164 ;; arguments, and authorizes the compiler to build in that assumption.
165 (defun discriminant (a b c
)
166 (declare (number a b c
))
167 "Compute the discriminant for a quadratic equation."
168 (- (* b b
) (* 4 a c
))) => DISCRIMINANT
169 (discriminant 1 2/3 -
2) => 76/9
171 ;; This function assumes its callers have not checked the types of the
172 ;; arguments, and performs explicit type checks before making any assumptions.
173 (defun careful-discriminant (a b c
)
174 "Compute the discriminant for a quadratic equation."
175 (check-type a number
)
176 (check-type b number
)
177 (check-type c number
)
178 (locally (declare (number a b c
))
179 (- (* b b
) (* 4 a c
)))) => CAREFUL-DISCRIMINANT
180 (careful-discriminant 1 2/3 -
2) => 76/9
186 (progn ;; FIXME: Regression modeling
188 ;; data setup in previous FIXME
189 (defparameter *m
* nil
191 ;; need to make vectors and matrices from the lists...
194 (def *m
* (regression-model (list->vector-like iron
)
195 (list->vector-like absorbtion
)))
197 (def m
(regression-model (list->vector-like iron
)
198 (list->vector-like absorbtion
) :print nil
))
202 (send m
:own-methods
)
203 ;; (lsos::ls-objects-methods m) ; bogus?
206 (def m
(regression-model (list->vector-like iron
)
207 (list->vector-like absorbtion
)))
209 (def m
(regression-model (listoflists->matrix-like
(list iron aluminum
))
210 (list->vector-like absorbtion
) :print nil
))
214 (send m
:sweep-matrix
)
215 (format t
"~%~A~%" (send m
:sweep-matrix
))
217 ;; need to get multiple-linear regression working (simple linear regr
218 ;; works)... to do this, we need to redo the whole numeric structure,
219 ;; I'm keeping these in as example of brokenness...
221 (send m
:basis
) ;; this should be positive?
222 (send m
:coef-estimates
) )
225 (progn ;; FIXME: Need to clean up data examples, licenses, attributions, etc.
226 ;; The following breaks because we should use a package to hold
227 ;; configuration details, and this would be the only package outside
228 ;; of packages.lisp, as it holds the overall defsystem structure.
229 (load-data "iris.lsp") ;; (the above partially fixed).
236 (progn ;; FIXME: read data from CSV file. To do.
239 ;; challenge is to ensure that we get mixed arrays when we want them,
240 ;; and single-type (simple) arrays in other cases.
243 (defparameter *csv-num
*
244 (cybertiggyr-dsv::load-escaped
245 #p
"/media/disk/Desktop/sandbox/CLS.git/Data/example-numeric.csv"
249 (nth 0 (nth 0 *csv-num
*))
251 (defparameter *csv-num
*
252 (cybertiggyr-dsv::load-escaped
253 #p
"/media/disk/Desktop/sandbox/CLS.git/Data/example-numeric2.dsv"
254 :field-separator
#\
:))
256 (nth 0 (nth 0 *csv-num
*))
259 ;; The handling of these types should be compariable to what we do for
260 ;; matrices, but without the numerical processing. i.e. mref, bind2,
261 ;; make-dataframe, and the class structure should be similar.
263 ;; With numerical data, there should be a straightforward mapping from
264 ;; the data.frame to a matrix. With categorical data (including
265 ;; dense categories such as doc-strings, as well as sparse categories
266 ;; such as binary data), we need to include metadata about ordering,
267 ;; coding, and such. So the structures should probably consider
269 ;; Using the CSV file:
271 (defun parse-number (s)
272 (let* ((*read-eval
* nil
)
273 (n (read-from-string s
)))
279 (parse-number " 34 ")
281 (+ (parse-number "3.4") 3)
282 (parse-number "3.4 ")
283 (parse-number " 3.4")
284 (+ (parse-number " 3.4 ") 3)
288 ;; (coerce "2.3" 'number) => ERROR
289 ;; (coerce "2" 'float) => ERROR
291 (defparameter *csv-num
*
292 (cybertiggyr-dsv::load-escaped
293 #p
"/media/disk/Desktop/sandbox/CLS.git/Data/example-numeric.csv"
295 :filter
#'parse-number
298 (nth 0 (nth 0 *csv-num
*))
300 (defparameter *csv-num
*
301 (cybertiggyr-dsv::load-escaped
302 #p
"/media/disk/Desktop/sandbox/CLS.git/Data/example-numeric2.dsv"
304 :filter
#'parse-number
))
306 (nth 0 (nth 0 *csv-num
*))
308 ;; now we've got the DSV code in the codebase, auto-loaded I hope:
309 cybertiggyr-dsv
:*field-separator
*
310 (defparameter *example-numeric.csv
*
311 (cybertiggyr-dsv:load-escaped
"Data/example-numeric.csv"
312 :field-separator
#\
,))
313 *example-numeric.csv
*
315 ;; the following fails because we've got a bit of string conversion
316 ;; to do. 2 thoughts: #1 modify dsv package, but mucking with
317 ;; encapsulation. #2 add a coercion tool (better, but potentially
319 #+nil
(coerce (nth 3 (nth 3 *example-numeric.csv
*)) 'double-float
)
321 ;; cases, simple to not so
322 (defparameter *test-string1
* "1.2")
323 (defparameter *test-string2
* " 1.2")
324 (defparameter *test-string3
* " 1.2 ")
329 (progn ;; experiments with GSL and the Lisp interface.
330 (asdf:oos
'asdf
:load-op
'gsll
)
331 (asdf:oos
'asdf
:load-op
'gsll-tests
)
333 ;; the following should be equivalent
334 (setf *t1
* (LIST 6.18d0
6.647777777777779d0
6.18d0
))
335 (setf *t2
* (MULTIPLE-VALUE-LIST
337 (gsll:make-marray
'DOUBLE-FLOAT
338 :INITIAL-CONTENTS
'(-3.21d0
1.0d0
12.8d0
)))
340 (gsll:MAKE-MARRAY
'DOUBLE-FLOAT
341 :INITIAL-CONTENTS
'(3.0d0
1.0d0
2.0d0
))))
342 (LET ((MEAN (gsll:MEAN VEC
)))
343 (LIST (gsll:ABSOLUTE-DEVIATION VEC
)
344 (gsll:WEIGHTED-ABSOLUTE-DEVIATION VEC WEIGHTS
)
345 (gsll:ABSOLUTE-DEVIATION VEC MEAN
))))))
348 ;; from (gsll:examples 'gsll::numerical-integration) ...
349 (gsll:integration-qng gsll
::one-sine
0.0d0 PI
)
351 (gsll:defun-single axpb
(x) (+ (* 2 x
) 3)) ;; a<-2, b<-3
352 (gsll:integration-qng axpb
1d0
2d0
)
356 (defun-single axpb2
(x) (+ (* a x
) b
)))
357 (gsll:integration-qng axpb2
1d0
2d0
)
360 ;; (gsll:integration-qng
363 ;; (defun-single axpb2 (x) (+ (* a x) b)))
366 ;; right, but weird expansion...
367 (gsll:integration-qng
370 (defun axpb2 (x) (+ (* a x
) b
))
371 (gsll:def-single-function axpb2
)
375 ;; Linear least squares
377 (gsll:gsl-lookup
"gsl_linalg_LU_decomp") ; => gsll:lu-decomposition
378 (gsll:gsl-lookup
"gsl_linalg_LU_solve") ; => gsll:lu-solve
384 (progn ;; philosophy time
386 (setf my-model
(model :name
"ex1"
387 :data-slots
(list w x y z
)
388 :param-slots
(list alpha beta gamma
)
389 :math-form
(regression-model :formula
'(= w
(+ (* beta x
)
393 :centrality
'median
; 'mean
400 (setf my-dataset
(statistical-table :table data-frame-contents
401 :metadata
(list (:case-names
(list ))
403 (:documentation
"string of doc"))))
405 (setf my-analysis
(analysis
408 :parameter-map
(pairing (model-param-slots my-model
)
409 (data-var-names my-dataset
))))
411 ;; ontological implications -- the analysis is an abstract class of
412 ;; data, model, and mapping between the model and data. The fit is
413 ;; the instantiation of such. This provides a statistical object
414 ;; computation theory which can be realized as "executable
415 ;; statistics" or "computable statistics".
416 (setf my-analysis
(analyze my-fit
417 :estimation-method
'linear-least-squares-regression
))
419 ;; one of the tricks here is that one needs to provide the structure
420 ;; from which to consider estimation, and more importantly, the
421 ;; validity of the estimation.
424 (setf linear-least-squares-regression
425 (estimation-method-definition
426 :variable-defintions
((list
427 ;; from MachLearn: supervised,
429 :data-response-vars list-drv
; nil if unsup
432 :data-predictor-vars list-dpv
433 ;; nil in this case. these
434 ;; describe "out-of-box" specs
435 :hyper-vars list-hv
))
436 :form
'(regression-additive-error
437 :central-form
(linear-form drv pv dpv
)
438 :error-form
'normal-error
)
439 :resulting-decision
'(point-estimation interval-estimation
)
440 :philosophy
'frequentist
441 :documentation
"use least squares to fit a linear regression
444 (defparameter *statistical-philosophies
*
445 '(frequentist bayesian fiducial decision-analysis
)
446 "can be combined to build decision-making approaches and
449 (defparameter *decisions
*
450 '(estimation selection testing
)
451 "possible results from a...")
452 ;; is this really true? One can embedded hypothesis testing within
453 ;; estimation, as the hypothesis estimated to select. And
454 ;; categorical/continuous rear their ugly heads, but not really in
457 (defparameter *ontology-of-decision-procedures
*
461 (list :maximum-likelihood
466 (list :maximum-likelihood
472 :bioequivalence-inversion
)
477 :partially-parametric
))
478 "start of ontology"))
489 :initial-contents
'((1d0 2d0
3d0
4d0
5d0
6d0
7d0
8d0
))))
495 :initial-contents
'((1d0 1d0
)
505 ;; so something like (NOTE: matrices are transposed to begin with, hence the incongruety)
506 (defparameter *xtx-2
* (m* (transpose *xv
+1*) *xv
+1*))
507 ;; #<LA-SIMPLE-MATRIX-DOUBLE 2 x 2
511 (defparameter *xty-2
* (m* (transpose *xv
+1*) (transpose *y
*)))
512 ;; #<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
516 (defparameter *rcond-2
* 0.000001)
517 (defparameter *betahat-2
* (gelsy *xtx-2
* *xty-2
* *rcond-2
*))
518 ;; *xtx-2* => "details of complete orthogonal factorization"
519 ;; according to man page:
520 ;; #<LA-SIMPLE-MATRIX-DOUBLE 2 x 2
521 ;; -119.33147112141039d0 -29.095426104883202d0
522 ;; 0.7873402682880205d0 -1.20672274167718d0>
524 ;; *xty-2* => output becomes solution:
525 ;; #<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
526 ;; -0.16666666666668312d0
527 ;; 1.333333333333337d0>
529 *betahat-2
* ; which matches R, see below
531 (documentation 'gelsy
'function
)
534 ;; (#<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
535 ;; -0.16666666666668312 1.333333333333337>
538 ;; ## Test case in R:
539 ;; x <- c( 1.0, 3.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0)
540 ;; y <- c( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
542 ;; ## => Call: lm(formula = y ~ x)
544 ;; Coefficients: (Intercept) x
551 ;; lm(formula = y ~ x)
554 ;; Min 1Q Median 3Q Max
555 ;; -1.833e+00 -6.667e-01 -3.886e-16 6.667e-01 1.833e+00
558 ;; Estimate Std. Error t value Pr(>|t|)
559 ;; (Intercept) -0.1667 1.1587 -0.144 0.89034
560 ;; x 1.3333 0.3043 4.382 0.00466 **
562 ;; Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
564 ;; Residual standard error: 1.291 on 6 degrees of freedom
565 ;; Multiple R-squared: 0.7619, Adjusted R-squared: 0.7222
566 ;; F-statistic: 19.2 on 1 and 6 DF, p-value: 0.004659
570 ;; which suggests one might do (modulo ensuring correct
571 ;; orientations). When this is finalized, it should migrate to
576 (defparameter *n
* 20) ; # rows = # obsns
577 (defparameter *p
* 10) ; # cols = # vars
578 (defparameter *x-temp
* (rand *n
* *p
*))
579 (defparameter *b-temp
* (rand *p
* 1))
580 (defparameter *y-temp
* (m* *x-temp
* *b-temp
*))
582 (defparameter *rcond
* (* (coerce (expt 2 -
52) 'double-float
)
583 (max (nrows *x-temp
*) (ncols *y-temp
*))))
584 (defparameter *orig-x
* (copy *x-temp
*))
585 (defparameter *orig-b
* (copy *b-temp
*))
586 (defparameter *orig-y
* (copy *y-temp
*))
588 (defparameter *lm-result
* (lm *x-temp
* *y-temp
*))
589 (princ (first *lm-result
*))
590 (princ (second *lm-result
*))
591 (princ (third *lm-result
*))
592 (v= (third *lm-result
*)
593 (v- (first (first *lm-result
*))
594 (first (second *lm-result
*))))
599 ;; Some issues exist in the LAPACK vs. LINPACK variants, hence R
600 ;; uses LINPACK primarily, rather than LAPACK. See comments in R
601 ;; source for issues.
604 ;; Goal is to start from X, Y and then realize that if
605 ;; Y = X \beta, then, i.e. 8x1 = 8xp px1 + 8x1
606 ;; XtX \hat\beta = Xt Y
607 ;; so that we can solve the equation W \beta = Z where W and Z
608 ;; are known, to estimate \beta.
610 ;; the above is known to be numerically instable -- some processing
611 ;; of X is preferred and should be done prior. And most of the
612 ;; transformation-based work does precisely that.
614 ;; recall: Var[Y] = E[(Y - E[Y])(Y-E[Y])t]
615 ;; = E[Y Yt] - 2 \mu \mut + \mu \mut
616 ;; = E[Y Yt] - \mu \mut
618 ;; Var Y = E[Y^2] - \mu^2
621 ;; For initial estimates of covariance of \hat\beta:
623 ;; \hat\beta = (Xt X)^-1 Xt Y
624 ;; with E[ \hat\beta ]
625 ;; = E[ (Xt X)^-1 Xt Y ]
626 ;; = E[(Xt X)^-1 Xt (X\beta)]
629 ;; So Var[\hat\beta] = ...
631 ;; and this gives SE(\beta_i) = (* (sqrt (mref Var i i)) adjustment)
637 (let ((*default-implementation
* :foreign-array
))
643 (rcond (* (coerce (expt 2 -
52) 'double-float
)
644 (max (nrows a
) (ncols a
))))
648 (list x
(gelsy a b rcond
))
649 ;; no applicable conversion?
650 ;; (m- (#<FA-SIMPLE-VECTOR-DOUBLE (10 x 1))
651 ;; (#<FA-SIMPLE-VECTOR-DOUBLE (10 x 1)) )
652 (v- x
(first (gelsy a b rcond
))))))
655 (princ *temp-result
*)
658 (let ((*default-implementation
* :lisp-array
))
664 (rcond (* (coerce (expt 2 -
52) 'double-float
)
665 (max (nrows a
) (ncols a
))))
669 (list x
(gelsy a b rcond
))
670 (m- x
(first (gelsy a b rcond
)))
672 (princ *temp-result
*)
678 :type
:row
;; default, not usually needed!
679 :initial-contents
'((1d0 3d0
2d0
4d0
3d0
5d0
4d0
6d0
))))
685 :initial-contents
'((1d0 2d0
3d0
4d0
5d0
6d0
7d0
8d0
))))
687 ;; so something like (NOTE: matrices are transposed to begin with, hence the incongruety)
688 (defparameter *xtx-1
* (m* *xv
* (transpose *xv
*)))
689 (defparameter *xty-1
* (m* *xv
* (transpose *y
*)))
690 (defparameter *rcond-in
* (* (coerce (expt 2 -
52) 'double-float
)
694 (defparameter *betahat
* (gelsy *xtx-1
* *xty-1
* *rcond-in
*))
696 ;; (#<LA-SIMPLE-VECTOR-DOUBLE (1 x 1)
697 ;; 1.293103448275862>
700 ;; ## Test case in R:
701 ;; x <- c( 1.0, 3.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0)
702 ;; y <- c( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
706 ;; lm(formula = y ~ x - 1)
719 (asdf:oos
'asdf
:load-op
'cl-plplot
)
725 (type-of #2A
((1 2 3 4 5)
728 (type-of (rand 10 20))
730 (typep #2A
((1 2 3 4 5)
734 (typep (rand 10 20) 'matrix-like
)
736 (typep #2A
((1 2 3 4 5)
740 (typep (rand 10 20) 'array
)