# --- 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")