#|

    class9.lsp     The EM algorithm 

    - Start with debugging examples, lowess simulation
    
    - Genetic example for EM
    
    - Mixtures

|#


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Debugging tools                                     ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;--- Stopping execution via (break) ... (continue)
(dotimes (i 10)
         (setf x (* i 10))
         (break)
         (setf y (+ 3 x))
         )


;--- Tracing, (trace/untrace f)
(defun f (n)
  (if (< 0 n)
      (+ 1 (f (1- n)))
      0))

(trace f)

(f 5)

(untrace f)


;--- Stepping (Next, Break, Help, Skip, Evaluate)

(step (f 2))


;--- (debug) and (baktrace)

(defun f (x)
  (+ x 3)
  (/ x 0)
  )

(f 4)

(debug)   

(baktrace)

(f 4)

(nodebug)


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Genetic example                                     ;;
;;                                                         ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(def y '(125 18 20 34))
(def x '(x0 x1 18 20 34))


;--- define some useful tools
  (defmacro WHILE (test &body body)
    `(do ()
         ((not ,test))
         ,@body))

  (defun NOT-CLOSE? (x y)
    (> (abs (- x y)) .00001))


;--- direct code

(defun ESTIMATE (p)
  ; Adjust the x counts
  (let ((sum (+ .5 (/ p 4)))  )
    (setf (select x 0) (* (/ .5 sum) (first y)))
    (setf (select x 1) (* (/ (/ p 4) sum) (first y)))
    )) 
  
(defun MAXIMIZE (x)
  ; compute new prob estimate
  (setf p^ (/ (+ (select x 1) (select y 3))
              (sum (select x 1) (select y '(1 2 3)))))
  )
 
;--- Do the calculation                 
(setf oldP 0)
(def p^ .5)
(while (not-close? oldP p^)
       (format t "P^ = ~a~%" p^)
       (setf oldP p^)
       (estimate p^)
       (maximize x)
       )



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;     Normal mixtures example                             ;;
;;                           (fun with values)             ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- First, some code to generate the mixtures
   (defun MIXTURE-RAND (n m s)
     "Generate a sample of n observations from a normal mixture
      with parameters in the lists m and s"
     (let* ((k (length m))  ; and s
            (z (truncate (* k (uniform-rand n))))   )
       (setf trueZ z)
       (+ (select m z)
          (* (select s z) (normal-rand n))
          )))

(setf p (plot-lines
         (kernel-dens
          (setf y (mixture-rand 100 '(0 2) '(1 1))))) ) )

(send p :add-lines (kernel-dens y))

*plot-symbols*

(send p :add-points y (repeat 0.01 (length y)))
(send p :point-symbol (iseq (length y)) 
      (select '(cross square) truez))


;--- EM steps

(setf y  '(0 1 2 3 4 5))
(setf sd '(1 1))

   (defun ESTIMATE-Z-FROM-Y (y m s)
     ; Compute prob for each population, the E step
     (transpose
      (mapcar #'(lambda (p) (/ p (sum p)))
              (transpose
               (mapcar #'(lambda (mean sd)
                           (normal-dens (/ (- y mean) sd)))
                       m s))
              )))

(setf z^ (estimate-z-from-y y '(-2 2) sd))

   (defun ESTIMATE-M-FROM-DATA (y z)
     ; Compute new estimates of the mean and standard deviation
     (mapcar #'(lambda (wts)
                 (/ (sum (* y wts)) (sum wts)))
             z))

(setf parm^ (estimate-m-from-data y z^))

      
;--- Do some iterations

(def means '(-2 2))
(dotimes (i 5)
         (setf z^ (estimate-z-from-y y means sd))
         (setf means (estimate-m-from-data y z^))
         (format t "~d --- ~%" i)
         (format t "Means = ~{ ~7,3f ~}~%" means)
         )
         
(plot-points (first z^) trueZ :variable-labels '("Prob in First" "Group"))