An intelligent and accurate modeling method is proposed combining near-infrared (NIR) spectroscopy for measuring organic matter content in soil samples. The proposed method uses Monte Carlo (MC) random sampling in the training set, where subsets were randomly selected from the samples and further selected using the butterfly optimization algorithm (BOA) to construct partial least squares (PLS) sub-models, named MC-BOA-PLS. Ultimately, the final prediction was obtained by averaging the predictions of these sub-models. The parameters of MC-BOA-PLS model were optimized, including the iteration number of BOA, the number of butterflies, and the number of PLS sub-models. Results show that MC-BOA-PLS exhibited superior predictive performance to predict organic matter content in soil compared with PLS and BOA-PLS.