# demonstrate simulation assessment of an estimator

nsim <- 1000        # number of simulations
nobs <- 30          # number of obs per data set
mu <- 2		  # population mean and s.d. and dl
sigma <- 0.3
dl <- 6

allmean <- allsd <- rep(0, nsim)  # vectors to store estimates

for (i in 1:nsim) {  # loop over simulations
  x <- exp(rnorm(nobs, mu, sigma))   # simulate lognormal
  cens <- (x < dl)   # true if x < dl, otherwise F
  x[cens] <- dl      # replace cens values by the dl

  temp <- ros(x,cens)  # lognormal ROS calculations

  allmean[i] <- mean(temp)  # extract and store mean
  allsd[i] <- sd(temp)      #  and sd
  }

# now have means and sds from 1000 simulated data sets
#  can do all sorts of analysis

plot(allmean,allsd)  # show the correlation

hist(allmean,prob=T)  # show the distribution
lines(density(allmean))

mean(allmean)  # estimate expected value
sd(allmean)    # and s.e.