3 ;;; Time-stamp: <2009-04-20 08:08: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
9 ;;; progns... This file contains the current challenges to
10 ;;; solve, including a description of the setup and the work
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.
20 ;;(asdf:oos 'asdf:load-op 'lisp-matrix)
21 ;;(asdf:oos 'asdf:compile-op 'lispstat :force t)
22 ;;(asdf:oos 'asdf:load-op 'lispstat)
24 (in-package :lisp-stat-unittests
)
26 ;; tests = 80, failures = 8, errors = 15
27 (run-tests :suite
'lisp-stat-ut
)
28 (describe (run-tests :suite
'lisp-stat-ut
))
30 ;; FIXME: Example: currently not relevant, yet
31 ;; (describe (lift::run-test :test-case 'lisp-stat-unittests::create-proto
32 ;; :suite 'lisp-stat-unittests::lisp-stat-ut-proto))
34 (describe (lift::run-tests
:suite
'lisp-stat-ut-dataframe
))
35 (lift::run-tests
:suite
'lisp-stat-ut-dataframe
)
39 :test-case
'lisp-stat-unittests
::create-proto
40 :suite
'lisp-stat-unittests
::lisp-stat-ut-proto
))
42 (describe 'lisp-stat-ut
)
46 (progn ;; FIXME: Regression modeling (some data future-ish).
49 ;; - confirm estimates for multivariate case,
50 ;; - pretty-print output
51 ;; - fix up API -- what do we want this to look like?
54 (regression-model (list->vector-like iron
) ;; BROKEN
55 (list->vector-like absorbtion
))
65 (covariance-matrix *m-fit
*)
68 (regression-model (transpose (listoflist->matrix-like
(list iron aluminum
)
69 :orientation
:row-major
))
70 (list->vector-like absorbtion
) ))
72 (defparameter *m3-fit
*
75 ;; Should the above look something like:
76 (defparameter *m3-fit
*
77 (spec-and-fit-model '(absorbtion = iron aluminum
)))
78 ;; in which case we split the list before/after the "=" character.
82 (covariance-matrix *m3-fit
*))
86 (progn ;; FIXME: Need to clean up data examples, licenses, attributions, etc.
87 ;; The following breaks because we should use a package to hold
88 ;; configuration details, and this would be the only package outside
89 ;; of packages.lisp, as it holds the overall defsystem structure.
90 (load-data "iris.lsp") ;; (the above partially fixed).
96 (progn ;; Importing data from DSV text files.
98 (defparameter *my-df-2
*
99 (make-instance 'dataframe-array
102 (cybertiggyr-dsv::load-escaped
103 "/media/disk/Desktop/sandbox/CLS.git/Data/example-mixed.csv"))
104 :doc
"This is an interesting dataframe-array"))
105 #|
:case-labels
(list "x" "y")
106 :var-labels
(list "a" "b" "c" "d" "e")
109 ;; a better approach is:
110 (asdf:oos
'asdf
:load-op
'rsm-string
)
111 (rsm.string
:file-
>string-table
112 "/media/disk/Desktop/sandbox/CLS.git/Data/example-mixed.csv")
114 (rsm.string
:file-
>number-table
115 "/media/disk/Desktop/sandbox/CLS.git/Data/example-numeric.csv")
117 (defparameter *my-df-2
*
118 (make-instance 'dataframe-array
121 (transpose-listoflist
122 (rsm.string
:file-
>string-table
123 "/media/disk/Desktop/sandbox/CLS.git/Data/example-mixed.csv")))
124 :doc
"This is an interesting dataframe-array"))
127 (defparameter *my-df-3
*
128 (make-instance 'dataframe-array
131 (transpose-listoflist
132 (rsm.string
:file-
>number-table
133 "/media/disk/Desktop/sandbox/CLS.git/Data/example-numeric.csv")))
134 :doc
"This is an interesting dataframe-array"))
144 (asdf:oos
'asdf
:load-op
'cl-plplot
)
147 (defparameter *gdev
* "xcairo")
148 ;; (defparameter *gdev* "xwin")
149 (cl-plplot::plsdev
*gdev
*)
151 ;; there is currently a loose pointer floating around that causes
152 ;; errors the 3rd time that we create a plot (and crashes SBCL the
153 ;; 4th time). Order independent.
162 (in-package :ls-user
)
166 ;; REVIEW: general Lisp use guidance
168 (fdefinition 'make-matrix
)
169 (documentation 'make-matrix
'function
)
171 #| Examples from CLHS
, a bit of guidance.
173 ;; This function assumes its callers have checked the types of the
174 ;; arguments, and authorizes the compiler to build in that assumption.
175 (defun discriminant (a b c
)
176 (declare (number a b c
))
177 "Compute the discriminant for a quadratic equation."
178 (- (* b b
) (* 4 a c
))) => DISCRIMINANT
179 (discriminant 1 2/3 -
2) => 76/9
181 ;; This function assumes its callers have not checked the types of the
182 ;; arguments, and performs explicit type checks before making any assumptions.
183 (defun careful-discriminant (a b c
)
184 "Compute the discriminant for a quadratic equation."
185 (check-type a number
)
186 (check-type b number
)
187 (check-type c number
)
188 (locally (declare (number a b c
))
189 (- (* b b
) (* 4 a c
)))) => CAREFUL-DISCRIMINANT
190 (careful-discriminant 1 2/3 -
2) => 76/9
196 (progn ;; experiments with GSL and the Lisp interface.
197 (asdf:oos
'asdf
:load-op
'gsll
)
198 (asdf:oos
'asdf
:load-op
'gsll-tests
)
200 ;; the following should be equivalent
201 (setf *t1
* (LIST 6.18d0
6.647777777777779d0
6.18d0
))
202 (setf *t2
* (MULTIPLE-VALUE-LIST
204 (gsll:make-marray
'DOUBLE-FLOAT
205 :INITIAL-CONTENTS
'(-3.21d0
1.0d0
12.8d0
)))
207 (gsll:MAKE-MARRAY
'DOUBLE-FLOAT
208 :INITIAL-CONTENTS
'(3.0d0
1.0d0
2.0d0
))))
209 (LET ((MEAN (gsll:MEAN VEC
)))
210 (LIST (gsll:ABSOLUTE-DEVIATION VEC
)
211 (gsll:WEIGHTED-ABSOLUTE-DEVIATION VEC WEIGHTS
)
212 (gsll:ABSOLUTE-DEVIATION VEC MEAN
))))))
215 ;; from (gsll:examples 'gsll::numerical-integration) ...
216 (gsll:integration-qng gsll
::one-sine
0.0d0 PI
)
218 (gsll:defun-single axpb
(x) (+ (* 2 x
) 3)) ;; a<-2, b<-3
219 (gsll:integration-qng axpb
1d0
2d0
)
223 (defun-single axpb2
(x) (+ (* a x
) b
)))
224 (gsll:integration-qng axpb2
1d0
2d0
)
227 ;; (gsll:integration-qng
230 ;; (defun-single axpb2 (x) (+ (* a x) b)))
233 ;; right, but weird expansion...
234 (gsll:integration-qng
237 (defun axpb2 (x) (+ (* a x
) b
))
238 (gsll:def-single-function axpb2
)
242 ;; Linear least squares
244 (gsll:gsl-lookup
"gsl_linalg_LU_decomp") ; => gsll:lu-decomposition
245 (gsll:gsl-lookup
"gsl_linalg_LU_solve") ; => gsll:lu-solve
251 (progn ;; philosophy time
253 (setf my-model
(model :name
"ex1"
254 :data-slots
(list w x y z
)
255 :param-slots
(list alpha beta gamma
)
256 :math-form
(regression-model :formula
'(= w
(+ (* beta x
)
260 :centrality
'median
; 'mean
267 (setf my-dataset
(statistical-table :table data-frame-contents
268 :metadata
(list (:case-names
(list ))
270 (:documentation
"string of doc"))))
272 (setf my-analysis
(analysis
275 :parameter-map
(pairing (model-param-slots my-model
)
276 (data-var-names my-dataset
))))
278 ;; ontological implications -- the analysis is an abstract class of
279 ;; data, model, and mapping between the model and data. The fit is
280 ;; the instantiation of such. This provides a statistical object
281 ;; computation theory which can be realized as "executable
282 ;; statistics" or "computable statistics".
283 (setf my-analysis
(analyze my-fit
284 :estimation-method
'linear-least-squares-regression
))
286 ;; one of the tricks here is that one needs to provide the structure
287 ;; from which to consider estimation, and more importantly, the
288 ;; validity of the estimation.
291 (setf linear-least-squares-regression
292 (estimation-method-definition
293 :variable-defintions
((list
294 ;; from MachLearn: supervised,
296 :data-response-vars list-drv
; nil if unsup
299 :data-predictor-vars list-dpv
300 ;; nil in this case. these
301 ;; describe "out-of-box" specs
302 :hyper-vars list-hv
))
303 :form
'(regression-additive-error
304 :central-form
(linear-form drv pv dpv
)
305 :error-form
'normal-error
)
306 :resulting-decision
'(point-estimation interval-estimation
)
307 :philosophy
'frequentist
308 :documentation
"use least squares to fit a linear regression
311 (defparameter *statistical-philosophies
*
312 '(frequentist bayesian fiducial decision-analysis
)
313 "can be combined to build decision-making approaches and
316 (defparameter *decisions
*
317 '(estimation selection testing
)
318 "possible results from a...")
319 ;; is this really true? One can embedded hypothesis testing within
320 ;; estimation, as the hypothesis estimated to select. And
321 ;; categorical/continuous rear their ugly heads, but not really in
324 (defparameter *ontology-of-decision-procedures
*
328 (list :maximum-likelihood
333 (list :maximum-likelihood
339 :bioequivalence-inversion
)
344 :partially-parametric
))
345 "start of ontology"))
356 :initial-contents
'((1d0 2d0
3d0
4d0
5d0
6d0
7d0
8d0
))))
362 :initial-contents
'((1d0 1d0
)
372 ;; so something like (NOTE: matrices are transposed to begin with, hence the incongruety)
373 (defparameter *xtx-2
* (m* (transpose *xv
+1*) *xv
+1*))
374 ;; #<LA-SIMPLE-MATRIX-DOUBLE 2 x 2
378 (defparameter *xty-2
* (m* (transpose *xv
+1*) (transpose *y
*)))
379 ;; #<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
383 (defparameter *rcond-2
* 0.000001)
384 (defparameter *betahat-2
* (gelsy *xtx-2
* *xty-2
* *rcond-2
*))
385 ;; *xtx-2* => "details of complete orthogonal factorization"
386 ;; according to man page:
387 ;; #<LA-SIMPLE-MATRIX-DOUBLE 2 x 2
388 ;; -119.33147112141039d0 -29.095426104883202d0
389 ;; 0.7873402682880205d0 -1.20672274167718d0>
391 ;; *xty-2* => output becomes solution:
392 ;; #<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
393 ;; -0.16666666666668312d0
394 ;; 1.333333333333337d0>
396 *betahat-2
* ; which matches R, see below
398 (documentation 'gelsy
'function
)
401 ;; (#<LA-SIMPLE-VECTOR-DOUBLE (2 x 1)
402 ;; -0.16666666666668312 1.333333333333337>
405 ;; ## Test case in R:
406 ;; x <- c( 1.0, 3.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0)
407 ;; y <- c( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
409 ;; ## => Call: lm(formula = y ~ x)
411 ;; Coefficients: (Intercept) x
418 ;; lm(formula = y ~ x)
421 ;; Min 1Q Median 3Q Max
422 ;; -1.833e+00 -6.667e-01 -3.886e-16 6.667e-01 1.833e+00
425 ;; Estimate Std. Error t value Pr(>|t|)
426 ;; (Intercept) -0.1667 1.1587 -0.144 0.89034
427 ;; x 1.3333 0.3043 4.382 0.00466 **
429 ;; Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
431 ;; Residual standard error: 1.291 on 6 degrees of freedom
432 ;; Multiple R-squared: 0.7619, Adjusted R-squared: 0.7222
433 ;; F-statistic: 19.2 on 1 and 6 DF, p-value: 0.004659
437 ;; which suggests one might do (modulo ensuring correct
438 ;; orientations). When this is finalized, it should migrate to
443 (defparameter *n
* 20) ; # rows = # obsns
444 (defparameter *p
* 10) ; # cols = # vars
445 (defparameter *x-temp
* (rand *n
* *p
*))
446 (defparameter *b-temp
* (rand *p
* 1))
447 (defparameter *y-temp
* (m* *x-temp
* *b-temp
*))
449 (defparameter *rcond
* (* (coerce (expt 2 -
52) 'double-float
)
450 (max (nrows *x-temp
*) (ncols *y-temp
*))))
451 (defparameter *orig-x
* (copy *x-temp
*))
452 (defparameter *orig-b
* (copy *b-temp
*))
453 (defparameter *orig-y
* (copy *y-temp
*))
455 (defparameter *lm-result
* (lm *x-temp
* *y-temp
*))
456 (princ (first *lm-result
*))
457 (princ (second *lm-result
*))
458 (princ (third *lm-result
*))
459 (v= (third *lm-result
*)
460 (v- (first (first *lm-result
*))
461 (first (second *lm-result
*))))
466 ;; Some issues exist in the LAPACK vs. LINPACK variants, hence R
467 ;; uses LINPACK primarily, rather than LAPACK. See comments in R
468 ;; source for issues.
471 ;; Goal is to start from X, Y and then realize that if
472 ;; Y = X \beta, then, i.e. 8x1 = 8xp px1 + 8x1
473 ;; XtX \hat\beta = Xt Y
474 ;; so that we can solve the equation W \beta = Z where W and Z
475 ;; are known, to estimate \beta.
477 ;; the above is known to be numerically instable -- some processing
478 ;; of X is preferred and should be done prior. And most of the
479 ;; transformation-based work does precisely that.
481 ;; recall: Var[Y] = E[(Y - E[Y])(Y-E[Y])t]
482 ;; = E[Y Yt] - 2 \mu \mut + \mu \mut
483 ;; = E[Y Yt] - \mu \mut
485 ;; Var Y = E[Y^2] - \mu^2
488 ;; For initial estimates of covariance of \hat\beta:
490 ;; \hat\beta = (Xt X)^-1 Xt Y
491 ;; with E[ \hat\beta ]
492 ;; = E[ (Xt X)^-1 Xt Y ]
493 ;; = E[(Xt X)^-1 Xt (X\beta)]
496 ;; So Var[\hat\beta] = ...
498 ;; and this gives SE(\beta_i) = (* (sqrt (mref Var i i)) adjustment)
504 (let ((*default-implementation
* :foreign-array
))
510 (rcond (* (coerce (expt 2 -
52) 'double-float
)
511 (max (nrows a
) (ncols a
))))
515 (list x
(gelsy a b rcond
))
516 ;; no applicable conversion?
517 ;; (m- (#<FA-SIMPLE-VECTOR-DOUBLE (10 x 1))
518 ;; (#<FA-SIMPLE-VECTOR-DOUBLE (10 x 1)) )
519 (v- x
(first (gelsy a b rcond
))))))
522 (princ *temp-result
*)
525 (let ((*default-implementation
* :lisp-array
))
531 (rcond (* (coerce (expt 2 -
52) 'double-float
)
532 (max (nrows a
) (ncols a
))))
536 (list x
(gelsy a b rcond
))
537 (m- x
(first (gelsy a b rcond
)))
539 (princ *temp-result
*)
545 :type
:row
;; default, not usually needed!
546 :initial-contents
'((1d0 3d0
2d0
4d0
3d0
5d0
4d0
6d0
))))
552 :initial-contents
'((1d0 2d0
3d0
4d0
5d0
6d0
7d0
8d0
))))
554 ;; so something like (NOTE: matrices are transposed to begin with, hence the incongruety)
555 (defparameter *xtx-1
* (m* *xv
* (transpose *xv
*)))
556 (defparameter *xty-1
* (m* *xv
* (transpose *y
*)))
557 (defparameter *rcond-in
* (* (coerce (expt 2 -
52) 'double-float
)
561 (defparameter *betahat
* (gelsy *xtx-1
* *xty-1
* *rcond-in
*))
563 ;; (#<LA-SIMPLE-VECTOR-DOUBLE (1 x 1)
564 ;; 1.293103448275862>
567 ;; ## Test case in R:
568 ;; x <- c( 1.0, 3.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0)
569 ;; y <- c( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
573 ;; lm(formula = y ~ x - 1)
584 (type-of #2A
((1 2 3 4 5)
587 (type-of (rand 10 20))
589 (typep #2A
((1 2 3 4 5)
593 (typep (rand 10 20) 'matrix-like
)
595 (typep #2A
((1 2 3 4 5)
599 (typep (rand 10 20) 'array
)