The impression that I got, based on several papers, books and articles that I've read, is that the recommended way of fitting a probability distribution on a set of data is by using maximum likelihood estimation (MLE). However, as a physicist, a more intuitive way is to just fit the pdf of the model to the empirical pdf of the data using least squares. Why then is MLE better than least squares in fitting probability distributions? Could someone please point me to a scientific paper/book that answers this question?
My hunch is because MLE does not assume a noise model and the "noise" in the empirical pdf is heteroscedastic and not normal.