In linear regression, are the noise terms independent of the coefficient estimators? In the Wikipedia article on the bias-variance tradeoff, the independence of the estimator $\hat f(x)$ and the noise term $\epsilon$ is used in a crucial way in the proof of the decomposition of the mean square error. No justification for this independence is given, and I can't seem to figure it out. For example, if $f(t)=\beta_0t + \beta_1$, $Y_i=f(x_i) + \epsilon_i$ ($i=1,\ldots,n$), and $\hat f(x)=\hat\beta_0 + \hat\beta_1 x$ as in simple linear regression, are the $\epsilon_i$ independent of $\hat\beta_0$ and $\hat\beta_1$? 
 A: $\hat f(x)$ and $\epsilon$ should not be independent --- only asymptotically so as the cardinality of the data set approaches infinity. Reference [7] in that Wikipedia article is a correct derivation. It, specifically Equation (8), uses the premise that $\mathbf E[\epsilon|x]=0$. The derivation in the article cited from reference [6] is wrong.
A: I've poured over the argument in the link Hans kindly provided and I think I (finally) understand what's going on.  I'd like to give a version of the argument in the simpler case where the $x$-values fixed quantities so that I can avoid conditional expectations. [Corrections/confirmation/feedback greatly appreciated!]
tl;dr The claim in the Wikipedia article that $\epsilon$ and $\hat f$ are independent is true when $\hat f$ is computed using training data, $(x_1,y_1),\ldots,(x_n,y_n)$, and $\epsilon=\epsilon_0$ is the noise term associated to an independent test measurement, $(x_0,y_0)$.
More formally, let $f(\theta, x)$ be a real-valued function, where $\theta\in\mathbb R^n$ and $x\in\mathbb R$.
Fix quantities $x_1,\ldots,x_n\in\mathbb R$. Let $\epsilon_1,\ldots,\epsilon_n$ be pairwise independent random variables with $E[\epsilon_i]=0$ and consider random variables
$$Y_i = f(\theta, x_i) + \epsilon_i.$$
Let $\hat\theta(Y_1,\ldots,Y_n)$ be an estimator of $\theta$.
For $x\in\mathbb R$, view $f(\hat\theta(Y_1,\ldots,Y_n), x)$ as an estimator of $f(\theta, x)$.
Let $x_0\in\mathbb R$, let $\epsilon_0$ be a random variable, independent of $\epsilon_1,\ldots,\epsilon_n$, with $E[\epsilon_0]=0$, and suppose
$$
Y_0 = f(\theta,x_0) + \epsilon_0.
$$
We want to determine the expected error in $f(\hat\theta(Y_1,\ldots,Y_n), x_0)$ as an estimate of the "new measurement" $Y_0$:

Theorem
$$E\left[\big(Y_0 - f(\hat\theta(Y_1,\ldots,Y_n), x_0)\big)^2\right] = E[\epsilon_0^2] + E\left[\big(f(\theta,x_0)-f(\hat\theta, x_0)\big)^2\right]$$
Proof
Write $\hat\theta$ for $\hat\theta(Y_1,\ldots,Y_n)$.
Doing some algebra and using the identity $Y_0=f(\theta,x_0)+\epsilon_0$, we get
$$
\begin{aligned}
E\left[\big(Y_0 - f(\hat\theta, x_0)\big)^2\right]
&= E\left[\big(Y_0 - f(\theta, x_0) + f(\theta, x_0) - f(\hat\theta, x_0)\big)^2\right]\\
&= E\left[\big(Y_0 - f(\theta, x_0)\big)^2\right] + E\left[\big(f(\theta, x_0) - f(\hat\theta, x_0)\big)^2\right]\\
&\qquad+ 2E\left[\big(Y_0 - f(\theta, x_0)\big)\big(f(\theta, x_0) - f(\hat\theta, x_0)\big)\right]\\
&= E\left[\epsilon_i^2\right] + E\left[\big(f(\theta, x_0) - f(\hat\theta, x_0)\big)^2\right]\\
&\qquad+ 2E\left[\epsilon_0\big(f(\theta, x_0) - f(\hat\theta, x_0)\big)\right].
\end{aligned}
$$
Since the random variable
$f(\theta, x_0) - f(\hat\theta, x_0)$
depends only on $\epsilon_1,\ldots,\epsilon_n$ (through $\hat\theta(Y_1,\ldots,Y_n)$) and $\epsilon_0$ is independent of $\epsilon_1,\ldots,\epsilon_n$,
it follows that $\epsilon_0$ and
$f(\theta, x_0) - f(\hat\theta, x_0)$ are independent.
Therefore,
$$
E\left[\epsilon_0\big(f(\theta, x_0) - f(\hat\theta, x_0)\big)\right]=
E[\epsilon_0]E\left[f(\theta, x_0) - f(\hat\theta, x_0)\right] = 0,
$$
since $E[\epsilon_0]=0$, by hypothesis.
The result follows.

