Theorem: Weierstrass Approximation Theorem

Let $f$ be a real function which is continuous on the closed interval $\Bbb I = \closedint a b$.
Then $f$ can be uniformly approximated on $\Bbb I$ by a polynomial function to any given degree of accuracy.

Proof:
{{WLOG}}, assume $\Bbb I = \closedint 0 1$.
For each $n \in \N$, let:
:$\ds \map {P_n} x := \sum_{k \mathop = 0}^n \map f {\dfrac k n } \dbinom n k x^k \paren {1 - x}^{n - k}$
We shall show that $\ds \lim_{n \mathop \to \infty} \norm { P_n - f}_\infty = 0$.
Let $\epsilon \in \R_{>0}$.
By Heine-Cantor Theorem, there is a $\delta \in \R_{>0}$ such that:
:$\forall x,y \in \Bbb I :  \size {x - y} \le \delta \implies \size {\map f x - \map f y} \le \epsilon $
Let $p \in \Bbb I$.
Let $Z_n$ be a random variable such that:
:$n Z_n \sim \Binomial n p$
where $\Binomial n p$ denotes the binomial distribution with parameters $n$ and $p$.
Observe that:
{{begin-eqn}}
{{eqn | l = \expect {\map f {Z_n} }
      | r = \sum_{k \mathop = 0}^n \map f {\dfrac k n} \map \Pr {Z_n = \dfrac k n}
      | c = {{Defof|Expectation of Discrete Random Variable}}
}}
{{eqn | r = \sum_{k \mathop = 0}^n \map f {\dfrac k n} \map \Pr {n Z_n = k}
}}
{{eqn | r = \sum_{k \mathop = 0}^n \map f {\dfrac k n} \dbinom n k p^k \paren {1 - p}^{n - k}
      | c = {{Defof|Binomial Distribution}}
}}
{{eqn | n = 1
      | r = \map {P_n} p
}}
{{end-eqn}}
Furthermore:
{{begin-eqn}}
{{eqn | l = \map \Pr {\size { Z_n - p} > \delta }
      | o = \le
      | r = \dfrac {\expect {\size {Z_n - p } ^2} } {\delta^2}
      | c = Bienaymé-Chebyshev Inequality
}}
{{eqn | r = \dfrac {\expect {\size {n Z_n - n p} ^2} } {\delta^2 n^2}
      | c = Expectation is Linear
}}
{{eqn | r = \dfrac {\expect {\size {n Z_n - \expect {n Z_n} } ^2} } {\delta^2 n^2}
      | c = Expectation of Binomial Distribution
}}
{{eqn | r = \dfrac {\var {n Z_n} } {\delta^2 n^2}
      | c = {{Defof|Variance}}
}}
{{eqn | r = \dfrac {p \paren {1 - p} } {\delta^2 n}
      | c = Variance of Binomial Distribution
}}
{{eqn | o = \le
      | r = \dfrac 1 {4 \delta^2 n}
      | c = Cauchy's Mean Theorem
}}
{{end-eqn}}
On the other hand:
{{begin-eqn}}
{{eqn | l = \size {\map f {Z_n} - \map f p}
      | o = \le
      | r = \epsilon + \size {\map f {Z_n} - \map f p } \chi_{\set {\size {Z_n - p} > \delta } }
}}
{{eqn | o = \le
      | r = \epsilon + 2 \norm f_\infty \chi_{\set {\size { Z_n - p } > \delta} }
}}
{{end-eqn}}
Therefore:
{{begin-eqn}}
{{eqn | l = \size {\map {P_n} p - \map f p}
      | r = \size {\expect {\map f { Z_n } } - \map f p}
      | c = by $(1)$
}}
{{eqn | r = \size {\expect {\map f {Z_n} } - \expect {\map f p} }
      | c = Expectation of Almost Surely Constant Random Variable 
}}
{{eqn | r = \size {\expect {\map f {Z_n} - \map f p} }
      | c = Expectation is Linear
}}
{{eqn | o = \le
      | r = \expect {\size {\map f {Z_n} - \map f p} }
      | c = Triangle Inequality for Integrals
}}
{{eqn | o = \le
      | r = \epsilon + 2 \norm f_\infty \map \Pr {\size {Z_n - p} > \delta}
}}
{{eqn | o = \le
      | r = \epsilon + \dfrac {\norm f_\infty} {2 \delta^2 n}
}}
{{end-eqn}}
Thus for all $n \in \N_{> 2 \delta^2 / \norm f_\infty}$
we have:
:$\size {\map {P_n} p - \map f p} \le 2 \epsilon$
As the above is true for all $p \in \Bbb I$, we have:
:$\forall n \in \N_{> 2 \delta^2 / \norm f_\infty} : \norm { P_n - f}_\infty \le 2 \epsilon$
{{qed}}