Why does $\underset{\theta}{\arg\max }p_{\text{model}}(X;{\theta}) =\underset{\theta}{\arg\max} \prod_{i=1}^{m}p_{\text{model}}({x}^{(i)} ;\theta)$? 
We have seen some definitions of common estimators and analyzed their properties. But where did these estimators come from? Rather than guessing that some function might make a good estimator and then analyzing its bias and variance, we would like to have some principle from which we can derive specific functions that are good estimators for different models.
The most common such principle is the maximum likelihood principle.
Consider a set of $m$ examples $\mathbb{X}=\left\{x^{(1)}, \ldots, x^{(m)}\right\}$ drawn independently from the true but unknown data-generating distribution $p_{\text {data }}(\mathbf{x})$.
Let $p_{\text {model }}(\mathbf{x} ; \boldsymbol{\theta})$ be a parametric family of probability distributions over the same space indexed by $\boldsymbol{\theta}$. In other words, $p_{\text {model }}(x ; \boldsymbol{\theta})$ maps any configuration $\boldsymbol{x}$ to a real number estimating the true probability $p_{\text {data }}(x)$.
The maximum likelihood estimator for $\theta$ is then defined as
$$
\begin{aligned}
\theta_{\mathrm{ML}} &=\underset{\theta}{\arg \max } p_{\text {model }}(\mathbb{X} ; \boldsymbol{\theta}) \\
&=\underset{\theta}{\arg \max } \prod_{i=1}^{m} p_{\text {model }}\left(\boldsymbol{x}^{(i)} ; \boldsymbol{\theta}\right)
\end{aligned}
$$
This product over many probabilities can be inconvenient for various reasons. For example, it is prone to numerical underflow. To obtain a more convenient but equivalent optimization problem, we observe that taking the logarithm of the likelihood does not change its arg max but does conveniently transform a product

I am a slow learner in Mathematics. I don't get the bold part as it is too abstract for me.

*

*I guess "$p_{\text {model }}(x ; \boldsymbol{\theta})$" which "maps any configuration $\boldsymbol{x}$ to a real number estimating the true probability $p_{\text {data }}(x)$" is a function with "$\theta$" parameter (or parameters) that more or less fit the data.


*I struggle with $\underset{\theta}{\arg \max }$. I guess that saying that the maximum likelihood estimator for $\theta$ is then defined as $\theta_{\mathrm{ML}} =\underset{\theta}{\arg \max } p_{\text {model }}(\mathbb{X} ; \boldsymbol{\theta})$ means that the maximum likelihood estimator for $\theta$ is actually looking for the "$\theta$" parameter (or parameters) that find the most likely function that explains observed data.


*I really don't get why $\underset{\theta}{\arg \max } p_{\text {model }}(\mathbb{X} ; \boldsymbol{\theta}) =\underset{\theta}{\arg \max } \prod_{i=1}^{m} p_{\text {model }}\left(\boldsymbol{x}^{(i)} ; \boldsymbol{\theta}\right)$. Why do we need to multiply the progbabilities to find the most likely function that explains observed data?
I took an example by gregmacfarlane from the answer to the question "Maximum Likelihood Estimation (MLE) in layman terms" to understand:

Maximum Likelihood Estimation (MLE) is a technique to find the most likely
function that explains observed data. I think math is necessary, but don't let it
scare you!
Let's say that we have a set of points in the $x,y$ plane, and we want to know
the function parameters $\beta$ and $\sigma$ that most likely fit the data
(in this case we know the function because I specified it to create this
example, but bear with me).
data   <- data.frame(x = runif(200, 1, 10))
data$y <- 0 + beta*data$x + rnorm(200, 0, sigma)
plot(data$x, data$y)


In order to do a MLE, we need to make assumptions about the form of the function.
In a linear model, we assume that the points follow a normal (Gaussian) probability
distribution, with mean $x\beta$ and variance $\sigma^2$:  $y = \mathcal{N}(x\beta, \sigma^2)$.  The equation of this probability density function is: $$\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left(-\frac{(y_i-x_i\beta)^2}{2\sigma^2}\right)}$$
What we want to find is the parameters $\beta$ and $\sigma$ that maximize this
probability for all points $(x_i, y_i)$. This is the "likelihood" function, $\mathcal{L}$
$$\mathcal{L} = \prod_{i=1}^n y_i = \prod_{i=1}^n \dfrac{1}{\sqrt{2\pi\sigma^2}}
                \exp\Big({-\dfrac{(y_i - x_i\beta)^2}{2\sigma^2}}\Big)$$
