;;;
;;;  March, 1999 - Bootstrap Methods in Linear Regression.
;;;  September, 2001 - modified for glms, nonlinear models, and gnms
;;;  Iain Pardoe and Sanford Weisberg.
;;;  Documentation at
;;;    http://www.stat.umn.edu/arc/bootmeth.pdf
;;;

;;;
;;;  Download this text file, and save it as boot.lsp in
;;;  the Extras directory in the Arc directory
;;;
;;; $Revision: 1.2 $
(defproto regression-bootstrap-proto '(model boots nboots levels 
                    method) () *object*)

(defmeth regression-bootstrap-proto :isnew (&key model 
     (nboots 999) (levels '(.95))(method :r-boot))
 (send self :model model)
 (send self :nboots nboots)
 (send self :levels levels)
 (send self :method method)
 (call-next-method))

;;; 
;;; slot accessors 
;;; 

(defmeth regression-bootstrap-proto :model (&optional (new nil set))
  (if set (setf (slot-value 'model) new))
  (slot-value 'model)) 

(defmeth regression-bootstrap-proto :nboots (&optional (new nil set))
  (if set (setf (slot-value 'nboots) new))
  (slot-value 'nboots))

(defmeth regression-bootstrap-proto :levels (&optional (new nil set))
  (if set (setf (slot-value 'levels) (if (consp new) new (list new))))
  (slot-value 'levels))

(defmeth regression-bootstrap-proto :method (&optional (new nil set))
  (if set (setf (slot-value 'method) new))  (slot-value 'method))

(defmeth regression-bootstrap-proto :boots (&key (verbose nil)) 
"Method args: (&key (verbose nil))
Returns or computes the bootstraps."
  (if (slot-value 'boots) 
      (slot-value 'boots)
      (setf (slot-value 'boots) 
            (send (send self :model) :compute-bootstrap))))
 
; this works for gnm's and its parents 
(defmeth regression-bootstrap-proto :names ()
   (if (send (send self :model) :intercept)
       (cons "Intercept" (send (send self :model) :predictor-names))
       (send (send self :model) :predictor-names))) 

;;; 
;;; auxillary functions 
;;; 

(defun boots (n)
"Function args: (n)
Returns a sample with replacement of size n from (iseq n)"
  (floor (* n (uniform-rand n))))

(defun permute (n &optional (size n))
"Function args: (n &optional (size n))
Returns a sample without replacement of size s from (iseq n)"
  (select (order (uniform-rand n)) (iseq size))) 

(defun incmissing (incmiss excmiss included) 
"Function args: (incmiss excmiss included) 
Replaces non-missing elements in incmiss with elements in excmiss"
   (if (arrayp incmiss)
       (let* ((incmiss1 (copy-array incmiss))
         (p (iseq (second (array-dimensions incmiss1)))))
    (setf (select incmiss1 included p) excmiss)
    incmiss1)
       (let* ((incmiss1 (copy-list incmiss)))
    (setf (select incmiss1 included) excmiss)
    incmiss1))) 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; regression-model-proto methods 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmeth regression-model-proto :bootstrap (&rest args)
"Method args: (&rest args)
Fills the slot 'bootstrap-object, saving all the args as slot values in 
this object.  The args can be :method, :nboots and :levels."
  (send self :add-slot 'bootstrap-object
   (apply #'send regression-bootstrap-proto :new :model self args))) 

(defmeth regression-model-proto :compute-bootstrap ()
  (let* ((b (slot-value 'bootstrap-object)))
   (send self (send b :method) :nboots (send b :nboots))))    

;;; for nonlinear models, the method returns the number of failed 
;;; reps as the second arg.  If this in non-zero, print it. 
(defmeth nonlin-model-proto :compute-bootstrap 
                      (&rest args &key (method :r-boot))
   (let* ((b (slot-value 'bootstrap-object))
          (ans (send self (send b :method) :nboots (send b :nboots)))
          (z (if ans (apply #'sum (last ans)))))
     (if (> z 0)
         (format t "Warning:  ~a bootstrap replications failed to converge~%"
           z))
    (first ans)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;   Types of bootstraps:
;;;    :r-boot  resample residuals
;;;    :c-boot  resample cases
;;;    :p-boot  parametric bootstrap
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; 
;;;  residual bootstrap.  Works for normal linear and nonlinear models
;;;  Not appropriate for other glms or gnms.
;;;  March, 1999 - computes bootstrap coefficient estimates,  
;;;  resampling centered modified residuals as defined in Davison & 
;;;  Hinkley (1997) p259. 
;;; 

;;; for linear regression only --- uses QR factorization 

(defmeth regression-model-proto :r-boot (&key (nboots 1000)) 
"Method args: (nboots 1000) 
Resamples from the centered modified residuals, see Davison & 
Hinkley (1997) p259, to get nboots pseudo samples, and for each 
sample the coefficient estimates are recomputed.  Returns a list of 
p lists of length nboots, containing the psuedo-estimates."
   (let* ((m (project (send self :x)
                      :included (send self :included)
                      :weights (send self :weights)))
          (f (send self :fit-values))
          (r (/ (send self :residuals)
            (sqrt (- 1 (send self :leverages)))))
          (inc (which (send self :included)))
          (rbar (mean (select r inc)))
          (e (- r rbar))
          (n (length inc))
          (sqrtw (if (send self :weights) (sqrt (send self :weights)) 1)))
    (transpose (mapcar #'(lambda (j)
           (send m :coef-estimates (+ f (/
            (incmissing e (select e (select inc (boots n))) inc)
            sqrtw))))
             (iseq nboots)))))     
        
;;; for normal nonlinear models, uses method's internal compute method. 
(defmeth nonlin-model-proto :r-boot (&key (nboots 1000) (verbose nil)) 
"Method args: (nboots 1000) 
Resamples from the centered modified residuals, see Davison & Hinkley 
(1997) p259, to get nboots pseudo samples, and for each sample the 
coefficient estimates are recomputed.  Returns a list of p lists of 
length nboots, containing the psuedo-estimates."
   (let* ((m (send self :make-clone))
          (f (send self :fit-values))
          (start (send self :slot-value 'theta-hat))
          (r (/ (send self :residuals) (sqrt (- 1 (send self :leverages)))))
          (inc (which (send self :included)))
          (r (select r inc))
          (rbar (mean r))
          (e (- r rbar))
          (num-fail 0)
          (n (length inc))
          (sqrtw (if (send self :weights) (sqrt (send self :weights)) 1)))
    (send m :verbose verbose)
    (let* ((a (transpose (mapcar #'(lambda (j)
                (let ((ans (send m :compute-r-boot e start f inc sqrtw)))
                 (if ans ans (error "Too many convergence failures"))))
                (iseq nboots)))))
     (list (butlast a) (apply #'sum (last a))))))

(defmeth nonlin-model-proto :compute-r-boot (r s f i w &optional (num-fail 0)) 
"Method args: (r s f i w &optional (num-fail 0)) 
r=residuals to be sampled, s=starting values f=fixed fitted vaues, 
i=included vector, w= square root of the weights.  r has been 
shortened to have only included cases.  If num-fail is too big, quit."  
  (when (> num-fail 25) (return-from :compute-r-boot nil))
   (send self :slot-value 'theta-hat s)
   (let ((z (repeat 0 (length f))))
     (setf (select z i) (select r (boots (length r))))
     (setf z (/ z w))     
     (send self :y (+ f z))
     (send self :yvar (send self :y))     
     (send self :coef-estimates)
     (if (> (slot-value 'termination-reason) 0)
         (send self :compute-r-boot r s f i w (1+ num-fail))
         (append (send self :coef-estimates) (list num-fail))))) 

;;; 
;;;  Sept, 2001 - parametric bootstrap
;;; 
(defmeth regression-model-proto :p-boot (&key (nboots 1000))
"Method args:  (nboots 1000)
Uses the parametric bootstrap."
  (let* ((r (send self :make-clone)))
    (send r :verbose nil)
    (transpose (mapcar #'(lambda (j)
                 (send r :p-boot1 self)) (iseq nboots)))))

(defmeth nonlin-model-proto :p-boot (&key (nboots 1000)) 
"Method args: (nboots 1000) 
Gets bootstrap pseudo-estimates from resampling cases in :p-boot1 using 
a clone of the regression model." 
  (let* ((r (send self :make-clone)))
    (send r :verbose nil)
    (let ((a (transpose (mapcar #'(lambda (j) (send r :p-boot1 self))
                                  (iseq nboots)))))
    (list (butlast a) (apply #'sum (last a))))))
                 
(defmeth regression-model-proto :p-boot1 (parent)
    (send self :y (send parent :random-response))
    (send self :yvar (send self :y))
    (send self :compute)
    (send self :coef-estimates)) 

(defmeth nonlin-model-proto :p-boot1 (parent &optional (num-fail 0))
  (when (> num-fail 25)
        (error "Too many convergence failures"))
;   reset starting values to converged values in the parent
   (send self :slot-value 'theta-hat
    (send parent :slot-value 'theta-hat))
   (let* ((ans (call-next-method parent))
          (reason (send self :slot-value 'termination-reason)))
    (if (> reason 0) (send self :p-boot1 parent (1+ num-fail))
        (append ans (list num-fail)))))
     
;;; 
;;;  March, 1999 - computes bootstrap coefficient estimates,  
;;;  resampling cases as defined in Davison & Hinckley (1997) p265. 
;;;  Modified January 13, 2001 to sample correctly for generalized linear 
;;;  models as well as linear models. 
;;; 

(defmeth regression-model-proto :c-boot (&key (nboots 1000)) 
"Method args: (nboots 1000) 
Gets bootstrap pseudo-estimates from resampling cases in :c-boot1 using 
a clone of the regression model." 
  (let* ((r (send self :make-clone)))
    (send r :verbose nil)
    (transpose (mapcar #'(lambda (j) (send r :c-boot1 self)) 
             (iseq nboots)))))

(defmeth nonlin-model-proto :c-boot (&key (nboots 1000)) 
"Method args: (nboots 1000) 
Gets bootstrap pseudo-estimates from resampling cases in :c-boot1 using 
a clone of the regression model." 
  (let* ((r (send self :make-clone)))
    (send r :verbose nil)
    (let ((a (transpose (mapcar #'(lambda (j) (send r :c-boot1 self))
                                  (iseq nboots)))))
    (list (butlast a) (apply #'sum (last a))))))

(defmeth regression-model-proto :c-boot1 (parent)
   (send self :c-boot2 parent))

(defmeth nonlin-model-proto :c-boot1 (parent &optional (num-fail 0))
  (when (> num-fail 25)
        (error "Too many convergence failures"))
;   reset starting values to converged values in the parent
   (send self :slot-value 'theta-hat
    (send parent :slot-value 'theta-hat))
   (let* ((ans (send self :c-boot2 parent))
          (reason (send self :slot-value 'termination-reason)))
    (if (> reason 0) (send self :c-boot1 parent (1+ num-fail))
        (append ans (list num-fail))))) 

;;; this method works for all models, but starting values may be poor
;;; for for gnms.
(defmeth regression-model-proto :c-boot2 (parent) 
"Method args: (parent) 
Resamples cases, see Davison & Hinkley (1997) p265, to get nboots 
pseudo samples, and for each sample the coefficient estimates are 
recomputed. Returns a list of p lists of length nboots, containing 
the psuedo-estimates."
   (let* ((x (send parent :x))
          (p (iseq (second (array-dimensions x))))
          (w (send parent :weights))
          (y (send parent :yvar))
          (inc (which (send parent :included)))
          (trials (when (eq (send self :type) :binomial) (send self :trials)))
          (n (length inc))
          (b (select inc (boots n))))
    (send self :x (incmissing x (select x b p) inc))
    (send self :y (incmissing y (select y b) inc))
    (send self :yvar (send self :y))
    (when trials
          (send self :trials (incmissing trials (select trials b) inc)))
    (if w (send self :weights (incmissing w (select w b) inc)))
    (send self :coef-estimates))) 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;  Summary methods
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;; 
;;;  March, 1999 - histograms of bootstrap coefficient estimates 
;;; 
;;;  eg (send r :histograms) 
;;;  or (send c :histograms) 
;;;  or (send reg :bootstrap-histograms) 
;;; 

(defmeth regression-bootstrap-proto :histograms ()
  (let* ((m (send self :model))
         (d (send m :data))
         (h (send d :m-hist :v (send self :boots)
                            :names (send self :names)
                            :title "Bootstrap coefficient estimates"
                :link-plot nil
                :plot-updating nil)))
   (send h :delete-overlay 
             (send h :find-overlay 'deletion-menu-control-proto))
   h))   

(defmeth regression-model-proto :bootstrap-histograms (&rest args)
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send (send self :slot-value 'bootstrap-object) :histograms)) 

;;; 
;;;  March, 1999 - probability plots of bootstrap coefficient estimates 
;;; 
;;;  eg (send r :probability-plots) 
;;;  or (send c :probability-plots) 
;;;  or (send reg :bootstrap-probability-plots)
;;; 
;;; change the quantile function of glms and gnms to normal. 
(defmeth regression-bootstrap-proto :probability-plots
     (&key (quantile-function
               #'(lambda (a) (t-quant a (send (send self :model) :df))))
           (quantile-name
               (format nil "t_~a quantiles" (send (send self :model)
                      :df))))  
  (let* ((m (send self :model))
         (d (send m :data))
         (coef (send m :coef-estimates))
         (coefse (send m :coef-standard-errors))
         (boot (send self :boots))
         (zboot (/ (- boot coef) coefse))
         (p (send d :probability-plot :values zboot
                    :names (send self :names)
                    :title "Bootstrap probability plots"
            :plot-updating nil
            :link-plot nil
                    :quantile-name quantile-name
                    :quantile-function quantile-function))
     (ov (send p :find-overlay 'deletion-menu-control-proto))
         (c (send bootstrap-control-proto :new 
         (send ov :location))))
    (send p :delete-overlay ov)
    (send p :linked nil)
    (send p :plot-updating nil)
    (send p :add-overlay c)
    (send p :add-slot 'boot self)
    p)) 
    

(defmeth regression-model-proto :bootstrap-probability-plots (&rest args)
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (apply #'send (send self :slot-value 'bootstrap-object) 
             :probability-plots args)) 

;  the following get the right quantile function depending on the proto
(defmeth glim-proto :bootstrap-probability-plots (&rest args)
  (apply #'call-next-method :quantile-function #'normal-quant 
                            :quantile-name "Normal quantiles" args))
(defmeth normalreg-proto :bootstrap-probability-plots (&rest args)
  (apply #'call-method regression-model-proto 
            :bootstrap-probability-plots args))

;;; moved to gnm.lsp
;(defmeth gnm-normal :bootstrap-probability-plots (&rest args)
;  (apply #'call-next-method args))
;(defmeth gnm-proto :bootstrap-probability-plots (&rest args)
;  (apply #'call-next-method :quantile-function #'normal-quant 
;                           :quantile-name "Normal quantiles" args))
;;; end move
(defmeth nonlin-model-proto :bootstrap-probability-plots (&rest args)
  (apply #'call-next-method args))

;;; 
;;;  March, 1999 - linear regression bootstrap standard errors. 
;;; 
;;;  eg (send r :bootstrap-standard-errors) 
;;;  or (send c :bootstrap-standard-errors) 
;;;  or (send reg :bootstrap-standard-errors) 
;;; 
(defmeth regression-bootstrap-proto :bootstrap-standard-errors () 
"Method args: () Computes bootstrap standard errors."
  (mapcar #'standard-deviation (send self :boots)))
    
(defmeth regression-model-proto :bootstrap-standard-errors (&rest args)   
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send (send self :slot-value 'bootstrap-object) :bootstrap-standard-errors)) 
;;; 
;;;  March, 1999 - linear regression bootstrap bias. 
;;; 
;;;  eg (send r :bootstrap-bias) 
;;;  or (send c :bootstrap-bias) 
;;;  or (send reg :bootstrap-bias) 
;;; 
(defmeth regression-bootstrap-proto :bootstrap-bias () 
"Method args: () 
Computes bootstrap bias, which is the difference between the observed
estimate and the average of the bootstrap replications."
  (mapcar #'(lambda (c b) (- c (mean b)))
            (send (send self :model) :coef-estimates)
        (send self :boots)))
        
(defmeth regression-model-proto :bootstrap-bias (&rest args)
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send (send self :slot-value 'bootstrap-object) :bootstrap-bias)) 
;;; 
;;;  March, 1999 - linear regression bootstrap confidence intervals. 
;;; 
;;;  eg (send r :display-confidence-intervals) 
;;;  or (send c :display-confidence-intervals) 
;;;  or (send reg :display-confidence-intervals) 
;;; 
(defmeth regression-bootstrap-proto :display-confidence-intervals 
  (&key levels)
"Method args: (&key levels)
If levels is set, the levels slot is changed. Prints confidence 
intervals using three methods for each level and for each regression 
coefficient."
 (if levels (send self :levels levels))
 (let* ((n-theory (send self :n-theory-conf-intervals))
        (percent-boot (send self :percentile-boot-intervals))
        (bca-boot (send self :BCa-bootstrap-intervals)))
   (format t "Level ~7tNormal Theory ~29tPercentile Bootstrap~52tBCa Bootstrap~%")
   (mapcar #'(lambda (name nt pb bcb)
          (format t "Coefficient: ~a~%" name)
          (mapcar #'(lambda (lev nt1 pb1 bcb1)
            (format t "~4f~7t~21a~29t~21a~52t~21a~%" 
               lev nt1 pb1 bcb1)) 
             (send self :levels) nt pb bcb))
      (send self :names) (third n-theory) (third percent-boot)
      (third bca-boot))
   t))

(defmeth regression-model-proto :bootstrap-confidence-intervals (&rest args)
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send self :bootstrap-description)
   (send (send self :slot-value 'bootstrap-object)
                 :display-confidence-intervals)) 
;;; different quantile function for gnms, otherwise should be OK. 
(defmeth regression-model-proto :n-theory-conf-intervals
   (&key (levels '(.95))) 
"Method args: (&key levels) Returns a list of three elements:
(list levels, term names, confidence levels), giving normal theory (t-based) 
confidence intervals for each term in the linear model, at each level specified.  
Levels should be a list."
  (let* ((levels (if (consp levels) levels (list levels)))
         (alpha (/ (- 1 levels) 2))
     (wald-mult (- (send self :wald-quantile alpha (send self :df))))
         (est (send self :coef-estimates))
         (se (send self :coef-standard-errors))
         (names (if (send self :intercept)
                    (cons "Intercept" (send self :predictor-names))
                    (send self :predictor-names)))
         (ci (mapcar #'(lambda (coef se)
                       (mapcar #'(lambda (a) 
                      (list (- coef (* a se))
                                    (+ coef (* a se)))) wald-mult)) 
              est se)))
   (list levels names ci)))

(defmeth regression-model-proto :wald-quantile (alpha df) (t-quant alpha df))
(defmeth glim-proto :wald-quantile (alpha df) (normal-quant alpha))
(defmeth normalreg-proto :wald-quantile (alpha df) (t-quant alpha df))
;;; moved to gnm.lsp
;(defmeth gnm-proto :wald-quantile (alpha df) (normal-quant alpha))
;(defmeth gnm-normal :wald-quantile (alpha df) (t-quand alpha df))
;;; end move

(defmeth regression-bootstrap-proto :n-theory-conf-intervals () 
"Method args: () 
Computes normal theory confidence intervals."
  (send (send self :model) :n-theory-conf-intervals
        :levels (send self :levels))) 
;; OK 
(defmeth regression-bootstrap-proto :percentile-boot-intervals () 
"Method args: () 
Computes percentile bootstrap confidence intervals."
   (let* ((boot (send self :boots))
          (alpha (/ (- 1 (send self :levels)) 2)))
    (list (send self :levels) (send self :names) (mapcar #'(lambda
    (coef) (mapcar #'(lambda (a)
    (quantile coef (list a (- 1 a)))) alpha)) boot)))) 

;;;; Check Efron and Tibshirani to see if this applies. 
(defmeth regression-bootstrap-proto :BCa-bootstrap-intervals ()
"Method args: ()
Computes BCa bootstrap confidence intervals, from Efron and
Tibshirani (1993), Section 14.3."
 (let* ((m (send self :model))
        (coef (send m :coef-estimates))
        (boot (send self :boots))
        (alpha (/ (- 1 (send self :levels)) 2))
        (bias-levels (transpose (mapcar #'(lambda (c) (transpose
          (send m :bias-corrected-percentiles boot :coverage c))) 
          (send self :levels)))))
   (list (send self :levels) (send self :names) (mapcar #'(lambda 
   (coef b) (mapcar #'(lambda (a) 
   (quantile coef a)) b)) boot bias-levels))))

(defmeth regression-model-proto :bias-corrected-percentiles 
  (boots &key (coverage .95))
"Method args: (boots &key (coverage .95))
Returns p_l and p_u such that an unbiased estimate of a confidence 
interval with coverage coverage is given by the interval from the p_l 
quantile of the bootstrap distribution to the p_u quantile.  From 
Efron and Tibshirani (1993), p185, equations 14.9 and 14.10."
  (let* ((alpha (/ (- 1 coverage) 2))
         (z0 (send self :bootstrap-bias-correction boots))
         (ahat (send self :std-kurtosis))
         (zalpha (normal-quant alpha))
         (zalpha1 (normal-quant (- 1 alpha)))
         (a1 (normal-cdf (+ z0 (/ (+ z0 zalpha)
                                  (- 1 (* ahat (+ z0 zalpha)))))))
         (a2 (normal-cdf (+ z0 (/ (+ z0 zalpha1)
                                  (- 1 (* ahat (+ z0 zalpha1))))))))
   (list a1 a2)))

(defmeth regression-model-proto :bootstrap-bias-correction (boots)
"Method args: (boots) 
Computes Phi^{-1}((# betahat_boot < betahat)/B), where B is the 
number of bootstraps. From Efron and Tibshirani (1993), equation 
14.14"
  (let* ((count (mapcar #'sum 
        (mapcar #'(lambda (a est) (if-else (< a est) 1 0))
           boots (send self :coef-estimates)))))
    (normal-quant (/ count (length (first boots))))))

(defmeth regression-model-proto :std-kurtosis ()
"Method args: ()
Computes (14.15) from Efron and Tibshirani (1993), which is a 
leave-one-out standardized kurtosis used in the BCa bootstrap."
  (let* ((basis (if (send self :intercept) (cons 0 (1+ 
      (send self :basis))) (send self :basis)))
         (xtxinv (send self :xtxinv))
         (inc (which (send self :included)))
         (w (if (send self :weights) (select (send self :weights) 
                    inc) 1))
         (x (select (send self :x-matrix) inc basis))
         (v (mapcar #'(lambda (a) (matmult xtxinv a))
              (* (row-list x) (select (send self :residuals) inc) w
                 (/ (- 1 (select (send self :leverages) inc))))))
         (vbar (/ (apply #'+ v) (length inc)))
         (vc (mapcar #'(lambda (v) (- v vbar)) v)))
  (/ (apply #'+ (^ vc 3)) (* 6 (^ (apply #'+ (^ vc 2)) 1.5)))))

;;;
;;;  March, 1999 - linear regression bootstrap confidence intervals 
;;;  for studentized estimates. 
;;; 
;;;  eg (send r :stud-percentile-boot-intervals) 
;;;  or (send c :stud-percentile-boot-intervals) 
;;; 
(defmeth regression-bootstrap-proto :stud-percentile-boot-intervals 
  (&key levels)
"Method args: (&key levels)
Computes percentile bootstrap confidence intervals for studentized 
estimates."
  (if levels (send self :levels levels))
  (let* ((m (send self :model))
         (coef (send m :coef-estimates))
         (boot (send self :boots))
         (bootse (mapcar #'standard-deviation boot))
         (zboot (/ (- boot coef) bootse))
         (alpha (/ (- 1 (send self :levels)) 2)))
   (format t "Level ~8tPercentile Bootstrap of (b*-bhat)/boot(se)~%")
   (mapcar #'(lambda (name pb)
    (format t "Coefficient: ~a~%" name)
    (mapcar #'(lambda (lev pb1)
     (format t "~4f~8t~21a~%" lev pb1)) (send self :levels) pb))
   (send self :names) (mapcar #'(lambda (z) (mapcar #'(lambda (a)
           (quantile z (list a (- 1 a)))) alpha)) zboot))
   t))
   
(defmeth regression-model-proto :bootstrap-percentile-intervals (&rest args)
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send self :bootstrap-description)
   (send (send self :slot-value 'bootstrap-object) :stud-percentile-boot-intervals)) 
;;; 
;;;  March, 1999 - linear regression bootstrap two-sided p-values 
;;;  for studentized estimates. 
;;; 
;;;  eg (send r :stud-boot-pvalues) 
;;;  or (send c :stud-boot-pvalues) 
;;; 
(defmeth regression-bootstrap-proto :stud-boot-pvalues () 
"Method args: () 
Computes bootstrap two-sided p-values for studentized estimates."
   (let* ((m (send self :model))
          (coef (send m :coef-estimates))
          (coefse (send m :coef-standard-errors))
          (zcoef (/ coef coefse))
          (boot (send self :boots))
          (bootse (mapcar #'standard-deviation boot))
          (zboot (/ (- boot coef) bootse))
          (nboot (send self :nboots)))
    (mapcar #'(lambda (z1 z2)
      (/ (sum (if-else (> (abs z1) (abs z2)) 1 0)) nboot)) zboot zcoef)))

(defmeth regression-model-proto :stud-boot-pvalues ()
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send (send self :slot-value 'bootstrap-object) :stud-boot-pvalues)) 
    
(defmeth regression-bootstrap-proto :description ()
  (let* ((meth (if (equal (send self :method) :r-boot)
                   "Bootstrap based on resampling residuals"
                   "Bootstrap based on resampling cases")))
   (format t "~a~%Number of replications = ~a~%" meth (send self :nboots))))
   
(defmeth regression-model-proto :bootstrap-description ()
 (when (send self :has-slot 'bootstrap-object)
  (send self :description)
  (send (slot-value 'bootstrap-object) :description)))
    
(defmeth regression-model-proto :bootstrap-Wald-tests (&rest args)
   (unless (and (null args) (send self :has-slot 'bootstrap-object))
           (apply #'send self :bootstrap args))
   (send self :bootstrap-description)
   (send (send self :slot-value 'bootstrap-object) :stud-boot-pvalues)) 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;; 
;;; Checked to here 9/15/01 
;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;;
;;;  March, 1999 - linear regression bootstrap p-value 
;;;  for testing null model nested in alternative.
;;; 
;;;  eg (send L1 :sigtest-boot-pvalue L2 :nboots 1000 :method :r-boot)
;;;

(defmeth regression-model-proto :sigtest-boot-pvalue 
  (submodel &key (nboots 999) (method :r-boot))
"Method args: (submodel &key (nboots 1000))
Computes bootstrap p-value for testing null model nested in 
alternative."
  (let* ((xmatrix (send self :permutex (send self :x) submodel))
    (proj (project xmatrix 
         :weights (send self :weights)
         :included (send self :included)
         :intercept (send self :intercept))) 
    (subp (length (send submodel :basis)))
    (inc (which (send self :included)))
    (yfull (send self :y))
    (fullfit (send self :fit-values))
    (fullres (send self :residuals))
    (fulldf (send self :df))
    (fullrms (/ (sum (select (^ fullres 2) inc)) fulldf))
    (zsample (/ (sum (select (^ (send proj :p2.1 yfull subp) 2) inc)) 
           fullrms))
    (fullresadj (/ fullres (with-missing #'sqrt 
                     (- 1 (send self :leverages)))))
    (rbar (mean (select fullresadj inc)))
    (e (- fullresadj rbar))
    (n (length inc))
    (sqrtw (if (send self :weights) (sqrt (send self :weights)) 1))
    (ind (if (equal method :r-boot)
        (if-else (> (mapcar #'(lambda (j) 
     (let* ((yboot (+ fullfit (/ (incmissing e (select e 
              (select inc (boots n))) inc) sqrtw)))
       (diff (- yboot yfull))) 
       (/ (sum (select (^ (send proj :p2.1 diff subp) 2) inc)) 
     (/ (sum (select (^ (send proj :residuals yboot) 2) inc)) fulldf)))) 
                  (iseq nboots)) zsample) 1 0)
        nil)))
   (if (equal method :r-boot) 
       (format t "Residual bootstrap~%p-value: ~,4f~%" (/ (sum ind) nboots))
       (format t "Case bootstrap not appropriate for~%testing null model ~
        nested in alternative~%"))
   t))

(defmeth regression-model-proto :permutex (x submodel)
"Method args: (x submodel)
Re-orders the columns of x so the submodel columns come first"
  (let* ((newx (copy-array x))
    (removed 0)
    (n (send self :num-included))
    (fullp (length (send self :basis)))
    (subp (length (send submodel :basis))))
    (dolist (j (iseq fullp))
      (cond ((member (select (send self :predictor-names) j) 
           (send submodel :predictor-names)) 
        (setf (select newx (iseq n) (- j removed)) 
         (select x (iseq n) j)))
       (t 
        (setf (select newx (iseq n) (+ subp removed)) 
         (select x (iseq n) j))
        (setf removed (1+ removed)))))
    newx))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; interface stuff
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defparameter *bootstrap-methods*
 (list (list :c-boot "Case resampling")
       (list :r-boot "Residual resampling")
       (list :p-boot "Parametric bootstrap")))

(defun bootstrap-dialog (&optional (methods *bootstrap-methods*))
"Gets necessary information to setup the bootstrap."
  (let* ((title (send text-item-proto :new (format nil 
"Setup the bootstrap.  This dialog will~%create graph with a menu and more options.")))
         (method (send text-item-proto :new "Bootstrap method"))
         (get-method (send choice-item-proto :new (mapcar #'second methods)))
         (num (send text-item-proto :new "Number of bootstraps"))
         (get-num (send edit-text-item-proto :new "999" 
                  :text-length (apply #'ms *text-width-short*)))
         (lev (send text-item-proto :new "Confidence levels"))
         (get-lev (send edit-text-item-proto :new ".90 .95 .99"
                  :text-length (apply #'ms *text-width-long*)))
         (cancel (send modal-button-proto :new "Cancel"))
         (get-answer #'(lambda () 
           (list :method (select (mapcar #'first methods) 
                         (send get-method :value))
                 :nboots
                  (with-input-from-string (s (send get-num :text)) (read s))
                 :levels
                  (let* ((a (string-trim " " (send get-lev :text)))
                         (c1 (cond ((> (length a) 0)(select a 0)) (t nil)))
                         (ans (cond ((or (equal c1 #\') (equal c1 #\()) ;)
                                     (with-input-from-string (s (send get-lev :text))
                                     (eval (read s))))
                                    ((equal a "") '(.95))
                                    (t (mapcar #'string-to-data 
                                       (divide-string (send get-lev :text)))))))
                  ans)
                 )))
         (ok (send modal-button-proto :new "OK" :action 
                     #'(lambda () (funcall get-answer))))
         (help (send modal-button-proto :new "Help" :action
                     #'(lambda () (bootstrap-dialog-help methods))))
         (dialog (send modal-dialog-proto :new
                           (list (list title)
                            (list method get-method)
                            (list num get-num)
                            (list lev get-lev)
                            (list ok cancel help)))))
        (send dialog :default-button ok)
        (send dialog :modal-dialog)))
        
(defun bootstrap-dialog-help (method)
(big-message-dialog
    "The bootstrap uses simulation to approximate properties of estimates.
Select the bootstrap method and the number of replications.  For each
replication, a \"bootstrap sample\" is generated, and the model is fit.  
The parameter estimates are saved, and available for further use."
    "The confidence levels are used in summarizing the bootstrap."
    )
    (bootstrap-dialog method))

(rc-menu-item 'bootstrap-item "Bootstrap..." :start-bootstrap)

(defun bootstrap-menu-item ()
 (let* ((pos (position 'prediction-dialog *arc-linear-regression-menu-items*))
        (in  (position 'bootstrap-item *arc-linear-regression-menu-items*))
        (pos1 (position 'emod-menu-item *arc-glm-regression-menu-items1*))
        (in1 (position 'bootstrap-item *arc-glm-regression-menu-items1*))
        (pos2 (position 'delta-method-menu-item 
                        *arc-nonlinear-regression-menu-items*))
        (in2 (position 'bootstrap-item *arc-nonlinear-regression-menu-items*)))
  (unless in
   (defparameter *arc-linear-regression-menu-items*
     (list-insert *arc-linear-regression-menu-items* 
                  'bootstrap-item (1+ pos)))) 
  (unless in2
   (defparameter *arc-nonlinear-regression-menu-items*
     (list-insert *arc-nonlinear-regression-menu-items* 
                  'bootstrap-item (1+ pos2))))
  (unless in1
   (defparameter *arc-glm-regression-menu-items1*
      (list-insert *arc-glm-regression-menu-items1*
                  'bootstrap-item (1+ pos1))))))

(bootstrap-menu-item)
 
(defmeth regression-model-proto :start-bootstrap ()
  (let* ((ans (bootstrap-dialog)))
    (when ans 
          (apply #'send self :bootstrap ans)
          (send self :bootstrap-probability-plots)
          (slot-value 'bootstrap-object))))
          
(defmeth normalreg-proto :start-bootstrap () 
  (call-method regression-model-proto :start-bootstrap))
          
(defmeth glim-proto :start-bootstrap ()
  (let* ((ans (bootstrap-dialog (select *bootstrap-methods* '(0 2)))))
    (when ans 
          (apply #'send self :bootstrap ans)
          (send self :bootstrap-probability-plots)
          (slot-value 'bootstrap-object))))


(defproto bootstrap-control-proto '() () popup-menu-control-proto)

(defmeth bootstrap-control-proto :isnew (loc)
  (call-next-method loc :title "Bootstrap"))

(defmeth bootstrap-control-proto :menu ()
"Method args : ()
Returns the bootstrap menu."
  (let* ((graph (send self :graph))
         (boot (send graph :slot-value 'boot)))
    (when (null (slot-value 'menu))
     (let  ((menu (send menu-proto :new 
              "Bootstrap methods menu"))
            (hist (send menu-item-proto :new "Show histograms"
                 :action #'(lambda () (send boot :histograms))))
            (disp (send menu-item-proto :new "Display summaries" 
                 :action #'(lambda () 
             (send boot :display-summaries))))
            (save (send menu-item-proto :new "Save as a dataset" 
                 :action #'(lambda () (send boot :bootstrap-save)))))
       (send menu :append-items disp hist save)
       (setf (slot-value 'menu) menu)))
       (slot-value 'menu)))

(defmeth regression-bootstrap-proto :bootstrap-save ()
"Method args: ()
Creates a dataset from a bootstrap."
 (let ((b self)
       (model (send self :model)))
  (arc :data (send b :boots)
       :data-names (send b :names)
       :description (format nil "Bootstrap summary from model ~a, dataset ~a"
               (send model :name)
               (send (send model :data) :name))
       :name (format nil "B~a" (send model :name) ))))

(defmeth regression-bootstrap-proto :display-summaries ()
  (let* ((m (send self :model))
         (name (mapcar #'(lambda (n) (substr n 0 10))  (send self :names)))
     (est  (send m :coef-estimates))
     (bias  (send self :bootstrap-bias))
     (se (send m :coef-standard-errors))
     (bse (send self :bootstrap-standard-errors))
     (theop (send m :coef-p-values))
     (bootp (send self :stud-boot-pvalues)))
  (send m :bootstrap-description)
  (format t "~13tObserved~25tBootstrap~40tModel")
  (format t "~51tBootstrap ~63tModel ~72tBoot~%")
  (format t "~13tEstimate~25tBias~40tSE")
  (format t "~51tSE~63tpvalue~72tpvalue~%")
  (mapcar #'(lambda (name est bias se bse theop bootp)
    (format t "~10a~13,5g~13,5g~13,5g~13,5g~8,4f~8,4f~%"
                     name est bias se bse theop bootp))
        name est bias se bse theop bootp)
  (format t "Coefficient confidence intervals, using three methods~%")
  (send self :display-confidence-intervals)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;  median bootstrap 
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun parametric-bootstrap (n &key (dist #'normal-rand)
                                    (statistic #'median)
                                (B 999))
"Function args: (n &key (dist #'normal-rand) (statistic #'median)
                        (B 999))
Generate samples of size n from the distribution dist.  On each sample,
compute the statistic statistic.  Repeat B times."
 (mapcar #'(lambda (j)
   (funcall statistic (funcall dist n))) (iseq B)))

(defun case-bootstrap (y &key (statistic #'median) (B 999))
"Functon args:  (y &key (statistic #'median) (B 999))
Performs B bootstraps reampling cases for the specified statistic."
  (parametric-bootstrap (length y)
     :dist #'(lambda (n) (select y (boots n)))
     :statistic statistic :B B))