Estimates of multiple linear regression parameters I have a question about the statistical and least squares interpretation of regression.  For the simple linear regression case, suppose we have:
$$
y= \alpha + \beta x + \epsilon
$$
We are interested in an estimate of $E(y|x)=\alpha+\beta x$.   Using the law of total expectations and some algebra, we get the following as the equation for $E(y|x)$.
$$
E(y|x) = \left(E(y)-\dfrac{Cov(y,x)}{Var(x)}E(x)\right) + \dfrac{Cov(y,x)}{Var(x)}x
$$
Since we don't know the probability distributions of $(y,x)$, we replace the variances, expectations, and covariances on the right side with their sample counterparts in order to estimate the conditional expectation, which gives us the familiar estimates that we get from minimizing least squares:
$$
{\displaystyle {\begin{aligned}{\hat {\alpha }}&={\bar {y}}-{\hat {\beta }}\,{\bar {x}},\\{\hat {\beta }}&={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}\\\end{aligned}}}
$$
I guess I have two questions here.  In the multiple linear regression case, we have the following true model:
$$
y = \mathbf{x}^\intercal\boldsymbol\beta + \epsilon
$$
with $K$ predictors in $\mathbf{x}$.
The first question is, how do we get an analogous equation for $E(y|\mathbf{x})$? It seems that it should be something like:
$$
E(y|\mathbf{x})=E(y)+Cov(y,\mathbf{x})(Var(\mathbf{x}))^{-1} \left(\mathbf{x}−E(\mathbf{x})\right) .
$$
Once we have the equation, how can we get, like in the simple linear regression case where we replaced first and second moments with their sample counterparts, the familiar OLS result
$$
\hat{E(y|x)} = X\hat\beta = X(X^\intercal X)^{-1}X^\intercal y
$$
 A: please forgive the change of notation - will not have time to rewrite the eqs
One "compact" way to get the normal equations 
$$ X^TX \theta = X^Ty $$ 
and from there:
$$ \theta = (X^TX)^{-1}X^Ty $$
is to start from the cost function
$$J(\theta) = \frac{1}{2} \sum_{i=1}^{m} ( h(x_i) - y_i )^2 $$
We know that $J(\theta)$ is convex and we can minimize it by setting its derivative wrt $\theta$ to zero. You can see the full proof here.
Pls note that $x_i$ is itself a vector while $X$ is a matrix where each row corresponds to $x_i^T$
If you are looking for a probabilistic interpretation:
For the sample $i$, we have:
$$ y_i = \theta^T x_i + \epsilon_i $$
If we make the  the assumption:


*

*$\epsilon_i$ are i.i.d and Normally distributed with mean zero and variance $\sigma^2$


$$ p(y_i | x_i; \theta) = \frac{1}{\sqrt{2 \pi} \sigma} \exp \bigg (-\frac{(y_i - \theta^T x_i)^2}{2 \sigma^2} \bigg )$$
The line of reason is clear: we got the mean of the distribution from $\theta^T x_i$ and the distribution around that depends on the error term.
If we look at our entire dataset, we have to get the likelihood given a specific $\theta$.
$$L(\theta) = \prod_{i=1}^{m} p(y_i | x_i; \theta) $$
Our best fit is the maximum likelihood, i.e. the $\theta$ that maximizes the log likelihood (we take the log to make the math easier).
Check the full proof here.
