How is probability distribution type decided before applying maximum likelihood estimation? What I can't grasp for a long time is how is the assumption about the probability distribution type made? Why do we assume that the human height is distributed exactly normally? Why not any other function that is numerically/shape-wise close to Gaussian? What is the reasoning behind it? I did not manage to find any criteria about how the probability function is selected. Why is function F over function G in a particular case? Do we just essentially guessume it?
 A: To choose a distribution, we must carefully consider the nature of the data being sampled. Some important criteria include (but are not limited to): 


*

*Discrete vs. continuous

*Non-negativity

*Symmetry

*Censoring and truncation.


Based on this, we can narrow down our choice for a reasonable distribution function. There are various tests that can be implemented to check of the choice of the distribution (e.g. $\chi^2$ goodness of fit test). 
In your height example, let's look down the list. Clearly height is a continuous variable, so we can rule out discrete distributions as our choices (Poisson, negative binomial, etc.). Height is a non-negative variable, but, I would guess, that 0 is sufficiently far from the mean (in terms of number of standard deviations) that fitting a normal distribution is not problematic. Whoever choose the normal distribution must have looked at a histogram of the data, and decided the distribution was very close to symmetric (there are tests for this as well). Finally, censoring or truncation is not an issue in height data, so we can disregard criteria 4. 
A: In some cases, the Gaussian assumption can be justified by the central limit theorem. For instance, if the data consist of averages of heights of students in each school, we can treat each average follows a Gaussian distribution. 
In general, the Gaussian assumption can be too restrictive so people tend to use the sub-Gaussian assumption which includes all distributions with fast-decaying tails including Gaussians, all bounded distributions and so on. In this case, you cannot use the maximum likelihood estimation since we do not specify any density function. But still we can prove many commonly used estimators such as sample mean, median, quantiles behave in some optimal manner. 
