Rainfall prediction is vital for the management of available water resources. Accordingly, this study used large lagged climate indices to predict rainfall in Iran's Sefidrood basin. A radial basis function neural network (RBFNN) and a multilayer perceptron (MLP) network were used to predict monthly rainfall. The models were trained using the naked mole rat (NMR) algorithm, firefly algorithm (FFA), genetic algorithm (GA), and particle swarm optimization (PSO) algorithm. Large lagged climate indices, as well as three hybrid models, i.e., inclusive multiple model (IMM)-MLP, IMM-RBFNN, and the simple average method (SAM), were then employed to predict rainfall. This paper aims to predict rainfall using large climate indices, ensemble models, and optimized artificial neural network models. Also, the paper considers the uncertainty resources in the modeling process. The inputs were selected using a new input selection method, namely a hybrid gamma test (GT). The GT was integrated with the NMR algorithm to create a new test for determining the best input scenario. Therefore, the main innovations of this study were the introduction of the new ensemble and the new hybrid GT, as well as the new MLP and RBFNN models. The introduced ensemble models of the current study are not only useful for rainfall prediction but also can be used to predict other metrological parameters. The uncertainty of the model parameters and input data were also analysed. It was found that the IMM-MLP model reduced the root mean square error (RMSE) of the IMM-RBFNN, SAM, MLP-NMR, RBFNN-NMR, MLP-FFA, RBFNN-FFA, MLP-PSO, RBFNN-PSO, MLP-GA, and RBFNN-GA, MLP, and RBFNN models by 12%, 25%, 31%, 55%, 60%, 62%, 66%, 69%, 70%, 71%, 72%, and 72%, respectively. The IMMs, such as the IMM-MLP, IMM-RBFNN, and SAM, outperformed standalone models. The uncertainty bound of the multiple inclusive models was narrower than that of the standalone MLP and RBFNN models. The MLP-NMR model decreased the RMSE of the RBFNN-NMR, RBFNN-FFA, RBFNN-PSO, and RBFNN models by 15%, 26%, 37%, 42%, and 45%, respectively. The proposed ensemble models were robust tools for combining standalone models to predict hydrological variables. [ABSTRACT FROM AUTHOR]