The air pollution problem has now amassed worldwide attention due to its multifaceted harm to human health. Exploring the concentration of air pollution and improving forecast have important consideration worldwide. In this research, we analyze the air pollution concentration of Southern Thailand and compare it with the central region. Also, we proposed a methodology based on the lag-dependent Gaussian process (LDGP), a Bayesian non-parametric machine learning model, with a stable optimization approach, which is a cluster-based multi-starter technique based on the Nelder-Mead optimizer. This model also provides the confidence band for forecasted values. We also used autoregressive deep neural network (AR-DNN), autoregressive random forest (AR-RF), gradient boosting (GB), and K-nearest neighbors (KNN) models. A comparison of the proposed methodology was performed on the daily air pollution data collected from the southern provinces and also from the capital of Thailand from 1 January 2018 to 31 December 2022. We used well-established performance evaluation measures to compare the performance of the models. To evaluate the bias due to overfit, we performed a tenfold cross-validation for all the pollutants in each region and compared the models to choose the best one. Moreover, we explored the concentration of air pollution in these regions. Results of descriptive analysis revealed that Bangkok had a much higher concentration of air pollution as compared to the southern region. However, the southern region had higher exposure to PM air pollutants as per WHO recommendations and also had higher exposure to O3 and CO levels. The proposed LDGP model outperformed the other machine learning models for forecasting all air pollutants. Hence, it is recommended to be used by experts for further research and studies with different kernel functions. This research is also expected to contribute to local government planning and prevention and worldwide use of the same methodology for the sustainability of public health.