The absence of key input parameters and use of least squares technique have resulted in most of the available empirical ground motion prediction models (GMPEs) in India, as a function of epicentral distances and magnitudes, resulting in large uncertainty in the predictive model. Here we apply a supervised Machine Learning technique (XGBoost) to design a robust improved GMPE of peak ground acceleration (PGA) for the Indian sub-continent and Nepal, utilizing independent input parameters viz., moment magnitudes, hypocentral distances, focal depths and site proxy (in terms of average seismic shear-wave velocity from the surface to a depth of 30m (Vs30)). This has been achieved using PGA dataset of nine earthquakes of Mw5.6-7.9 from Nepal Himalaya, Uttarakhand Himalaya, Kachchh (Gujarat), and Koyna (Maharashtra). Our dataset also includes 225 strong-motion records of 30 significant Bhuj aftershocks of Mw3.3-5.0 (during 2002–2008) with epicentral distances ranging from 1.0 to 288 km. We also test the ground motion predictability of our supervised ML model using the whole dataset and randomized test dataset (30% of total data points). The correlation coefficient of observed and predicted PGA values for the whole dataset is found to be ±0.994 while our predictability on the randomized test dataset suggests a correlation coefficient of 0.944. Also, the model predictability of our supervised ML model suggests a good prediction of the observed PGA dataset for twelve Indian and Nepal earthquakes of Mw5.6-7.9 (viz., the 2001 Bhuj (Gujarat) earthquake sequence, the 2015 Nepal earthquake sequence, the 1999 Chamoli earthquake, the 1967 Koyna earthquake, the 2009 Mw6.2 Bhutan event and two Mw5.9 Myanmar-India border earthquakes of September 2009 and January 2013), suggesting that our ML GMPE model would work for earthquakes occurring in India and Nepal, within the range of Mw [3.3–7.9], focal depths [3–94 km], epicentral distances [1-565 km] and Vs30 [177–760 m/s].