From jbond Wed Jun 15 15:26 PDT 1994
Return-Path: <jbond>
Received: from oak.math.ucla.edu by laplace.math.ucla.edu 
	(Sendmail 5.0/1.10) id AA28125; Wed, 15 Jun 1994 15:26:00 +0800
Received: by oak.math.ucla.edu 
	(Sendmail 4.1/1.10) id AA23074; Wed, 15 Jun 94 15:25:56 PDT
Date: Wed, 15 Jun 94 15:25:56 PDT
From: Jason Bond <jbond>
Message-Id: <9406152225.AA23074@oak.math.ucla.edu>
To: afsharto
Content-Type: text
Content-Length: 5111
Status: RO


(defun start (data num)
"Args: DATA.
Data must be in the form of a matrix with variables as columns and cases
as rows.  This function"
 (let* (
       (collist (column-list data))
       (rowlist (row-list data))
       (l1 (length (first rowlist)))
       (l2 (length (first collist)))
       (presind (mapcar #'which collist))
       (presindt (mapcar #'which rowlist))
       (missindt (mapcar #'(lambda (x) (which
                           (mapcar #'not (combine x)))) rowlist))
       (meanlist (mapcar #'(lambda (x y) (mean (select x y))) collist presind))
       (collist (mapcar #'(lambda (x y) (if-else x x y)) collist meanlist))
       (sigma (apply #'covariance-matrix collist))
       (temp (make-array (list l2 l1)))
       (mu (map 'vector #'mean collist))
       (newsigma 0)
       (likehist ())
       (outprod 0)
       )
       (dotimes (j num)

         (dotimes (i l2)
           
           (setf torf (not (not (and (elt missindt i) (elt presindt i)))))
           (setf params (param-est sigma collist (elt missindt i)
                                     (elt presindt i) mu i torf l1))
           (setf newsigma (+ newsigma (first params) (second params)))
           (setf outprod (+ outprod (first params)))
           (cond (torf (setf rowlisti (impute params (elt missindt i)
                                         (elt presindt i) mu)))
                 ((not (not (elt missindt i))) (setf rowlisti (combine mu)))
                 (t (setf rowlisti (combine (last params)))))
           (mapcar #'(lambda (x y) (setf (aref temp i x) y))
                                   (iseq l1) rowlisti)
         )
         (setf collist (column-list temp))
         (setf sigma (/ newsigma l2)) 
         (setf mu (map 'vector #'mean collist))

         (def sinv (inverse sigma)) 

         (setf loglike (list (- (* l2 (determinant sinv))
                          (sum (diagonal (matmult sinv outprod))))))

         (setf likehist (append likehist loglike))
;;         (print loglike)
         (setf newsigma 0)
         (setf outprod 0)

       )

        (print likehist)
        (print "Mu:")
        (print mu)
        (print "Sigma:")
        (print-matrix sigma)
        (print "Data:")
        (print-matrix temp)
   (plot-points (iseq (length likehist)) likehist)
   data
 )
)


(defun param-est (sigma data missind presind mu i torf len)
"Args: (DATA MISSIND PRESIND)
Takes in the data, missing and present indicies and estimates the parameters
of a multivariate normal distribution.  Data must be in the form of a list
of vectors.  Missind must be in the form of a list for case 1.  Presind must
also be in the form of a list for case 2."
(cond (torf
      (let* (
            (sig (inverse sigma))
            (smm-1 (inverse (apply #'bind-rows
                   (mapcar #'(lambda (x) (mapcar #'(lambda (y) (aref sig
                              x y)) missind)) missind))))

            (sooinv (inverse (apply #'bind-rows (mapcar #'(lambda (x)
                 (mapcar #'(lambda (y) (aref sigma x y)) presind)) presind))))
            (smo (apply #'bind-rows (mapcar #'(lambda (x) (mapcar
                 #'(lambda (y) (aref sigma x y)) presind)) missind)))
;;            (uglymat (- smm (matmult smo sooinv (transpose smo))))
            (vi (make-vi smm-1 missind len))
            (datai (mapcar #'(lambda (x) (elt x i)) data))
            (yminusmu (- datai mu))
            )
          (list (outer-product yminusmu yminusmu) vi 
                                            smo sooinv datai)
     )
     )
      (t   (let* (
                   (datai (mapcar #'(lambda (x) (elt x i)) data))
                   (yminusmu (- datai mu))
                 )
                   (cond ((not (not missind))
                          (list (outer-product yminusmu yminusmu) 
                                    sigma))
                   (t    (list (outer-product yminusmu yminusmu) 0 datai))
                   )
           )
      )
)
)




(defun impute (params missind presind mu)
"Args: (DATA PARAMS MISSIND PRESIND)
Imputes the missing data vector for case i indexed my missind by
regressing the missing values on the non-missing values.  Returns the
list of imputed and observed values for case i."
 (let* (
       (lparams (combine (last params)))
       (newvals (+ (select mu missind)
                   (inner-product
                      (matmult (select params 2) (select params 3))
                      (- (select lparams presind)
                         (select mu presind)))))
       )
  (setf (select lparams missind) newvals)
  lparams
 )
)


(defun make-vi (uglymat missind p)
"Args: Uglymat missind p.
Returns a matrix of zeros with the ugly covariance matrix in place of
the missing variables."

 (let* (
       (vi (make-array (list p p) :initial-element 0))
       )
       (mapcar #'(lambda (x y) (mapcar #'(lambda (z w)
                                  (setf (aref vi x z) (aref uglymat y w)))
                                   missind (iseq (length missind))))
                                   missind (iseq (length missind)))
  vi
 )
)