R/ezsim.r

In ezsim: provide an easy to use framework to conduct simulation

#' Create an ezsim object from 4 important arguments(dgp,estimator,true,parameter_def). 
#' @name ezsim
#' @aliases ezsim dgp
#' @title Create an ezsim Object.
#' @param m The number of times you want to simulate.
#' @param estimator A function takes the return value of dgp as argument and return estimates of the dgp
#' @param dgp A function defines the data generating process(dgp). It returns an object (usually it is a \code{data.frame}) which can be used as argument of \code{estimator}. It will be evaluated under an environment generated by a set of parameter generated by parameter_def.
#' @param parameter_def A parameter_def object will be used in this simulation.
#' @param true_value A function defines the true value of estimators(TV). Similar to dgp, but it returns the true value of estimates. It is necessary for computing the bias and rmse of the estimator. The length of its return value must be the same as the lenght of \code{estimator}. If it is missing, True Value, Bias and rmse will be NA in the summary of ezsim.
#' @param display_name Display name for the name of parameter and estimator. see \code{\link{plotmath}} for details.
#' @param estimator_parser A function to parse the return value of estimator function
#' @param auto_save number of auto save during the simulation, default is 0.
#' @param run Whether the simulation will be ran right after the creation.
#' @param run_test Whether to perform a test before the simulation.
#' @param use_seed The seed to be used in the simulation. If \code{use_core=1}, \code{set.seed(use_seed)} will be called. If \code{use_core=>1} and \code{cluster=NULL}, \code{clusterSetRNGStream(cluster,use_seed)} will be used. Ignored if \code{use_core=>1} and \code{cluster} is provided. 
#' @param use_core Number of cpu core to be used.
#' @param cluster_packages Names of the packages to be loaded in the cluster. 
#' @param cluster cluster for parallelization. If it is NULL, a cluster with \code{use_code} cores will be created automatically (will be removed after the simulation) 
#' @return An ezsim object.
#' @author TszKin Julian Chan \email{ctszkin@@gmail.com}
#' @seealso \code{\link{createParDef}} \code{\link{setBanker}},\code{\link{setSelection}} \code{\link{summary.ezsim}}
#' @export
#' @examples         
#' \dontrun{
#' ## Example 1
#' ezsim_basic<-ezsim(
#'     m             = 100,
#'     run           = TRUE,
#'     display_name  = c(mean_hat="hat(mu)",sd_mean_hat="hat(sigma[hat(mu)])"),
#'     parameter_def = createParDef(list(n=seq(20,80,20),mu=c(0,2),sigma=c(1,3,5))),
#'     dgp           = function() rnorm(n,mu,sigma),
#'     estimator     = function(x) c(mean_hat = mean(x), 
#'                                  sd_mean_hat=sd(x)/sqrt(length(x)-1)),
#'     true_value    = function() c(mu, sigma / sqrt(n-1))
#' )
#' 
#' ## Test whether an ezsim object is valid. 
#' ## Print the result of the test and dont return the name of estimator.
#' test(ezsim_basic,print_result=TRUE,return_name=FALSE)
#'  
#' ## Summary of an ezsim object
#' summary(ezsim_basic)
#' 
#' ## Summary of a subset of ezsim object
#' summary(ezsim_basic,subset=list(estimator='mean_hat',n=c(20,40),sigma=c(1,3)))
#' 
#' ## More Summary Statistics
#' summary(ezsim_basic,simple=FALSE,subset=list(estimator='mean_hat',n=c(20,40),sigma=c(1,3)))
#' 
#' ## Customize the Summary Statistics
#' summary(ezsim_basic,stat=c("q25","median","q75"),Q025=quantile(value_of_estimator,0.025), 
#'   Q975=quantile(value_of_estimator,0.975),subset=list(estimator='mean_hat',n=c(20,40),sigma=c(1,3)))
#' 
#' ## Plot an ezsim object
#' plot(ezsim_basic)
#' ## Subet of the Plot
#' plot(ezsim_basic,subset=list(estimator="sd_mean_hat",mu=0))
#' plot(ezsim_basic,subset=list(estimator="mean_hat",sigma=3))
#' ## Parameters Priority of the Plot
#' plot(ezsim_basic,subset=list(estimator="sd_mean_hat",mu=0),parameters_priority=c("sigma","n"))
#' plot(ezsim_basic,subset=list(estimator="mean_hat",sigma=c(1,3)),parameters_priority="mu")
#' 
#' ## Density Plot
#' plot(ezsim_basic,'density')
#' plot(ezsim_basic,"density",subset=list(estimator="mean_hat",sigma=3),parameters_priority="n",
#'   benchmark=dnorm)
#' plot(ezsim_basic,"density",subset=list(estimator="mean_hat",mu=0),parameters_priority="n" ,
#'   benchmark=dnorm)
#' 
#' ## Plot the summary ezsim
#' plot(summary(ezsim_basic,c("q25","q75")))
#' plot(summary(ezsim_basic,c("q25","q75"),subset=list(estimator='mean_hat')))
#' plot(summary(ezsim_basic,c("median"),subset=list(estimator='sd_mean_hat')))
#' 
#' ## Example 2
#' ezsim_ols<-ezsim(
#'     m             = 100,    
#'     run           = TRUE,
#'     display_name  = c(beta_hat='hat(beta)',es='sigma[e]^2',xs='sigma[x]^2',
#'	                      sd_beta_hat='hat(sigma)[hat(beta)]'),
#'     parameter_def = createParDef(selection=list(xs=c(1,3),beta=c(0,2),n=seq(20,80,20),es=c(1,3))),
#'     dgp           = function(){
#'                         x<-rnorm(n,0,xs)
#'                         e<-rnorm(n,0,es)
#'                         y<-beta * x + e
#'                         data.frame(y,x)
#'                     },
#'     estimator     = function(d){
#'                         r<-summary(lm(y~x-1,data=d))
#'                         out<-r$coef[1,1:2]
#'                         names(out)<-c('beta_hat','sd_beta_hat')
#'                         out
#'                     },
#'     true_value    = function() c(beta, es/sqrt(n)/xs) 
#' )
#' summary(ezsim_ols)
#' plot(ezsim_ols)
#' plot(ezsim_ols,subset=list(beta=0))
#' 
#' plot(ezsim_ols,'density')
#' plot(ezsim_ols,'density',subset=list(es=1,xs=1))
#' 
#' 
#' ## example 3
#' ezsim_powerfun<-ezsim(
#'     run           = TRUE,   
#'     m             = 100,
#'     parameter_def = createParDef(selection=list(xs=1,n=50,es=c(1,5),b=seq(-1,1,0.1))),
#'     display_name  = c(b='beta',es='sigma[e]^2',xs='sigma[x]^2'),
#'     dgp           = function(){
#'                         x<-rnorm(n,0,xs)
#'                         e<-rnorm(n,0,es)
#'                         y<-b * x + e
#'                         data.frame(y,x)
#'                     },
#'     estimator     = function(d){
#'                         r<-summary(lm(y~x-1,data=d))
#'                         stat<-r$coef[,1]/r$coef[,2]
#' 
#'                         # test whether b > 0
#'                         # level of significance : 5%
#'                         out <- stat > c(qnorm(.95), qt(0.95,df=r$df[2]))
#'                         names(out)<-c("z-test","t-test")
#'                         out
#'                     }
#' )
#' plot(ezsim_powerfun,'powerfun')
#' }

