The shape of earthquake source spectra, traditionally fit by physics-based models, contains important parameters to constrain rupture dimension, duration, and geometry. Here we apply machine learning (ML) to derive single-variable and double-variable data-driven models of source spectra from 3675 Mw>5.5 global earthquakes, assuming that the Fourier transform of source time functions well represent earthquake source spectra below 1 Hz. The single-variable ML model, in the same degree of freedom as the Brune model, improves the goodness of fit by 8.5%. Specifically, the ML model fits the data without systematic bias, whereas the Brune model tends to underestimate at intermediate frequencies and overestimate at high frequencies. The latter discrepancy cannot be modelled by increasing the fall-off exponent in the Brune-type or the Boatwright-type models. The double-variable ML model is compared to existing double-corner-frequency models and is found to capture the second-order features such as the subtle curvature differences around the corner. Our results demonstrate that unsupervised machine learning can extract hidden global characteristics of high-dimensional data and provide observational evidence to amend existing physical models.