Linear and nonlinear dynamic analyses become popular in seismic design and seismic performance evaluation procedures. The accuracy of the results of such analyses depends not only on accurate analytic model of structures but also a proper selection of input ground motions. The purpose of this study is to develop a computationally efficient and accurate algorithm for selecting ground motions considering the mean and variance of target response spectrum. The proposed algorithm adopts a simulation based process for generating a desired number ($=n_g$) of simulated response spectra using the Monte Carlo Simulation with information of the mean, variance, and correlation structure of target response spectrum. The accuracy of selecting ground motions using the proposed algorithm is improved by generating multiple (=100) sets of simulated response spectra. Ten sets are chosen from the 100 sets, which best match the mean and standard deviation of a target response spectrum. Ground motions are then selected for each of the top ten sets. The set of ground motions having the mean and standard deviation of their response spectra that best match the target response spectrum is the one to be selected. To verify the accuracy and efficiency of the algorithm, numerical examples are provided.