ezsim <-
function(m,estimator,dgp,parameter_def,true_value=NULL,
		display_name=NULL,estimator_parser=NULL,
		auto_save=0,run=TRUE,run_test=TRUE, use_seed = round(runif(1)*1e+5) , use_core=1,cluster_packages=NULL, cluster=NULL ){
	
	out<-list(m=m,estimator=estimator,true_value=true_value,dgp=dgp,parameter_def=parameter_def,display_name=display_name,estimator_parser=estimator_parser,auto_save=auto_save,use_seed=use_seed,use_core=use_core,cluster_packages=cluster_packages, cluster=cluster )
	class(out)<-"ezsim"
	
	# check estimator_parser
	if (!is.function(out$estimator_parser)){
		out$estimator_parser <- function(x) x 
	}
	
	# set parallel=FALSE if number_of_workers==1 
	# if (out$number_of_workers==1 & out$parallel){
	# 	warning("number_of_workers==1, set parallel=FALSE")
	# 	out$parallel=FALSE
	# }
	
	## Decide whether we need to stop the cluster after the test
	create_cluster_flag <- FALSE

	## create workers for parallel computing
	if (out$use_core > 1 & is.null(out$cluster) & (run | run_test)){
		out$cluster <- makeCluster(out$use_core)
		if (!is.null(out$cluster_packages)){
			for (i in 1:length(out$cluster_packages))
				eval(substitute( 
					clusterEvalQ(out$cluster , require(w, character.only=TRUE)   )  ,
					list( w=out$cluster_packages[i] )
				))
		}
		clusterSetRNGStream(out$cluster,out$use_seed)
		create_cluster_flag<-TRUE
	}

	if (out$use_core == 1) 
		set.seed(out$use_seed)

	## using tryCatch to make sure the cluster is stopped and removed
	tryCatch({
		# test ezsim
		if (run_test){
			tryCatch({
				test(out)
			}, error = function(e){
				cat("\n")
				print(e)
				cat("\n")
				stop("Error in testing ezsim\n")
			})
		}
		# run ezsim
		if (run){
			tryCatch({
				out<-run(out)
			}, error = function(e){
				cat("\n")
				print(e)
				cat("\n")
				stop("Error in running ezsim\n")
			})
		}
	}, finally = function(e){
		if (create_cluster_flag){
			tryCatch({
				stopCluster(out$cluster)
			}, finally = {
				out$cluster<-NULL
			})
		}	
	})
	return(out)
}