bc.ci<-function(res, est, alpha){ ## function to obtain BC CI
  iB<-length(res)
  alpha<-c((1-alpha)/2, (1+alpha)/2)
  z0<-qnorm(sum(res<est)/iB)
  zalpha<-qnorm(alpha)
  bcalpha<-pnorm(2*z0+zalpha)
  BC<-quantile(res,bcalpha)
  c(est, sd(res), BC)
}

boot.mplus<-function(B, datafile, runmplus, estfile, alpha=.95){
  temp.data<-read.table(datafile)
  N<-dim(temp.data)[1]
  write.table(temp.data, 'data.txt', row.names=F, col.names=F)
  system(runmplus,show.output.on.console=F)
  res<-scan(estfile)

  boot.res<-NULL
  for (j in 1:B){
    index<-sample(1:N, replace=T)
    dset<-temp.data[index, ]
    write.table(dset, 'data.txt', row.names=F, col.names=F)
    system(runmplus,show.output.on.console=F)
    temp<-scan(estfile)
    boot.res<-rbind(boot.res, temp)
    cat(j, '\n')
    }

  res<-res[c(19:21,24,25,27,31,33,36)]
  res<-c(res, res[4]*res[5])
  boot.res<-boot.res[,c(19:21,24,25,27,31,33,36)]
  boot.res<-cbind(boot.res, boot.res[,4]*boot.res[,5])
  BC.matrix<-NULL
  for (i in 1:10){
    BC.matrix<-rbind(BC.matrix, bc.ci(boot.res[,i], res[i], alpha))  
    }
  colnames(BC.matrix)<-c('Estimate','S.E.', '2.5%', '97.5%')
  rownames(BC.matrix)<-c('iM','iY','iX','a','b',"c'",'eM2','eY2','eX2','ab')
  list(res=res, boot.res=boot.res, BC=BC.matrix)
}

runmplus<-"c:/programs/mplus/mplus MODEL.ACTIVE.AUX.inp"
allres<-boot.mplus(1000,'active.txt',runmplus,'est.txt', .95)
allres