Understanding maximum likelihood estimation I am told

the method of maximum likelihood says we should use the model that
assigns the greatest probability to the data we have observed;
formally, the maximum likelihood estimator is found by solving 
$\hat{\theta}= \arg_{\theta}\max\{p(x|\theta)\}$  where $p(x|\theta)$
is called the likelihood function.

Am I correct reading $p(x|\theta)$ in English as "the probability of the data given the parameters"?
I am confused because we at first seem to be told that MLE is about the probability of the parameters given the data.
[Update]
I am still confused about whether the likelihood function returns probability or probability density. Because Wikipedia says

The likelihood function (often simply called the likelihood)
represents the probability of random variable realizations conditional
on particular values of the statistical parameters

I am a programmer. When I write a function in code I expect it to return a value of a known type. I want to understand the type that the likelihood function returns.
If the type can be either probability or probability density, why does Wikipedia not make that clear at the start?
 A: Two years too late but here's my go at it :)
DISCLAIMER: this follows Bishop's amazing book Pattern Recognition and Machine Learning. Some sentences are directly copied (particularly in the deterministic approach section). I have not used quotations because it could break up the flow of the words, as I often make additional statements in between. I'm not sure if that's not good practice, so please let me know if that is not the way to do it.
Setup
We observe a real-valued input variable $x$ and we wish to use this scalar observation to predict the value of a real-valued target variable $t$. To this end, we have gathered a set of $N$ realizations $x_i$ of $x$ together with a corresponding set of realizations $t_i$ of $t$
\begin{equation}
\mathbf{x} = (x_1, x_2,...x_N)^T \; , \; \mathbf{t} = (t_1, t_2,...t_N)^T
\end{equation}
In this synthetic curve fitting example, the data is generated by the function $\text{sin}(2 \pi x)$, with random noise included in the target values so we can write
\begin{equation}
y = \text{sin}(2 \pi x) + \mathcal{N}(0, \beta^{-1}).
\end{equation}
Our goal is to exploit the training set in order to make predictions of the value $\hat{t}$ of the target variable for some new value $\hat{x}$ of the input variable. This involves implicitly trying to find the underlying function which we know to be $\text{sin}(2 \pi x)$.
Deterministic Approach: Least Squares
Model Setup
In the deterministic--or least squares--approach, we simply consider our output $t$ to be a parameterized function of the input variable $x$ and weights $\mathbf{w}$. For this example, we will consider a hypothesis set of all polynomial functions of the following form:
\begin{equation}
y(x, \mathbf{w}) = w_0 + w_1x + w_2x^2 + \dots + w_Mx^M = \sum_{j=0}^M w_jx^j \tag{1}
\end{equation}
where $M$ is the order of the polynomial, $x^j$ denotes $x$ raised to the power of $j$, and the polynomial coefficients $w_0,...,w_M$ are collectively denoted by the vector $\mathbf{w}$. Although $y(x, \mathbf{w})$ is a nonlinear function of $x$, it is a linear function of the coefficients $\mathbf{w}$. Functions that are linear in the unknown weights have special properties and are called linear models.
Inference
We want to find $\mathbf{w}$ so that $(1)$ fits the data. To do so, we find a $\mathbf{w}$ that minimizes some notion of "error". The fit to the data is constrained by the family of functions we've decided to consider: polynomial functions of order $M$.
Our notion of error can be mathematically expressed as an error function that measures the misfit between the function $y(x, \mathbf{w})$ for a given value of $\mathbf{w}$, and the training data points. A widely used --and simple-- error function is given by the sum-of-squares of the errors between predictions $y(x_n, \mathbf{w})$ and the corresponding target values $t_n$, so that we minimize
\begin{equation}
E(\mathbf{w}) = \frac{1}{2} \sum_{n=1}^N \left( y(x_n, \mathbf{w})- t_n \right)^2 \tag{2}
\end{equation}
where $\frac{1}{2}$ is for convenience later on. A Visual representation of the sum-of-squares error function in $(2)$ is shown here:

