(defun start (data num)
"Args: (DATA NUM)
This function starts the imputation and parameter estimation using   
multivariate normal distributions.  Data must be in the form of a matrix
with variables as columns and cases as rows.  NUM should be the number
of iterations that are to be carried out...careful..this could take a while."
 (let* (
       (collist (column-list data))
       (rowlist (row-list data))
       (l1 (length collist))
       (l2 (length rowlist))
       (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))
       (siginv (identity-matrix l1)) 
       (mu (map 'vector #'mean collist))
       (likehist ())
       (bool nil)
       )
       (dotimes (j num)
;;       (do
        (setf newsigma 0)
        (setf outprod 0)
        (setf outprodsum 0)
        (setf rowlisti 0)
        (setf params 0)
        (dotimes  (i l2)
           (setf torf (not (not (and (elt missindt i) (elt presindt i)))))
           (setf params (param-est siginv collist (elt missindt i)
                                     (elt presindt i) mu i torf l1))
           (cond (torf (setf rowlisti (impute params (elt missindt i)
                                         (elt presindt i) mu collist i)))
                 ((not (not (elt missindt i))) (setf rowlisti (combine mu)))
                 (t (setf rowlisti (combine (last params)))))
           (mapcar #'(lambda (x y) (setf (elt (elt collist x) i) y))
                                              (iseq l1) rowlisti)
           (setf yminusmu (- rowlisti mu))
           (setf outprod (outer-product yminusmu yminusmu))
           (setf outprodsum (+ outprodsum outprod))
           (setf newsigma (+ newsigma (first params) outprod))
         )
        (setf loglike (- (* l2 (determinant siginv))
                           (sum (diagonal (matmult siginv outprodsum)))))
        (setf likehist (append likehist (list loglike)))
;;        (if (< loglike (last likehist)) (setf bool t))
;;        (if (not bool) 
        (setf siginv (inverse (/ newsigma l2)))
;;        (if (not bool) 
        (setf mu (map 'vector #'mean collist))
;;       (until bool))
       )
    (setf data (apply #'bind-columns collist))
    (setf results (list data mu (inverse siginv) (butlast likehist)))
    (print-results results)
 )
)
(defun print-results (results)
#|    (format t "The Imputed Data Matrix is: ~%")
    (print-matrix (first results))
|#
    (format t "~% The Estimated Variable means are: ~%")
    (print (second results))
    (format t "~%~% The Estimated Covariance Matrix is: ~%")
    (print-matrix (third results))
    (format t "~%~% The Log-Likelihood History is: ~%")
    (print (elt results 3))
    (format t "~% All of the above are in the variable results ~%")
    (def pp (plot-points (iseq (length (combine (elt results 3)))) 
                                         (combine (elt results 3))))
)
(defun param-est (sig data missind presind mu i torf l1)
"Args: (SIGMA DATA MISSIND PRESIND MU I TORF L1)
Takes in the data, missing and present indicies and estimates the parameters
of a multivariate normal distribution.  SIG, the inverse covariance matrix
should be in the form of a matrix.  Data must be in the form of a list of
vectors MISSIND must be in the form of a list for case i.  PRESIND must also
be in the form of a list for case i.  MU must be in the form of a vector.  I
is the case on which the imputation is being performed.  TORF is a boolean 
that is true if there are at least one missing observation and at least one 
present observation for case i.  l1 is the dimensionality of the multivariate 
normal."
(cond (torf
      (let* (
            (smm-1 (inverse (apply #'bind-rows
                   (mapcar #'(lambda (x) (mapcar #'(lambda (y) (aref sig
                              x y)) missind)) missind))))
            (smo (apply #'bind-rows (mapcar #'(lambda (x) (mapcar
                          #'(lambda (y) (aref sig x y)) presind)) missind)))
            (vi (make-vi smm-1 missind l1))
            )
       (list vi smm-1 smo)
     )
     )
      (t   (cond ((not (not missind)) (list (inverse sig)))
                  (t    (list 0 (mapcar #'(lambda (x) (elt x i)) data)))
           )
      )
)
)

(defun impute (params missind presind mu data i)
"Args: (PARAMS MISSIND PRESIND MU DATA I).
Imputes the missing data vector for case i indexed by missind by
regressing the missing values on the non-missing values.  PARAMS is
a list of matricies and lists.  MISSIND are the indicies of the missing
variables for case i.  PRESIND are the indicies of the present variables
for case i.  MU is the current mean vector.  DATA is a column list of
the data.  This funtion is only called if there are at least one present 
variable and at least one missing variable for case i.  Returns the list 
of imputed and observed values for case i."
 (let* (
       (dati (mapcar #'(lambda (x) (elt x i)) data)) 
       (newvals (+ (select mu missind) (inner-product
                      (matmult (second params) (third params))
                      (- (select dati presind)
                         (select mu presind)))))
       )
  (setf (select dati missind) newvals)
  dati
 )
)
(defun make-vi (smm-1 missind l1)
"Args:  smm-1 missind l1.
Returns a matrix of size l1xl1 with smm-1 in place of the missing variables
and zeros elsewhere."
 (let* (
       (va (make-array (list l1 l1) :initial-element 0))
       (ra (mapcar #'(lambda (x y) (mapcar #'(lambda (z w)
                                  (setf (aref va x z) (aref smm-1 y w)))
                                   missind (iseq (length missind))))
                                   missind (iseq (length missind))))
       )
   va
 )
)