;;; Bootstrap lack of fit tests, modified 15 May 2000
;;; Fixed bug so this now works with missing values, 3 Oct 2000.

;;; This duplicates code in updates.lsp and overlay4.lsp
(defmeth polynomial-proto :get-menu-items ()
  (let* ((control self)
         (graph (send self :graph))
         (d (send graph :point-coordinate '(0 1) 
                     (repeat (list (iseq (send graph :num-points))) 2)))
         (extract (send menu-item-proto :new "Extract mean..."
                 :action #'(lambda () (send self :extract))))
         (fbg (send menu-item-proto :new "Fit by marks--general"
                 :action #'(lambda () (send self :fit-by-marks 'general))))
         (fbgp (send menu-item-proto :new "Fit by marks--parallel"
                 :action #'(lambda () (send self :fit-by-marks 'parallel))))
         (fbgc (send menu-item-proto :new "Fit by marks--equal intercept"
                 :action #'(lambda () (send self :fit-by-marks 'concurrent))))
         (npt (if (send self :has-method :bootstrap1)
	          (send menu-item-proto :new "Lack-of-fit test"
                      :action #'(lambda () (send self :bootstrap1)))))
         (add (send menu-item-proto :new "Add power to slidebar..."
                :action #'(lambda () (send self :add-to-slidebar-dialog))))
         (items (remove nil (mapcar #'(lambda (s) 
                        (send menu-item-proto :new (send s :title) 
                            :action #'(lambda () (send self :setup s)))) 
                                        (send self :smoothers)))))
    (defmeth extract :update ()
       (send extract :enabled (and (> (send control :index) 0)
                                   (null (send control :fit-by-marks)))))
    (mapcar #'(lambda (c m) (defmeth m :update () 
       (send m :enabled (send c :enable-menu-item d)))) 
       (send control :smoothers) items)
    (defmeth fbg :update () 
     (let ((owner (when (send graph :has-slot 'owner)
                        (send graph :slot-value 'owner))))
       (send self :mark (equal (send control :fit-by-marks) 'general))
       (cond ((null owner) (send self :enabled nil))
             (t (send self :enabled (send owner :mark))))))
    (defmeth fbgp :update () 
     (let ((owner (when (send graph :has-slot 'owner)
                        (send graph :slot-value 'owner))))
       (send self :mark (equal (send control :fit-by-marks) 'parallel))
       (cond ((null owner) (send self :enabled nil))
             (t (send self :enabled (send owner :mark))))))
    (defmeth fbgc :update () 
     (let ((owner (when (send graph :has-slot 'owner)
                        (send graph :slot-value 'owner))))
       (send self :mark (equal (send control :fit-by-marks) 'concurrent))
       (cond ((null owner) (send self :enabled nil))
             (t (send self :enabled (send owner :mark))))))
    (defmeth add :update ()
       (send self :enabled (equal (send control :title) "Power curv")))
    (if npt
        (defmeth npt :update () 
         (let* ((np (send graph :find-overlay 'smoother-proto)))
          (send self :enabled (and (> (send control :index) 0)
                                   (null (send control :fit-by-marks))
                                   (null (send np :fit-by-marks))
                                   (> (send np :index) 0))))))
    (append (list extract fbg fbgp fbgc) 
            (append items (if npt (list npt) nil))
            (cons add (send self :other-menu-items)))))
;;;
;;;  Bootstrap nonparametric test for a parametric model
;;;  Based on Hart (1997), Nonparametric Smoohting and Lack-of-Fit Tests
;;;  p. 128 and p. 150.  Works only without point marking.
;;;

;;; duplicate in boot.lsp
(defun boots (n)
"Function args: (n)
Returns a sample with replacement of size n from (iseq n)"
  (floor (* n (uniform-rand n))))

;;; new code follows
(defun np-lof (x res smoother &key (verbose nil))
"Function args: (res smoother)
Computes the nonparametric lack of fit test R_s given on page 150 of 
Hart, J (1997), Nonparametric Smoothing and Lack of Fit Tests, Springer. 
x is a list of design points. res are the residuals from the parametric
model.  smoother is a function such that (funcall #'smoother x res)
returns the smooth.  The denominator of the test is from page 128.  The
normalizing constant a_n is correct for simple linear regression but a 
bit too big otherwise."
  (let* ((s2 (np-var x res))
         (g (funcall smoother res)))
   (if verbose (format t "s = ~a ~15tNum = ~a ~34tR = ~a~%" 
            (sqrt s2) (sum (^ g 2)) (/ (sum (^ g 2)) s2)))
   (/ (sum (^ g 2)) s2)))

(defun np-var (x res)
"Function args: (x res)
Computes a nonparametric estimate of sigmasq, assuming constant variance,
based on the residuals.  A normalizing constant is required, and the one
computed here assumes the parametric fit was simple linear regression.  See
J. Hart (1997), Nonparametric Smoothing and Lack of Fit Tests, Springer, p.
128."
  (/ (sum (^ (- (select res (iseq (1- (length res))))
             (select res (iseq 1 (1- (length res))))) 2))
     (an-slr x)))

; revised 2/5/04 by sw
(defun an-slr (x)
"Function args: (x) computes the constant a_n, page 128 of Hart, for simple
linear regression with design x."
;This quantity is 2(n-1) trace(HQQ') where
;H is a tridiagonal matrix, Q is from the QR factorization of (1, x).  
;Now H1=0, so we need only compute the second column of Q, which is just
;a scaled standardized version of x.
  (let* ((n (length x)) ; sample size
         (z (- x (mean x)))
	 (z (/ z (sqrt (sum (^ z 2))))) ; standardized , column of Q
	 (Hz (+ (- z (combine (first z) (select z (iseq (- n 1)))))
	        (- z (combine (select z (iseq 1 (- n 1))) (last z)))))
	;(k (break "1"))
	 (traceHzz (sum (* Hz z))))
   (- (* 2 (- n 1)) traceHzz)))
   


(defun np-lof-boot (x res smoother &key (nboot 49))
"Function args: (x res smoother &key (nboot 49))
Computes a bootstrap approximation to get a significance level for the
lack of fit test comparing a parametric to a nonparametric model.  From
Hart, J (1997), Nonparametric Smoothing and Lack of Fit Tests, Springer,
the bootstrap is outlined in Section 5.4.3 using the statistic R_s defined
on page 150."
  (let* ((Rs (np-lof x res smoother))
         (boot (mapcar #'(lambda (j)
                  (np-lof x (select res (boots (length res))) smoother))
                             (iseq nboot))))
   (- 1 (/ (position Rs (sort-data (combine Rs boot)))
           (length (combine Rs boot))))))

(arc-setting (QUOTE *nptest-num-bootstraps*)
(QUOTE 99)
"Sets the number of bootstraps to use in computing
a nonparametric test of goodness of fit on 2D graphs."
(QUOTE 99)
)

(defmeth polynomial-proto :bootstrap1 () 
 (format t "Computing bootstrap significance level...~%")
 (send self :bootstrap))
(defmeth polynomial-proto :bootstrap (&key (nboot *nptest-num-bootstraps*))
"Method args: (&key (nboot *nptest-num-bootstraps*))
Computes a lack of fit test for the parametric model currently specified by
the slider, using the smoother currently specified on the nonparametric
slider."
  (let* ((graph (send self :graph))
         (np (send graph :find-overlay 'smoother-proto))
         (showing (send graph :points-showing))
         (n (send graph :num-points))
         (c (first (first (send self :compute))))
         (x1 (send graph :point-coordinate 0 showing))
         (y1 (send graph :point-coordinate 1 showing))
         (nb nboot)
         (id (order x1))
; This method deletes points that are not showing from the computations, but
; the :compute method in the next function also deletes points not showing.
; Consequently, we must expand the list res to the right length before calling
; the compute method.  (3 October 2000)
         (sm #'(lambda (res) 
                (let* ((a (repeat 0 n)))
                 (setf (select a showing) res)
                 (second (first (first (send np :compute :vertical a)))))))
         (ans (np-lof-boot (first c)
                           (- (select y1 id) (second c))
                           sm :nboot nb)))
   (format t "~a~%Lack of fit p-value = ~a based on ~a bootstraps.~%"
     (send graph :title) ans nb)
   (send graph :redraw)
     ans))