Let $Y_{1},Y_{2},\ldots,Y_{n}$ be independent random variables with expected values $\mu_{1},\mu_{2},\ldots,\mu_{n}$, respectively. Suppose that the $\mu_{i}$'s are functions of the parameter vector that we want to estimate $\beta = [\beta_{1},\beta_{2},\ldots,\beta_{n}]^{T}$, $p < n$. Thus \begin{align*} E(Y_{i}) = \mu_{i}(\beta) \end{align*}
The simplest form of the method of least squares consists of finding the estimator $\hat{\beta}$ that minimizes the sum of squares of the differences between $Y_{i}$'s and their expected values \begin{align*} S = \sum[Y_{i}-\mu_{i}(\beta)]^{2} \end{align*}
My question is quite simple and naive, but I still would like to make it: what motivates the application of this method in such general context where there is no prior distribution related to $Y_{i}$ neither an explicit relation between the $Y_{i}$'s and the $\mu_{i}$'s?
I am doing this question because, in the normal linear model, we have the assumptions \begin{align*} \textbf{Y} = \textbf{X}\beta + \epsilon \end{align*}
where $\epsilon\sim\mathcal{N}(\textbf{0},\sigma^{2}\textbf{I})$ and $\mu = \textbf{X}\beta$. Since $\epsilon_{i} = Y_{i} - \mu_{i}$, it sounds reasonable to apply the LSE since we are interested in minimizing the corresponding error.
Hopefully, I have been clear enough at posing my doubts. Any help is appreciated. Thanks!