For various reasons, it's easier to use the log of the likelihood function:
$$\log(\mathcal{L}) = \sum_{i = 1}^n-\frac{n}{2}\log(2\pi) -\frac{n}{2}\log(\sigma^2) -
      \frac{1}{2\sigma^2}(y_i - x_i\beta)^2$$

But I still don't get why we need to compute this multiplication. Did I get wrong the first two points as well?
I am a slow learner in Mathematics and Statistics. Don't hesitate to explain it slowly, or with examples like generic_user's answer to the same question

Since you give a lot of information, I will focus on the question you posted as a title. Bear in mind, that I have problems myself understanding the concepts you bolded and three points you indicate (see my question here).
But let us return to your question:
Why does $\underset{\theta}{\mathrm{arg\,max}}\, p_{model}(\mathbb{X}; \boldsymbol{\theta}) \stackrel{?}{=} \underset{\theta}{\mathrm{arg\,max}} \prod^m_{i=1}p_{model}(\boldsymbol{x}^{(i)}; \boldsymbol{\theta})$?
First, let us ignore the $\underset{\theta}{\mathrm{arg\,max}}$ because it is at both sides, and it is not relevant to the question. Additionally, the argmax affects only $\theta$ so we can ignore it too.
So now the question reduces to:
Why does $p_{model}(\mathbb{X}) \stackrel{?}{=} \prod^m_{i=1}p_{model}(\boldsymbol{x}^{(i)})$?
Let us abstract ourselves from the context of $p_{data}$. What matters is that $p_{data}$ is a probability distribution, so I will call it $p_{dis}$. Now we have:
Why does $p_{dis}(\mathbb{X}) \stackrel{?}{=} \prod^m_{i=1}p_{dis}(\boldsymbol{x}^{(i)})$?
Our goal is to reach the right part, starting from the left part. So, $\mathbb{X} = \{\boldsymbol{x}^{(i)}, \cdots, \boldsymbol{x}^{(m)}\}$ per definition. Substituting in the left part, we have:
$p_{dis}(\{\boldsymbol{x}^{(i)}, \cdots, \boldsymbol{x}^{(m)}\})$.
Observe that this is a joint probability of all examples, i.e.:
$p_{dis}(\{\boldsymbol{x}^{(i)}, \cdots, \boldsymbol{x}^{(m)}\}) = p_{dis}(\boldsymbol{x}^{(1)} \cap \boldsymbol{x}^{(2)} \cap \cdots \cap \boldsymbol{x}^{(m)})$.
Now comes the part with the independence assumption:
Per definition (see here for example) "Two events A and B are independent if and only if their joint probability equals the product of their probabilities", i.e., $P(A \cap B) = P(A)P(B)$. Because all exampes $x^{(i)}$ are mutually independent (the (in)famous i.i.d.), we have:
$p_{dis}(\{\boldsymbol{x}^{(i)}, \cdots, \boldsymbol{x}^{(m)}\}) = p_{dis}(\boldsymbol{x}^{(1)}) p_{dis}(\boldsymbol{x}^{(2)}) \cdots p_{dis}(\boldsymbol{x}^{(m)})$, which is a product over the probability of every example:
$p_{dis}(\{\boldsymbol{x}^{(i)}, \cdots, \boldsymbol{x}^{(m)}\}) = \prod^m_{i=1} p_{dis}(\boldsymbol{x}^{(i)})$ and so we arrived at the right part.
If we like, we can go backwards:

*

*Go back to $p_{model}$: $p_{model}(\{\boldsymbol{x}^{(i)}, \cdots,
    \boldsymbol{x}^{(m)}\}) = \prod^m_{i=1}
    p_{model}(\boldsymbol{x}^{(i)})$


*Go back to the set of observed examples $\mathbb{X}$:
$p_{model}(\mathbb{X}) = \prod^m_{i=1}
p_{model}(\boldsymbol{x}^{(i)})$


*Incorporate argmax and the only parameter affected by it ($\theta$):
$\underset{\theta}{\mathrm{arg\,max}}\, p_{model}(\mathbb{X};
\boldsymbol{\theta}) = \underset{\theta}{\mathrm{arg\,max}}
\prod^m_{i=1}p_{model}(\boldsymbol{x}^{(i)}; \boldsymbol{\theta})$
That being said, we considered that $p_{model}$ was a probability distribution. But it seems to me that $p_{model}$ is nothing else but the likelihood function (LF) and if so, it does not have to integrate to one, therefore must not be considered as a probability distribution (see above my mentioned question). Nevertheless, I hope that the independence notion holds as well for general distributions like $p_{model}$, which is, according to the authors, "a parametric family of probability distributions".