The term
$$
\mathop{\text{MSE}}\big(f(\hat\theta, x_0), f(\theta,x_0)\big)
=E\left[\big(f(\theta,x_0)-f(\hat\theta, x_0)\big)^2\right]
$$
decomposes, further, as
$$
\left(E\big[f(\hat\theta, x_0)\big] - f(\theta, x_0)\right)^2 + \mathop{\text{Var}}\, f(\hat\theta, x_0),
$$
giving the bias-variance decomposition. (See below.)
Note, also, that we haven't assumed that the $\epsilon_i$ have equal variances.

Added in response to Hans's comment:
Theorem
$$
E\left[\big(f(\theta,x_0)-f(\hat\theta, x_0)\big)^2\right]
=
\left(E\big[f(\hat\theta, x_0)\big] - f(\theta, x_0)\right)^2 + \mathop{\text{Var}}\, f(\hat\theta, x_0),
$$
Proof
The idea is to let $E[f(\hat{\theta},x_0)]$ mediate between $f(\hat{\theta},x_0)$ itself and $f(\theta,x_0)$
$$
\begin{aligned}
E\left[\big(f(\theta,x_0)-f(\hat\theta, x_0)\big)^2\right]
&= E\left[\big(\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}-\left\{f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]\right\}\big)^2\right]\\
&= E\left[\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}^2\right]\\
&\qquad + E\left[\left\{f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]\right\}^2\right]\\
&\qquad - 2E\left[\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}\left\{f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]\right\}\right]\\
&=: A + B - 2C
\end{aligned}
$$
Since $\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}^2$ is just a number, it's its own expected value:
$$
A = E\left[\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}^2\right]
=\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}^2 =
\mathop{\text{Bias}}(f(\hat{\theta},x_0), f(\theta,x_0))^2.
$$
By definition,
$$
B = E\left[\left\{f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]\right\}^2\right]
= \mathop{\text{Var}}f(\hat{\theta},x_0).
$$
Since $f(\theta,x_0) - E[f(\hat{\theta},x_0)]$ is just a number and $f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]$ has expected value $0$,
$$
\begin{aligned}
C &= E\left[\left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}\left\{f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]\right\}\right]\\
&= \left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}E\left[\left\{f(\hat\theta, x_0) - E[f(\hat{\theta},x_0)]\right\}\right]\\
&= \left\{f(\theta,x_0) - E[f(\hat{\theta},x_0)]\right\}\cdot 0\\
&= 0.
\end{aligned}
$$
A: When doing ordinary linear regression then $\hat f(x)$ is the vector projection of $f(x)$ onto the column space of the model matrix. 
If $f(x)$ is distributed as a multivariate normal distribution with equal variance and no correlation then the components $\hat f(x)$ and the residual $\epsilon$ in this projection will be independent. 
But this is certainly not always the case when the variance is not equal and/or when there is correlation and/or when the error distribution is not Gaussian.
Example
To make the plotting easier we use the simplest of all models (estimating the mean) and only two data points, but you can extend this to more dimensions (in which case it becomes a projection onto a hyperplane instead of a projection onto a line).
We test use the model below for fitting.
$$f(x_i) = \beta_0+\epsilon_i$$
We generate a hundred data sets with the 'true' model: $\beta_0 = 5$, $\epsilon_1 \sim N(0,0.25)$, and $\epsilon_2 \sim N(0,4)$. Below you see this plotted. 
In the first graph you see the pairs of points $\lbrace x_1,x_2 \rbrace$ projected onto the line $x_1=x_2$ (depicted by the gray dotted lines). The projected point correponds to the sample mean $\bar{x}$.
In the second graph you see the plot of the estimator model $\hat{f}(x) = \lbrace \bar{x}, \bar{x} \rbrace$ versus the residual of the first coordinate $\epsilon_1 = x_1-\bar{x}$. 
It is obvious that the two are not independent.

The image below is for independent normal distributed variables with equal variance. $\epsilon_i \sim N(0,2)$

