In linear regression, each predicted value is assumed to have been picked from a normal distribution of possible values. See below.
But why is each predicted value assumed to have come from a normal distribution? How does linear regression use this assumption? What if possible values are not normally distributed?
Linear regression by itself does not need the normal (gaussian) assumption, the estimators can be calculated (by linear least squares) without any need of such assumption, and makes perfect sense without it.
But then, as statisticians we want to understand some of the properties of this method, answers to questions such as: are the least squares estimators optimal in some sense? or can we do better with some alternative estimators? Then, under the normal distribution of error terms, we can show that this estimators are, indeed, optimal, for instance they are "unbiased of minimum variance", or maximum likelihood. No such thing can be proved without the normal assumption.
Also, if we want to construct (and analyze properties of) confidence intervals or hypothesis tests, then we use the normal assumption. But, we could instead construct confidence intervals by some other means, such as bootstrapping. Then, we do not use the normal assumption, but, alas, without that, it could be we should use some other estimators than the least squares ones, maybe some robust estimators?
In practice, of course, the normal distribution is at most a convenient fiction. So, the really important question is, how close to normality do we need to be to claim to use the results referred to above? That is a much trickier question! Optimality results are not robust, so even a very small deviation from normality might destroy optimality. That is an argument in favour of robust methods. For another tack at that question, see my answer to Why should we use t errors instead of normal errors?
Another relevant question is Why is the normality of residuals "barely important at all" for the purpose of estimating the regression line?
EDIT
This answer led to a large discussion-in-comments, which again led to my new question: Linear regression: any non-normal distribution giving identity of OLS and MLE? which now finally got (three) answers, giving examples where non-normal distributions lead to least squares estimators.
There is no deep reason for it, and you are free to change the distributional assumptions, moving to GLMs, or to robust regression. The LM (normal distribution) is popular because its easy to calculate, quite stable and residuals are in practice often more or less normal.
As any regression, the linear model (=regression with normal error) searches for the parameters that optimize the likelihood for the given distributional assumption. See here for an example of an explicit calculation of the likelihood for a linear model. If you take the log likelihood of a linear model, it turns out to be proportional to the sum of squares, and the optimization of that can be calculated quite conveniently.
If you want to fit a model with different distributions, the next textbook steps would be generalized linear models (GLM), which offer different distributions, or general linear models, which are still normal, but relax independence. Many other options are possible. If you just want to reduce the effect of outliers, you could for example consider robust regression.
This discussionWhat if residuals are normally distributed, but y is not? has well addressed this question.
In short, for a regression problem, we only assume that the response is normal conditioned on the value of x. It is not necessary that the independent or response variables are independent.
Let me stick to the case of a one variable regression. The details are the same, but the notation is more cumbersome in the case of a multivariate regression. Given any data set $(x_i,y_i)$ one can find the 'least squares line' $ y = \beta x +c$ , that is find $\beta$ so that $\sum_i (y_i - \sum_i \beta x_i - c)^2$ is minimized. That is pure mathematics. However under the assumption that the residuals $ \eta_i = y_i - (\beta x_i +c) $ are independent identically distributed gaussian variables with a common variance, then one can get statistical estimates of how accurate the point estimate $\beta$. In particular, one can construct the 95% confidence interval for $\beta$. After all we are assuming that we are sampling from the underlying (true) distribuion and hence if we sampled again, we should expect to get a, possibly just slightly, different answer. In particular, the p-value is the probability of observing the given $\beta$ under the hypothesis that the true value of $\beta$ is zero. So the statistics comes about as information about how accurate is the point estimate $\beta$ . What to do in the case one doesn't have statistical properties of the error term ? With apologies to "The Graduate" - one word bootstrap.
