# --- MUST source ss1.R and EM1.R 


y=matrix(scan("/mydata/Jones.dat"),ncol=3)
num=nrow(y)       
   # measurement matrices are A = I or 0 matrix
   A=array(0, dim=c(3,3,num))  # creates num 3x3 zero matrices
   for(k in 1:num){
      if (y[k,1]>0) A[,,k]=diag(1,3)  # if obs in 1st column is observed, they're all observed
                }
# initial values
 mu0=matrix(0,3,1) 
 Sigma0=diag(c(.1,.1,1),3)
 Phi=diag(1,3)
 cQ=diag(c(.1,.1,1),3)
 cR=diag(c(.1,.1,1),3)
 
# Ups, Gam, and input are =0 because they're not used
# the iteration number is printed to the screen but
#   you have to scroll the screen manually to see it
(em=EM1(num,y,A,mu0,Sigma0,Phi,0,0,cQ,cR,0,100,.01))


#  graph smoother
 ks=Ksmooth1(num,y,A,em$mu0, em$Sigma0,em$Phi,0,0,chol(em$Q),chol(em$R),0) 
# note that y_t^n = x_t^n for t=1,...,n 
  y1smoo=ks$xs[1,,]
  y2smoo=ks$xs[2,,]
  y3smoo=ks$xs[3,,]
  p1=2*sqrt(ks$Ps[1,1,])
  p2=2*sqrt(ks$Ps[2,2,])
  p3=2*sqrt(ks$Ps[3,3,])
#
  par(mfrow=c(3,1))
  plot.ts(y[,1],type="p", ylim=c(1,5), main="log(white blood count)" )
  lines(y1smoo)
  lines(y1smoo+p1,lty="dashed")
  lines(y1smoo-p1,lty="dashed")
  plot.ts(y[,2], type="p",ylim=c(3,6), main="log(platelet)" )
  lines(y2smoo)
  lines(y2smoo+p2,lty="dashed")
  lines(y2smoo-p2,lty="dashed")
  plot.ts(y[,3], type="p", ylim=c(20,40), main="HCT" )
  lines(y3smoo)
  lines(y3smoo+p3,lty="dashed")
  lines(y3smoo-p3,lty="dashed")