Now the residual term and fitted term are independent. 
A: No, they're not independent: In multiple linear regression the OLS coefficient estimator can be written as:
$$\begin{equation} \begin{aligned}
\hat{\boldsymbol{\beta}} 
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} (\mathbf{x}^\text{T} \mathbf{y}) \\[6pt]
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} (\mathbf{x} \boldsymbol{\beta} + \boldsymbol{\varepsilon}) \\[6pt]
&= \boldsymbol{\beta} + (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{\varepsilon}. \\[6pt]
\end{aligned} \end{equation}$$
In regression problems we analyse the behaviour of the quantities conditional on the explanatory variables (i.e., conditional on the design matrix $\mathbf{x}$).  The covariance between the coefficient estimators and errors is:
$$\begin{equation} \begin{aligned}
\mathbb{Cov} ( \hat{\boldsymbol{\beta}}, \boldsymbol{\varepsilon} |\mathbf{x}) 
&= \mathbb{Cov} \Big( (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{\varepsilon}, \boldsymbol{\varepsilon} \Big| \mathbf{x} \Big) \\[6pt]
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbb{Cov} ( \boldsymbol{\varepsilon}, \boldsymbol{\varepsilon} | \mathbf{x} ) \\[6pt]
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbb{V} ( \boldsymbol{\varepsilon} | \mathbf{x} ) \\[6pt]
&= \sigma^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \boldsymbol{I} \\[6pt]
&= \sigma^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T}. \\[6pt]
\end{aligned} \end{equation}$$
In general, this covariance matrix is a non-zero matrix, and so the coefficient estimators are correlated with the error terms (conditional on the design matrix).

Special case (simple linear regression): In the special case where we have a simple linear regression with an intercept term and a single explanatory variable we have design matrix:
$$\mathbf{x} 
= \begin{bmatrix}
1 & x_1 \\[6pt]
1 & x_2 \\[6pt]
\vdots & \vdots \\[6pt]
1 & x_n \\[6pt]
\end{bmatrix},$$
which gives:
$$\begin{equation} \begin{aligned}
(\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} 
&= \begin{bmatrix}
n & & \sum x_i \\[6pt]
\sum x_i & & \sum x_i^2 \\[6pt]
\end{bmatrix}^{-1} 
\begin{bmatrix}
1 & 1 & \cdots & 1 \\[6pt]
x_1 & x_2 & \cdots & x_n \\[6pt]
\end{bmatrix} \\[6pt]
&= \frac{1}{n \sum x_i^2 - (\sum x_i)^2} \begin{bmatrix}
\sum x_i^2 & & -\sum x_i \\[6pt]
-\sum x_i & & n \\[6pt]
\end{bmatrix}
\begin{bmatrix}
1 & 1 & \cdots & 1 \\[6pt]
x_1 & x_2 & \cdots & x_n \\[6pt]
\end{bmatrix} \\[6pt]
&= \frac{1}{n \sum x_i^2 - (\sum x_i)^2} \begin{bmatrix}
\sum x_i(x_i-x_1) & \cdots & \sum x_i(x_i-x_n) \\[6pt]
-\sum (x_i-x_1) & \cdots & -\sum (x_i-x_n) \\[6pt]
\end{bmatrix}. \\[6pt]
\end{aligned} \end{equation}$$
Hence, we have:
$$\begin{equation} \begin{aligned}
\mathbb{Cov}(\hat{\beta}_0, \varepsilon_k) 
&= \sigma^2 \cdot \frac{\sum x_i(x_i-x_k)}{n \sum x_i^2 - (\sum x_i)^2}, \\[10pt]
\mathbb{Cov}(\hat{\beta}_1, \varepsilon_k) 
&= - \sigma^2 \cdot \frac{\sum (x_i-x_k)}{n \sum x_i^2 - (\sum x_i)^2}. \\[10pt]
\end{aligned} \end{equation}$$
We can also obtain the correlation, which is perhaps a bit more useful.  To do this we note that:
$$\mathbb{V}(\varepsilon_k) = \sigma^2 \quad \quad \quad
\mathbb{V}(\hat{\beta}_0) = \frac{\sigma^2 \sum x_i^2}{n \sum x_i^2 - (\sum x_i)^2} \quad \quad \quad
\mathbb{V}(\hat{\beta}_1) = \frac{\sigma^2 n}{n \sum x_i^2 - (\sum x_i)^2}.$$
Hence, we have correlation:
$$\begin{equation} \begin{aligned}
\mathbb{Corr}(\hat{\beta}_0, \varepsilon_k) 
&= \frac{\sum x_i(x_i-x_k)}{\sqrt{(\sum x_i^2)(n \sum x_i^2 - (\sum x_i)^2)}}, \\[10pt]
\mathbb{Corr}(\hat{\beta}_1, \varepsilon_k) 
&= - \frac{\sum (x_i-x_k)}{\sqrt{n(n \sum x_i^2 - (\sum x_i)^2)}}. \\[10pt]
\end{aligned} \end{equation}$$