The sum-of-squares error function in $(2)$ is computed by taking one half the sum of the squared distances of each data point from the function $y(x, \mathbf{w})$. These displacements are shown in red.
Solving for $\mathbf{w}$ in this setting is fairly straightforward. Because the error function $(2)$ is a quadratic function of the coefficients $\mathbf{w}$, its derivative with respect to $\mathbf{w}$ will be linear in the elements of $\mathbf{w}$. Therefore, the minimization of the error function has a unique solution, denoted $\mathbf{w}^{\star}$ which can be found in closed form. The resulting polynomial function is then $y(x, \mathbf{w}^{\star})$.
Maximum Likelihood Approach
In the Least Squares approach we viewed curve fitting purely as an error minimization problem. Now we take a first step in viewing it from a probabilistic perspective so that we can express uncertainty in our predictions.
Model Setup
As in the least-squares approach we could consider $x$ to be related to $t$ through a parameterized function $y(x, \mathbf{w})$. However, we may not want to make such a strong statement as saying $t$ is exactly equal to $y(x, \mathbf{w})$. This could be because we think there is noise in the observations $\mathbf{t}$, for example due to measurement error. To introduce such uncertainty, we need to place a distribution over the target variable $t$. A sensible distributional assumption is to place a Gaussian distribution over $t$, with its mean given by the parameterized function $y(x, \mathbf{w})$, and its variance being fixed and unknown. This is visualized here:

Illustration of a Gaussian conditional distribution over $t$ given $x$, where the mean is given by some function of $x$, and the variance is fixed.
A more general way to think about this is through a data-generative process. In the least-squares approach, we related $t$ and $x$ through the equation $t = y(x, \mathbf{w})$. We can instead think about $t$ as being generated by some process involving $x$ and parameters $\mathbf{w}$, where the process is simply the deterministic function $y(x, \mathbf{w})$.
In this section, we extend the process by making the assumption that each $t$ is the result of $y(x, \mathbf{w})$ and some additive uncertainty. This additive uncertainty takes the form of a zero-mean Gaussian distribution with unknown variance. This is to say that we are assuming, to have gotten a particular instance of $t$:
▪ the world took an instance of $x$. 
▪ the world used this instance of $x$ to get the output of the parameterized function $y(x, \mathbf{w})$. 
▪ the world sampled from a zero-mean Gaussian, with fixed and unknown variance, and added it to the output.
This leads to the following model
\begin{equation}
p(t | x, \mathbf{w}, \beta) = y(x, \mathbf{w}) + \mathcal{N}(0, \beta^{-1}) = \mathcal{N}\left( y(x, \mathbf{w}), \beta^{-1} \right) \tag{4}
\end{equation}
where we've used the scaling property of the Gaussian distribution's mean, and we've defined $\beta = \frac{1}{\sigma^2}$. $(4)$ is referred to as the Gaussian noise model.
Inference
As the name suggests, in order to use the training data $( \mathbf{x}, \mathbf{t})$ to determine the values of the unknown parameters $\mathbf{w}$ and $\beta^{-1}$, we will search for a setting of $\mathbf{w}$ that maximizes the likelihood of the data $t$. In other words, we've defined a generative process, and we want to find the correct setting of the parameters $\mathbf{w}$ so that the likelihood of our process having created our observed $\mathbf{t}$ from our observed $\mathbf{x}$ is maximized.
Assuming the data $(\mathbf{x}, \mathbf{t})$ were independently sampled from $(4)$, the likelihood function is simply the product of each conditional distribution and is evaluted for a particular setting of $\mathbf{w}$:
\begin{equation}
p(\mathbf{t}|\mathbf{x}, \mathbf{w}, \beta) = \prod_{n=1}^N \mathcal{N}(t_n|y(x_n, \mathbf{w}), \beta^{-1}) \tag{5}
\end{equation}
Each time we choose a setting for $\mathbf{w}$ and plug it into our model, we are defining a conditional distribution--in particular the one in $(4)$. This conditional distribution may agree with the data we have, or it may not. Examples of agreement and disagreement are shown in Figure 3.

