This is a follow-up question to 1 and 2. So we have the normal linear model \begin{align*} \textbf{Y} = \textbf{X}\beta + \epsilon \end{align*}
where $\epsilon\sim\mathcal{N}(\textbf{0},\sigma^{2}\textbf{I})$, $\mu_{i} = \beta_{0} + \beta_{1}x_{i1} + \ldots + \beta_{p}x_{ip}$ and $\mu = \textbf{X}\beta$. As far as I have understood, we take $n$ observations \begin{align*} Y_{1} & = \beta_{0} + \beta_{1}x_{11} + \ldots + \beta_{p}x_{1p} + \epsilon_{1}\\ Y_{2} & = \beta_{0} + \beta_{1}x_{21} + \ldots + \beta_{p}x_{2p} + \epsilon_{2}\\ &\vdots\\ Y_{n} & = \beta_{0} + \beta_{1}x_{n1} + \ldots + \beta_{p}x_{np} + \epsilon_{n}\\ \end{align*}
and apply the least square method, for instance, to obtain $\hat{\beta} = (\textbf{X}^{T}\textbf{X})^{-1}\textbf{X}^{T}\textbf{Y}$.
The problem which concerns me is the interpretation of such process. Let us suppose, for example, that $Y$ represents the income, $x_{1}$ indicates the gender and $x_{2}$ stands for the age. Thus we draw someone from target population and obtain the first observation
\begin{align*} Y_{1} = \beta_{0} + \beta_{1}x_{11} + \beta_{2}x_{12} + \epsilon_{1} \end{align*}
After so, we draw another person (with replacement) from the target population and obtain the second observation \begin{align*} Y_{2} = \beta_{0} + \beta_{1}x_{21} + \beta_{2}x_{22} + \epsilon_{2} \end{align*}
We repeat such process until $n$ observations are made. Once we have $\hat{\beta}$ at hand, we can estimate $\mu$ according to $\hat{\mu} = \textbf{X}\hat{\beta}$. Moreover, we can also estimate the variance $\sigma^{2}$ through the estimator \begin{align*} S^{2} = \frac{(\textbf{Y} - \textbf{X}\hat{\beta})^{T}(\textbf{Y} - \textbf{X}\hat{\beta})}{n-p-1} \end{align*}
where it is assumed that $\operatorname{Rank}(\textbf{X}) = p+1$ and $\textbf{X}$ has full rank.
My first question is: am I describing the observation process rightly?
My second question is: how should we interpret the distribution $\textbf{Y} = \mathcal{N}(\textbf{X}\beta,\sigma^{2}\textbf{I})$?
The last question may be confusing me because, in the context of inference, we normally assume the sample consists in independent identically distributed random variables and in the multiple linear regression problem we just assume independence. In other words, the distribution of $\textbf{Y}$ corresponds to the distribution of the sample $(Y_{1},Y_{2},\ldots,Y_{n})$ and the means $\mu_{i}$ do not need to be the same. Is it correct?
Any help is appreciated. Thanks in advance!