Why does Maximum Likelihood estimation maximize probability density instead of probability I am trying to understand Maximum likelihood estimation but it looks like I am missing something rather elementary.


suppose we have an iid random sample $X_1, X_2,..., X_n$ for which the probability density function of each $X_i$ is $f(x_i; \theta)$ where $\theta$ is an unknown parameter. Then, the joint probability density function of $X_1, X_2,..., X_n$ is given by:
$f(X_1=x_1,X_2=x_2,\cdots,X_n=x_n)=\prod_{i=1}^{n} f(x_i, \theta)$


In Maximum Likelihood estimation, we try to maximize $f$ as a function of $\theta$.
Question:
Why do we maximize probability density instead of probability? In what way does it make sense?
Edit: As explained in the answers, the probability of choosing a finite number of sample points from a continuous probability distribution is zero, so maximizing probability doesn't make sense, but how does maximizing probability density make sense? 
 A: $f(x_i, \theta)$ may not be a probability, it is a density function. In general statistics, we don't want to have to make special exceptions for continuous versus discrete random variables all the time, especially since there is a field of mathematics that gives us a unified approach yet allows us to be rigorous about such things.
The rationale for maximizing the product of the densities of a sample, or the likelihood, is much like the rationale for an integral in calculus. Take height, it is a continuous value. And suppose I have some belief about a "normal, maximum entropy Gaussian" spread to underlie this distribution in a population, and it is parametrized by a mean and standard deviation. My height is measured with error, and even if I knew it to an atomic level I could never actually find a probability associated with that single value. The probability that my height is between 5'10" and 5'11" is small, but between 5'10.25" and 5'10.75" is even smaller, and if I squeeze and squeeze this range into an $\epsilon$-ball, the associated probability goes to 0, even if my height happens to be the mean, mode, and median of the population sample. So how is it that this value which is highly characteristic of the population shows such a small probability? A zen answer might be: the infinitessimal differences make up the whole. By look at the density, or the differential of probability, you actually find that a random observation achieving a mean, mode, median is actually very characteristic: it achieves the highest likelihood of any other value in that density.
A: Your question applies only to continuous random variables. In the case of discrete random variables you do use probabilities and not densities. For a continuous random variable, the probability of each point (one value of the variable) is 0, and only intervals have positive probabilities obtained by integrating the density function over the interval. Since the sample consists of points, you cannot multiply probabilities (the result will always be 0) and you must multiply densities (which are in some sense a "representative" of the probability but cannot be called probability).
To be even more specific: "probability density" and "density" are one and the same - two names for the same function. To understand what the density function means you should have a knowledge of calculus. The density function f(x) can be explained as the "slope" of the probability at point x. f(x)dx can be explained as the probability of the point x, which on one hand is equal to 0 (because dx is equal to 0), but on the other hand becomes greater than 0 when integrated over an interval. So f(x) only represents how "dense" the probability is at point x, but is not the probability, still can be used as a "proxy" to probability.
A: I read the question as: why do we start from the density function $f(\boldsymbol{x}|\theta)$ (with $\theta$ constant) to change point of view and interpret it as a function of $\theta$ (with $\boldsymbol{x}$'s constant) that we want to maximize?
Intuitively and absolutely not rigorously, if we consider an infinitesimal interval $d\boldsymbol{x}$ around $\boldsymbol{x}$, then $f(\boldsymbol{x}|\theta)d\boldsymbol{x}$ can be thought as the infinitesimal probability of getting inside that infinitesimal interval, so in a sense it is a probability (i.e. when "summed" over all possible infinitesimal intervals it yields 1 as you would expect from a probability. This summation is called integration in calculus).
Now you want to maximize against $\theta$, so you want to find that value $\hat{\theta}$ such that:
$$\forall \theta: f(\boldsymbol{x}|\hat{\theta})d\boldsymbol{x}\geq f(\boldsymbol{x}|\theta)d\boldsymbol{x}$$
Now... assuming we trust that we can divide by $d\boldsymbol{x}$ on both sides, we obtain:
$$\forall \theta: f(\boldsymbol{x}|\hat{\theta})\geq f(\boldsymbol{x}|\theta)$$
i.e. $\hat{\theta}$ is the value of $\theta$ that maximizes $f(\boldsymbol{x}|\theta)$.
Again, this is not rigorous but I hope it gives you the gist of it. If these "infinitesimals" disturb you, try to think in terms of finite probabilities of falling inside finite intervals $\Delta\boldsymbol{x}$ and then evaluate the limit for this interval's amplitude that goes to 0 in all different $\boldsymbol{x}$'s...
A: The key idea here is to consider that although Point probability is not defined for the continuous probability distribution but we can easily see that probability that the random variable('X') is "around" x is equals to f(X=x)dx. Therefore when its multiplied for all of the points the Likelihood function would not be affected by all of these dx's so that we could ignore those and just maximize the product of f(X=x_i) for all data points.
Hope that helps, Cheers!
Reference: https://ocw.mit.edu/courses/mathematics/18-05-introduction-to-probability-and-statistics-spring-2014/readings/MIT18_05S14_Reading10b.pdf