A Gaussian noise model shown for a handful of $x$, with two different settings for $\mathbf{w}$ and $\beta$. On the left is a setting of $\left( \mathbf{w}, \beta \right)$ that yields a model that disagrees with our observed data. On the right is a setting of $\left( \mathbf{w}, \beta \right)$ that yields a model that agrees much better with our observed data. Maximum likelihood looks for the setting of $\left( \mathbf{w}, \beta\right)$ that best agrees with our observed data.
We begin by finding the maximum likelihood estimates for $\mathbf{w}$. For this example, this amounts to taking the derivative of the likelihood $(5)$, setting it equal to zero, and then solving for $\mathbf{w}$. So again, finding the maximum likelihood setting for $\mathbf{w}$ can be found in closed form. It is common to instead maximize the log likelihood instead of $(5)$ for numerical stability and convenience, This which can be written as
\begin{equation}
\text{ln} p(\mathbf{t}|\mathbf{x}, \mathbf{w}, \beta^{-1}) = -\frac{\beta}{2} \sum_{n=1}^N \left(y(x_n, \mathbf{w}) - t_n \right)^2 + \frac{N}{2} \text{ln} \beta - \frac{N}{2} \text{ln}(2 \pi) \tag{6}
\end{equation}
For the purpose of taking the derivative of $(6)$ with respect to $\mathbf{w}$, we can omit the last two terms as they do not depend on $\mathbf{w}$. We can also replace the coefficient $\frac{\beta}{2}$ with $\frac{1}{2}$ since scaling $(6)$ by a constant won't chance the location of the maximum with respect to $\mathbf{w}$. Lastly, we can equivalently minimize the negative log likelihood. This leaves us with minimizing the following:
\begin{equation}
\frac{1}{2} \sum_{n=1}^N \left(y(x_n, \mathbf{w}) - t_n \right)^2 \tag{7}
\end{equation}
and so we see that the sum-of-squares error function has arisen as a consequence of maximizing the likelihood under the assumption of a Gaussian noise distribution.
Once we've found the maximum likelihood estimate for $\mathbf{w}$, which we will denote $\mathbf{w}_{\text{ML}}$, we can use it to find the setting for the precision parameter $\beta$ of the Gaussian conditional distribution. Maximizing $(6)$ with respect to $\beta$ gives
\begin{equation}
\frac{1}{\beta_{\text{ML}}} = \frac{1}{N} \sum_{n=1}^N \left(y(x_n, \mathbf{w}_{\text{ML}}) - t_n \right)^2
\end{equation}
and so we see that the maximum likelihood procedure yields a variance $\sigma^2$ being the average squared deviation between the observed data points and the fitted $y(x, \mathbf{w}_{\text{ML}})$.
A: Parameters are not random variables but fixed unknowns (at least in the likelihood approach to inference) calibrating the distribution of the observables/observations. Outside a Bayesian setup, it is thus incorrect to talk of a probability distribution on the parameters.
The MLE is the numerical value of the parameters that makes the actual observations the most likely to have occurred
$$p(x|\hat\theta) = \max_\theta p(x|\theta)$$ The quantity $p(x|θ)$ in the above thus reads as

*

*the probability of observing the realisation $x$ of the random variable $X$ when the parameter value (indexing its distribution) is equal to $θ$ in a discrete setting and

*the density of $X$ at the value $x$ when the parameter value is equal to $θ$ in a continuous setting.

(I would even avoid the term given since this could be interpreted as a conditional probability or density, which does not make sense if $\theta$ is not a random variable, i.e., outside a Bayesian framework.)
To quote the very originator of the notion of likelihood, R.A. Fisher in ‘On the mathematical foundations of theoretical statistics’ (1922):

I suggest that we may speak without confusion of the likelihood of one value of p being thrice the likelihood of another (…) likelihood is not here used loosely as a synonym of probability, but simply to express the relative frequencies with which such values of the hypothetical quantity p would in fact yield the observed sample.

Note the stress made in the discussed Wikipedia page

The fact that the likelihood function can be defined in a way that includes contributions that are not commensurate (the density and the probability mass) arises from the way in which the likelihood function is defined up to a constant of proportionality, where this "constant" can change with the observation $x$, but not with the parameter $\theta$.

which reinforces the point that the likelihood function is not a probability density function (for its argument $\theta$) and can take numerical values above $1$ (or any other bound).
