R/RunningMCMC.R

In EditImputeCont: Simultaneous Edit-Imputation for Continuous Microdata

createModel <- function(data.obj,K){

	if ( class(data.obj) != "EditIn.data" ){
		stop("data.obj should be 'EditIn.data' format generated by 'readData' function")
	}
	
	Y = as.matrix(data.obj$Y.input)

	EIN = min( min(Y[Y>0]), 0.01) ; Y[Y==0] = EIN

	logBound_L_U = log(data.obj$Bound.LU)

  model <- new(bei,Y,data.obj$Edit.matrix,logBound_L_U,data.obj$n.balance,K)    # new(bei,D_obs,Edit,logB, # of balance edits,  K )
	
  # Careful with the no. of balance edits
  
	current_time = Sys.time(); print(paste("Current time,",current_time))
	print("Initializing S and D .....................")
  model$InitializeSandD()

	# current_time = Sys.time(); print(paste("Current time,",current_time))
	print("Building Feaibility Map...................")
	model$BuildFeasibilityMap()

	# current_time = Sys.time(); print(paste("Current time,",current_time))
	print("Simulating log(Uniform)...")
	model$Simulate_logUnif_case2(100) # no. of simulation for calculating Unif_case2
	
	return(model)
	
} # Initialize <- function

multipleEI <- function(model.obj,n.burnin,m,int.btw.EI,show.iter=TRUE){
	
	if ( class(model.obj)[[1]] != "Rcpp_bei" ){
		stop("model.obj needs to be prepared by 'createModel' function")
	}
	
	n_sample = dim(model.obj$D.Initial)[[1]]
	n_var = dim(model.obj$D.Initial)[[2]]
	
	X_EI_Mat = array(0,c(m,n_sample,n_var))
	
	total_iter = n.burnin + m*int.btw.EI
	
  S_small_init = model.obj$S.Compact.Initial
  sum_s_i = apply(S_small_init,1,sum)
  
	if (model.obj$msg.level >= 1){
		print(paste("avg_sum_s_ij of initial S_small= ", round(mean(sum_s_i),2), sep=""))
  }
	
  # burn-in		
	
	if (show.iter==TRUE){
		current_time = Sys.time(); print(paste("Current time,",current_time)) ; prev_time = current_time
	}
	
	print( paste("Total iteration=", (total_iter) ) )
		
	print("Burn-in ..................................")
	
	if ( n.burnin > 20 ){
		
		n_repeat_burnin = floor( n.burnin / 20 )
		resid_n = n.burnin - n_repeat_burnin * 20 
		
		for (i_repeat in 1:n_repeat_burnin){
			model.obj$Run(20)
				current_time = Sys.time()
			if (show.iter==TRUE){
				int_time = round( as.numeric(difftime(current_time, prev_time, units = "mins")), 1 )
				# print(paste("Current time,",current_time))
				print(paste("Iter=",(i_repeat*20),", Burn-in.",int_time," Min for prev. 20 iter",sep=""))
			}
				prev_time = current_time
		}
		
		if (resid_n>0){
			model.obj$Run(resid_n)
		}
		
	} else {
		model.obj$Run(n.burnin)
	}
		
	# Store EI datasets

	print("Drawing Edit-imputed datasets ............")
	
  for ( i_EI in 1:m ){
		
    model.obj$Run(int.btw.EI)
    X_EI_Mat[i_EI,,] = model.obj$Y.edited;
		
			current_time = Sys.time()
			int_time = round( as.numeric(difftime(current_time, prev_time, units = "mins")), 1 )
			# print(paste("Current time,",current_time))
			message(paste("Iter=",(n.burnin+i_EI*int.btw.EI),", EI ", i_EI, "/",m,". ", int_time," Min for prev. ",int.btw.EI," iter",sep=""))
			prev_time = current_time
	
	} # for ( i_EI in 1:no_EI )
		
	if (min(model.obj$D.Initial)>=0.01){
		X_EI_Mat = round(X_EI_Mat,1)
	}
		
	return(X_EI_Mat)
		
} # multipleEI <- function