#|

    class8.lsp     Simulation Programming in LispStat 


    - Random number generators

    - Setting the random seed

    - Simulation methods, generic "simulator"

    - Bootstrap resampling

|#

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;  Random number generators                               ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(apropos '-rand)

(apropos 'poisson)



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;  Setting the random seed                                ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


*random-seed*

(setf saveSeed *random-seed*)

(uniform-rand 4)

(setf saveSeed (make-random-seed *random-seed*))

(uniform-rand 4)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;  Simulating the mean and standard deviation             ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Standard setup  (might want to add some printed output)

(setf xbar ())
(setf sd   ())
(setf seed (make-random-state *random-state*))
(dotimes (i 200)
         (let ((x (normal-rand 25)) )   ; avoid global x variable  
           (push (mean x) xbar)
           (push (standard-deviation x) sd)
           ))

(plot-points xbar sd)



;--- Generic format is quite clear
(setf stat1 ())
(setf stat2 ())
(setf seed (make-random-state *random-state*))
(dotimes (i nReps)
         (let ((x ( xxx-rand n))  )
           (push (compute stat1 x ) stat1)
           (push (compute stat2 x ) stat2))
         )



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;       Generic simulator                                 ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


(defun SIMULATE (stat nReps gen)
  "Args (stat nReps gen): Evaluate function stat nReps times
   for samples created by function gen"
  (let ((results ())  )
    (dotimes (i nReps)
             (push (funcall stat (funcall gen))
                   results))
    results))


(simulate #'mean 20 #'(lambda () (normal-rand 25)))


;--- What do you do if you have more than two statistics? (pair vs two-sample)

(defun SIMULATE (stat nReps gen)
  "Args (stat nReps gen): Evaluate function(s) in stat nReps times
   for samples created by function gen"
  (if (listp stat) ; recursion for list of arguments
      (simulate #'(lambda (x) 
                    (mapcar #'(lambda (s) (funcall s x)) stat))
                nReps gen)
      (let ((results ())  )
        (dotimes (i nReps)
                 (push (funcall stat (funcall gen))
                       results))  
        results)
      ))


;--- Gaussian structure
  (def simres
       (simulate (list #'mean #'standard-deviation)
                 20
                 #'(lambda () (normal-rand 30))))

  (plot-points (transpose simres)) 


;--- For Poisson

  (def simresp
       (simulate (list #'mean #'standard-deviation)
                 20
                 #'(lambda () (poisson-rand 30 2))))

(plot-points (transpose simresp))  ; don't look so independent now!



;--- Function that builds a simulator of list of statistics
;    This is does what the other does, just a but differently.
;    It returns a function that is used to simulate the stat.

(defun SIMULATOR (stats)
  "Args (stats): Returns a function of two arguments (gen nReps)
   gen generates the samples for each of nReps calls.  Calls
   the functions in the list stats on samples generated by gen."
  (lambda (gen nReps)
    (let ((results ())
          )
      (dotimes (i nReps)
               (let ((x (funcall gen))  )
                 (push 
                  (mapcar #'(lambda (s) (funcall s x)) stats)
                  results)))
      results)))

(def sim (simulator (list #'mean #'standard-deviation)))
     

;--- Do it for Gaussian and plot
(def simres
     (funcall sim 
              #'(lambda () (normal-rand 20))
              200))

(plot-points (transpose simres))   ; often need list of cols, not rows


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;  Bootstrap calculculations                              ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- sampler

(sample (iseq 10) 10)     ; without replacement

(sample (iseq 10) 10 t)   ; WITH replacement



;--- generic bootstrap program uses our prior simulation routine

(defun BOOTSTRAP (stat data nReps)
  (let ((n (length data))  )
    (simulate stat nReps 
              #'(lambda () (sample data n t))
              )))
  

(bootstrap #'mean (iseq 10) 30)