1) If n rv's (not one) are independent, then their joint density is the product of their individual densities (this extends the notion of the probability of intersection of n independent events in the case of rv's). If, moreover, they are (assumed) identically distributed, then
$$ L(X_1,\ldots,X_n| \theta)=\prod_{i=1}^nf(X_i|\theta)=f(X_1|\theta)\ldots f(X_n|\theta)$$
is the joint density.
If we view this expression as a function of $\theta$, then it becomes
$$ L(\theta|x_1,\ldots,x_n)=\prod_{i=1}^nf(\theta|x_i)=f(\theta|x_1)\ldots f(\theta|x_n)$$
and now it is the likelihood function of the sample. Is there a difference? Big one: In the joint density, the argument/variable of the function is the X's while the $\theta$'s are treated as constants. In the likelihood function, the arguments/variables are the $\theta$'s while the x's are treated as constants (changing from uppercase to lowercase for the x's is a usual -and good- mnemonic) So the likelihood function is not the joint density - they may have the same functional form, but they are functions of different variables.
2) In view of 1), we do not select the $\theta$'s "so as the rv's assume a particular value...etc" as you write. We select the $\theta$'s given the specific values that the r.v's have already taken (values that constitute our sample), so as the resulting value of the function is maximum. Then we flip our point of view again and say (very loosely speaking), that this value of $\theta$ "maximizes the probability that the sample comes from the assumed distribution".
3) The uniform distribution is special in that the variable itself does not appear in the density function. The specific case we study is even more special, because not even the limits of the interval into which the rv ranges appear in it: the density is equal to unity for the specified range, and zero elsewhere (even though the range of the rv is not [0,1]). Instead of writing a piece-wise density we write compactly $$ f_U (X_i)=I_{[\theta-\frac{1}{2},\;\theta+\frac{1}{2}]} (X_i)$$ Then the joint density for the case of n i.i.d rv's is $$ L(X_1,\ldots,X_n| \theta)=\prod_{i=1}^nf(X_i|\theta)=\prod_{i=1}^n I_{[\theta-\frac{1}{2},\;\theta+\frac{1}{2}]} (X_i) =\min\{I_{[\theta-\frac{1}{2}\le X_1\le\theta+\frac{1}{2}]} ,...,I_{[\theta-\frac{1}{2}\le X_n\le\theta+\frac{1}{2}]}\}$$ where the last equality comes from basic properties of indicator functions.
Consider now the joint density as a likelihood function of $\theta$, with given x's. Then we should maximize this likelihood function? No. We are applying here maximum likelihood estimation, and the ML estimators are M-estimators, which maximize the expected value of the likelihood function (or its sample analogue, whatever is feasible). So our estimation problem here can be stated as
$$ \max_{\theta}E\min\{I_{[\theta-\frac{1}{2}\le x_1\le\theta+\frac{1}{2}]} ,...,I_{[\theta-\frac{1}{2}\le x_n\le\theta+\frac{1}{2}]}\}$$
Let's start the optimization procedure. Note that if we choose some $\hat \theta$ such that some of the realized values in the sample fall outside the specified range, the min{} expression will be equal to zero. But if the chosen $\hat \theta$ has a value that leaves all realized values of the rv's inside the resulting common interval, then the min{} expression will be equal to unity. I adopt the notation $x_1=x_{min}$ and $x_n=x_{max}$. So a first conclusion is that we must have $$ x_1\ge \theta-\frac 12 \text{and} \ x_n\le \theta+\frac 12 \\ \Rightarrow x_n-\frac 12\le \theta\le x_1+\frac 12$$
(Parenthesis: For $\hat\theta$ to be able to obtain a value in the reals this interval must exist. So we must have $$x_n-\frac 12\le x_1+\frac 12 \Rightarrow x_n-x_1 \le 1$$
In other words, if in the sample at hand the difference between the maximum realized value and the minimum realized value is greater than unity, then one of the initial assumptions is wrong: either the n rv's do not have the postulated distribution, or they are not identically distributed). Assume that $x_n-x_1 \le 1$ holds. Then our maximization problem has been transformed (after this first step of optimization) as
$$ \max_{\theta}EI_{[x_n-\frac 12\ \, \le \, \theta \, \le \, x_1+\frac 12]}$$
The $\min$ operator has disappeared because all indicator functions have become identical.
Now by the properties of indicator functions,
$$ EI_{[x_n-\frac 12\ \, \le \, \theta \, \le \, x_1+\frac 12]} = P\left(\left[x_n-\frac 12\ \, \le \, \theta \, \le \, x_1+\frac 12\right]\right)$$
where P() is the probability of the interval. Next we can write the upper bound for $\theta$ as a constraint and solve the constrained maximization problem
$$ \max_{\theta}P\left(\left[x_n-\frac 12\ \, , \, \theta\right]\right) \qquad \text{s.t.}\ \theta\le x_1+\frac 12$$
Obviously this probability increases with the length of the interval, and so
$$ \hat \theta : [x_n-\frac 12\ , \hat \theta] = \text{maximum}$$
which obviously results in
$$ \hat \theta_{\text{ML}} = x_1+\frac 12$$
Then the maximum likelihood estimate of the common distribution of the n rv's is $$U\left(x_1\, , \, x_1+1\right).$$
It is here I am stuck. The individual probabilities of P(Xi=xi) are 0 for continuous pdf without a continuity correction. So how do we justify doing above step? What is the notion I am missing here?
You need to use probability density instead of probabilities. Definition of MLE:
The likelihood function is always defined as a function of the parameter ${\displaystyle \theta }$ equal to (or sometimes proportional to) the density of the observed data with respect to a common or reference measure, for both discrete and continuous probability distributions. So, your step is completely justified.
True that the individual probabilities are zero and to get probabilities from pdfs you need to integrate over some interval. Have a look at this answer on CV.
$L(\theta_1,\theta_2) = p(X_1=x_1;X_2=x_2;\cdots;X_m=x_m|\theta_1,\theta_2) = \prod_{i=1}^{m} \dfrac{1}{\sqrt{2\pi\theta_2}}{\text{exp}}{\Big[ -\dfrac{ (x_i-\theta_1)^2 }{2\theta_2} \Big]} \tag{1}$
where $\theta_1,\theta_2$ are mean and variance respectively.
And $p(X_1=x_1;X_2=x_2;\cdots;X_m=x_m)$ is the joint density.
And, the likelihood will be:
$L(\theta_1,\theta_2) = p(X_1=x_1;X_2=x_2;\cdots;X_m=x_m|\theta_1,\theta_2) =\dfrac{1}{\sqrt{(2\pi\theta_2)^m}}{\text{exp}}{-\frac{1}{2\theta_2}}{\sum_{i = 1}^{m}{(x_i-\theta_1)^2 }} \tag{2}$
To get the values of $\theta_1,\theta_2$, maximize the likelihood by taking derivatives w.r.t the desired variable.
Instead we are interested in a collective probability density function of all samples' individual probability densities.
I doubt this statement. You are interested in pdf of the sample from a normal population. But, you are missing that sample points are realizations of a distribution. They are not random variables but rather a realization of random variables.
From wiki:
A sample concretely represents the results of n experiments in which the same quantity is measured. For example, if we want to estimate the average height of members of a particular population, we measure the heights of n individuals. Each measurement is drawn from the probability distribution F characterizing the population, so each measured height ${\displaystyle x_{i}} $ is the realization of a random variable ${\displaystyle X_{i}}$ with distribution F. Note that a set of random variables (i.e., a set of measurable functions) must not be confused with the realizations of these variables (which are the values that these random variables take). In other words, ${\displaystyle X_{i}}$ is a function representing the measurement at the i-th experiment and ${\displaystyle x_{i}=X_{i}(\omega )}$ is the value obtained when making the measurement.
Best Answer
Some discussion here and here.
Long story short, the likelihood function should be thought of as some measure of goodness of fit of a model to some data, and this quantity should be comparable for different models for the same data. See this Wikipedia section for a description of how a likelihood function is generally defined. In the special case of discrete distributions, this general notion happens to coincide with "the probability of seeing my data under this model," but more generally it is not a probability, as you note. It is probably best not to think of likelihoods as probabilities.