# creation for random data set
x <- rnorm(100,17,24)
y <- rnorm(100,7,11)
A <- matrix(c(10,25,8,40),nrow=2,ncol=2)

z <- rbind(x, y) # now 2 x 100
w <- A %*% z     # 2 x 100

x <- w[1,]
y <- w[2,]


         print("meanX      meanY       varX       varY       covar") 


         # mean of the variates

         u1 <- mean(x)
	 u2 <- mean(y)
          
         # Variances of the variates
         varX <- var(x)
	 varY <- var(y)
         
         # coVariances of the variates
         covXY <- cov(x,y)
         

         # printing mean , var , covar


         print(c(u1,u2,varX,varY,covXY))



# now replace randomly with NA
x[sample(1:length(x), 10)] <- NA
y[sample(1:length(y), 10)] <- NA

# just a variable
sample <- data.frame("x"=x,"y"=y)


# Vector for missing values
missing <- is.na(sample$x) | is.na(sample$y)
xNA <- array(sample$x[missing], dim=c(1,length(sample$x[missing])))
yNA <- array(sample$y[missing], dim=c(1,length(sample$y[missing])))


# find the  mean of the non missing values in the array

mX <- mean(sample$x,na.rm=TRUE)
mY <- mean(sample$y,na.rm=TRUE)


# Imputing the missing values with current estimated means

 x<-array(ifelse(is.na(sample$x), mX, sample$x),dim=c(1,length(sample$x)))
 y<-array(ifelse(is.na(sample$y), mY, sample$y),dim=c(1,length(sample$y)))


         print("      meanX      meanY       varX       varY       covar") 

# algorithm function to find out new estimate values

algoResult <- function(parrX,parrY,parrXNA,parrYNA,ictr) {

	# Variables for New X and Y
	newx <- array(parrX, dim=c(1,length(parrX)))
	newy <- array(parrY, dim=c(1,length(parrY)))

         # mean of the variates

         u1 <- mean(parrX)
	 u2 <- mean(parrY)
          
         # Variances of the variates
         varX <- var(parrX[1,])
	 varY <- var(parrY[1,])
         
         # coVariances of the variates
         covXY <- cov(parrX[1,],parrY[1,])
         

         # printing mean , var , covar


         print(c(u1,u2,varX,varY,covXY))


          # Expected or updated values for the missing values         

       
          # loop for finding positon of missing vectors
	  for (i in 1:length(missing)) {


	     if (missing[i]==TRUE) 


	        {
          
	         if (is.na(sample$x[i]== TRUE))
		          {
		            newx[i] <- u1 + covXY/varY * (sample$y[i] -u2)

		          } 		
                     if (is.na(sample$y[i]== TRUE))
		          {

		            newy[i]<- u2 + covXY/varX * (sample$x[i] - u1)

 
                          } 

	        }
	     }

                            l <- list(c(newx),c(newy)) 
	                    names(l)<-c("X1","Y1"); as.data.frame(l)          
			    return(l)	       

       }


# the iteration loop to do the algorithm
for (ctr in 1:50) {

   Output <- algoResult(x,y,xNA,yNA,ctr); Output

   x <- array(Output$X1, dim=c(1,length(Output$X1)))
   y <- array(Output$Y1, dim=c(1,length(Output$Y1)))


   }