;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;    Code for simulations in Arc
;;    Documentation available at http://www.stat.umn.edu/arc/simulate.pdf
;;
;;        Save this file to your disk, and put it in the Extras directory
;;        in the directory from which you start Arc.
;;
;;  $Revision: 1.2 $
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun sim1 (&key 
            (n 100)   ; sample size
            (B 1000)  ; number of simulations
            (eta (list 0 0)) ; values of the regression parameters
            (x-dist #'normal-rand) ; distribution to generate x
            (error-dist #'(lambda (n) (- (uniform-rand n) .5))))
"Function args:  (&key (n 100) (m 1000) (eta (list 0 0))
            (x-dist #'normal-rand)
            (error-dist #'(lambda (n) (- (uniform-rand n) .5))))
Does B simulations of n observations on the simple regression of 
y on x with x sampled from the x-dist distribution.  The true 
values of the parameters is given by eta, and the error 
distribution is given by error-dist."
; Set up initial conditions
  (let* ((eta0 (first eta))
         (eta1 (second eta))
         (x (funcall x-dist n))
         (yhat (+ eta0 (* eta1 x)))
         (e (funcall error-dist n))
         (d (make-dataset :data (list x yhat (+ yhat e))
                          :data-names (list "x" "yhat" "y") 
                          :name "sim")))
; define a regression model, named in this case L1:
      (send d :make-glm :predictors (list "X") 
                        :response "Y" 
                        :type :normal 
                        :name "L1"
                        :weights nil 
                        :offset nil 
                        :trials nil 
                        :mean-function "Identity" 
                        :intercept t
                        :display nil)
; do the simulation
     (transpose (mapcar #'(lambda (j)
          (send L1 :yvar (+ yhat (funcall error-dist n)))
          (/ (send L1 :coef-estimates) 
             (send L1 :coef-standard-errors)))
       (iseq B)))
     ))

(defun sim2 (&key 
            (n 100)   ; sample size
            (B 1000)  ; number of simulations
            p         ; number of predictors
            eta       ; values of the regression coef, including intercept
            (x-dist #'(lambda () (normal-rand (repeat n p)))) 
                      ; distribution to generate multivariate x
            (error-dist #'(lambda (n) (- (uniform-rand n) .5))))
"Function args:  (&key (n 100) (m 1000) eta p
            (x-dist #'normal-rand)
            (error-dist #'(lambda (n) (- (uniform-rand n) .5)))
            (x-dist #'(lambda () (normal-rand (repeat n p)))))
Does B simulations of n observations on the simple regression of 
y on x with x sampled from the x-dist distribution.  The true 
values of the parameters is given by eta, and the error 
distribution is given by error-dist."
; Set up initial conditions
  (let* ((x (funcall x-dist n))
         (yhat (+ (first eta) (sum (* (rest eta) x))))
         (e (funcall error-dist n))
         (r (make-linear :x x :y (+ yhat e))))
; do the simulation
     (transpose (mapcar #'(lambda (j)
          (send r :yvar (+ yhat (funcall error-dist n)))
          (send r :simulated-values))
       (iseq B)))
     ))

(defun sim3 (&key 
            (n 100)   ; sample size
            (B 1000)  ; number of simulations
            p         ; number of predictors
            eta       ; values of the regression coef, including intercept
            (x-dist #'(lambda () (normal-rand (repeat n p)))) 
                      ; distribution to generate multivariate x
            (error-dist #'(lambda (n) (- (uniform-rand n) .5))))
"Function args:  (&key (n 100) (m 1000) eta p
            (x-dist #'normal-rand)
            (error-dist #'(lambda (n) (- (uniform-rand n) .5)))
            (x-dist #'(lambda () (normal-rand (repeat n p)))))
Does B simulations of n observations on the simple regression of 
y on x with x sampled from the x-dist distribution.  The true 
values of the parameters is given by eta, and the error 
distribution is given by error-dist."
; Set up initial conditions
  (let* ((x (funcall x-dist n))
         (yhat (+ (first eta) (sum (* (rest eta) x))))
         (e (funcall error-dist n))
         (r (make-phd :x x :y (+ yhat e))))
; do the simulation
     (transpose (mapcar #'(lambda (j)
          (send r :yvar (+ Yhat (funcall error-dist n)))
          (send r :simulated-values))
       (iseq B)))
     ))

(defmeth inverse-regression-model-proto :simulated-values ()
  (fourth (first (send self :phd-test))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                                  ;;
;;  Functions for simulations and bootstraps in Arc                 ;;
;;  April, 1999                                                     ;;
;;                                                                  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; The purpose of these functions is to allow creating model objects in Arc
; without first creating a dataset.  This will be useful primarily for 
; simulation studies.  None of these functions use any graphical methods,
; like menus, graphs and dialogs.
;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  GLM objects, including normal linear models                     ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;  In each of the functions below, the following arguments may appear:
;    :x gives the list of lists of predictors  (REQUIRED)
;    :y gives the response  (REQUIRED)
;    :trials (binomial only) gives the number of trials  (REQUIRED)
;  The rest are optional:
;    :weights gives the weights
;    :intercept t or nil for an intercept
;    :offset a list of numbers if there is an offset
;    :link gives the link (e.g., probit-link) if not the canonical link
;  You can get a list of the link function names using the function
;          (link-functions :normal)
;          (link-functions :poisson)
;          (link-functions :gamma)
;          (link-functions :binomial)
;  Of course you can write your own link functions as well

(defun make-linear (&rest args) (apply #'make-glm :normal args))

(defun make-binomial (&rest args) (apply #'make-glm :binomial args))

(defun make-poisson (&rest args) (apply #'make-glm :poisson args))

(defun make-gamma (&rest args) (apply #'make-glm :gamma args))

(defun make-glm (type &rest args)
 (apply #'send (type-to-proto type) :new :print nil args))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  Inverse regression objects                                      ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  These are similar to the GLM functions, but they return an inverse
;  regression object.  The arguments are:
;   :x list of lists of predictors REQUIRED
;   :y list of response REQUIRED
;   :nslices, SIR/SAVE only
;  Other args are optional, like
;   :weights
;   :method-response either :yvar or :residuals
;   anything else passed to the :isnew .

(defun make-ir (&rest args &key method info y x
        (method-response :yvar)
        (nslices (floor (/ (length y) 5))))
"Method args: (&rest args)
Method for creating an inverse regression model without using menus."
  (apply #'send inverse-regression-model-proto :new :x x :y y :nslices nslices
    :method-response method-response
    :add-dirs-to-dataset nil :method method :info info args))

(defun make-sir (&rest args)
  (apply #'make-ir :method :sir :info (select *ir-meths* 0) args))

(defun make-phd (&rest args)
  (apply #'make-ir :method :phd :info (select *ir-meths* 1) 
      :method-response :residuals args))

(defun make-save (&rest args)
  (apply #'make-ir :method :save :info (select *ir-meths* 0) args))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;;     nonlinear regression
;; REQUIRED arguments:
;;   :x predictors
;;   :y response
;;   :mean-function a STRING giving the mean function
;;   :starting-values the starting values for the iteration
;; OPTIONAL arguments:
;;   any other args permitted for nonlinear regression
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun make-nonlinear (&rest args &key x y mean-function starting-values)
 (apply #'nrcode :pred-vars x :y y 
   :mean-function (string-to-function mean-function)
   :exact-first-derivative (ignore-errors (string-to-fun-deriv mean-function))
   :theta starting-values
   args))