From f7d1922e7963d9862a53552b7a98a47077890889 Mon Sep 17 00:00:00 2001 From: AJ Rossini Date: Tue, 10 Mar 2009 08:50:03 +0100 Subject: [PATCH] implement and export xtxinv. include within regression calculations. Signed-off-by: AJ Rossini --- src/packages.lisp | 5 ++--- src/stat-models/regression.lisp | 24 +++++++++++++++++++----- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/src/packages.lisp b/src/packages.lisp index 6e11fc9..3598d27 100644 --- a/src/packages.lisp +++ b/src/packages.lisp @@ -1,6 +1,6 @@ ;;; -*- mode: lisp -*- -;;; Time-stamp: <2009-02-17 08:22:21 tony> +;;; Time-stamp: <2009-03-10 08:43:49 tony> ;;; Creation: <2008-03-11 19:18:34 user> ;;; File: packages.lisp ;;; Author: AJ Rossini @@ -310,8 +310,7 @@ total-sum-of-squares residual-sum-of-squares predictor-names response-name case-labels ;; functions for helpers - lm xtxinv - )) + lm xtxinv )) (defpackage :lisp-stat (:documentation "Experimentation package for LispStat. Serious diff --git a/src/stat-models/regression.lisp b/src/stat-models/regression.lisp index 12dbcee..570051d 100644 --- a/src/stat-models/regression.lisp +++ b/src/stat-models/regression.lisp @@ -52,6 +52,23 @@ + +(defun xtxinv (x) + "In: X + Out: (XtX)^-1 + +X is NxP, so result is PxP. Represents Var[\hat\beta], the vars for +\hat \beta from Y = X \beta + \eps. Done by Cholesky decomposition, +using LAPACK's dpotri routine to invert, after factorizing with dpotrf. + + + (let ((m1 (rand 7 5))) + (xtxinv m1)) +" + (check-type x matrix-like) + (minv-cholesky (m* (transpose x) x))) + + ;; might add args: (method 'gelsy), or do we want to put a more ;; general front end, linear-least-square, across the range of ;; LAPACK solvers? @@ -538,12 +555,9 @@ basis." (coerce (compound-data-seq (select m (1+ n) indices)) 'list))) ;; ERROR (defmeth regression-model-proto :xtxinv () -"Message args: () + "Message args: () Returns ((X^T) X)^(-1) or ((X^T) W X)^(-1)." - (let ((indices (if (send self :intercept) - (cons 0 (1+ (send self :basis))) - (1+ (send self :basis))))) - (select (send self :sweep-matrix) indices indices))) + (xtxinv (send self x))) (defmeth regression-model-proto :coef-standard-errors () "Message args: () -- 2.11.4.GIT