dat <- data.frame(Response=c(rlnorm(10, 2.9, 0.3), rlnorm(14, 2.88, 0.38), rlnorm(19, 2.44, 0.08)),TestNum=rep(c(1,4,9), times=c(10,14,19)))
dat$TestNum<-factor(dat$TestNum)
dat

dat_fit1 <- with(dat,
        by(dat[,1], TestNum, fitdist, "lnorm"))
dat_fit2 <-t(sapply(dat_fit1, coef))

draw_dat <- with(dat,
        by(dat[,1], TestNum, function(x){
           out<-replicate(1000, rlnorm(150, meanlog=dat_fit2[,1], sdlog=dat_fit2[,2]), simplify=FALSE)}))

perce0.5<-rapply(draw_dat, quantile, probs=0.5, classes = "numeric", how = "replace")
perce0.5<-do.call(rbind, perce0.5)
perce0.5<-matrix(unlist(perce0.5), nrow(perce0.5), dimnames = dimnames(perce0.5))

mm0.5 <-rowMeans(perce0.5)
mlu0.5<- apply(perce0.5, 1, quantile, probs = c(0.025, 0.975))
mlu<-t(rbind(mm0.5, mlu0.5))
mlu

mm0.5     2.5%    97.5%
1 15.61274 14.25022 17.06212
4 15.51227 14.15273 16.95342
9 15.51304 14.17514 16.96741

for(i in 1:n) plot(dat_fit1[[2]])