The beam output of double scattering proton system is difficult to be accurately modeled by treatment planning system (TPS). This study aims to design an empirical method using the analytical and machine learning (ML) models to estimate proton output in a double scattering proton system. Three analytical models and three ML models using Gaussian process regression (GPR) were generated on a training dataset consisting of 1544 clinical measurements, and the accuracy of each model was validated against additional 241 clinical measurements as testing dataset. Two most robust models (polynomial model and the ML GPR model with exponential kernel) were selected, and these two independent models agreed with less than 2% deviation using the testing dataset. The minimum number of samples needed for either model to achieve sufficient accuracy (± 3%) was determined by evaluating the mean average percentage error (MAPE) with increasing sample number, and the differences between the estimated outputs using the two models were also compared for 1000 proton beams with randomly generated range, and modulation for each option. These two models can be used as an independent output prediction tool for a double scattering proton beam, and a secondary output check tool for cross check between themselves.