[proofplan]
Fix $n \in \mathbb{N}$; with Androma's convention $\mathbb{N}=\{1,2,\dots\}$, the Bernstein nodes $k/n$ are defined for every $k \in \{0,\dots,n\}$. We write the Bernstein polynomial in the Bernstein basis and compute its first and second derivatives in terms of first and second finite differences of the sampled values $f(k/n)$. The Bernstein basis functions are non-negative on $[0,1]$, so the sign of each derivative is controlled by the corresponding finite differences. Monotonicity of $f$ gives non-negative first differences, while convexity of $f$ gives non-negative second differences at equally spaced nodes.
[/proofplan]
[step:Introduce the Bernstein basis and its derivative identity]
For each integer $m \geq 0$ and each $j \in \{0,\dots,m\}$, define the Bernstein basis polynomial $p_{m,j}: [0,1] \to \mathbb{R}$ by
\begin{align*}
p_{m,j}(x) = \binom{m}{j}x^j(1-x)^{m-j}.
\end{align*}
For notational boundary terms, set $p_{m,j}: [0,1] \to \mathbb{R}$ to be the zero function whenever $j < 0$ or $j > m$.
For $m \geq 1$ and $j \in \{0,\dots,m\}$, differentiating the product $x^j(1-x)^{m-j}$ gives, for $x \in (0,1)$,
\begin{align*}
p_{m,j}'(x) = m\bigl(p_{m-1,j-1}(x)-p_{m-1,j}(x)\bigr).
\end{align*}
Indeed, using $j\binom{m}{j}=m\binom{m-1}{j-1}$ and $(m-j)\binom{m}{j}=m\binom{m-1}{j}$, the derivative is exactly the displayed expression.
[/step]
[step:Express the first derivative as a positive combination of first differences]
Define the sampled values $a_k \in \mathbb{R}$ by
\begin{align*}
a_k = f\left(\frac{k}{n}\right)
\end{align*}
for $k \in \{0,\dots,n\}$. Since
\begin{align*}
B_n f = \sum_{k=0}^{n} a_k p_{n,k},
\end{align*}
the derivative identity gives, for $x \in (0,1)$,
\begin{align*}
(B_n f)'(x) = n\sum_{k=0}^{n} a_k\bigl(p_{n-1,k-1}(x)-p_{n-1,k}(x)\bigr).
\end{align*}
Reindexing the first sum by replacing $k$ with $k+1$, and using the convention that out-of-range basis functions vanish, yields
\begin{align*}
(B_n f)'(x) = n\sum_{k=0}^{n-1}(a_{k+1}-a_k)p_{n-1,k}(x).
\end{align*}
[guided]
The purpose of this computation is to make monotonicity visible. A derivative of a polynomial is easy to compute, but we need it in a form whose sign can be read from the values of $f$ at the Bernstein nodes.
Define $a_k \in \mathbb{R}$ by
\begin{align*}
a_k = f\left(\frac{k}{n}\right)
\end{align*}
for $k \in \{0,\dots,n\}$. Then
\begin{align*}
B_n f = \sum_{k=0}^{n} a_k p_{n,k}.
\end{align*}
Using the derivative identity for the basis polynomials, we obtain for every $x \in (0,1)$,
\begin{align*}
(B_n f)'(x) = \sum_{k=0}^{n} a_k p_{n,k}'(x)
\end{align*}
and therefore
\begin{align*}
(B_n f)'(x) = n\sum_{k=0}^{n} a_k\bigl(p_{n-1,k-1}(x)-p_{n-1,k}(x)\bigr).
\end{align*}
Now separate the two sums. In the term containing $p_{n-1,k-1}$, replace the index $k$ by $k+1$. The convention that $p_{n-1,j}$ is the zero function outside $j \in \{0,\dots,n-1\}$ removes the boundary terms. Hence
\begin{align*}
\sum_{k=0}^{n} a_k p_{n-1,k-1}(x) = \sum_{k=0}^{n-1} a_{k+1}p_{n-1,k}(x).
\end{align*}
Similarly,
\begin{align*}
\sum_{k=0}^{n} a_k p_{n-1,k}(x) = \sum_{k=0}^{n-1} a_k p_{n-1,k}(x).
\end{align*}
Subtracting these two expressions gives
\begin{align*}
(B_n f)'(x) = n\sum_{k=0}^{n-1}(a_{k+1}-a_k)p_{n-1,k}(x).
\end{align*}
This is the desired form: the derivative is a weighted sum of first finite differences of the sampled values of $f$.
[/guided]
[/step]
[step:Use non-negative first differences to prove monotonicity]
Assume that $f$ is increasing on $[0,1]$. For every $k \in \{0,\dots,n-1\}$,
\begin{align*}
\frac{k}{n} \leq \frac{k+1}{n},
\end{align*}
so
\begin{align*}
a_{k+1}-a_k = f\left(\frac{k+1}{n}\right)-f\left(\frac{k}{n}\right) \geq 0.
\end{align*}
Also $p_{n-1,k}(x) \geq 0$ for every $x \in [0,1]$. Therefore
\begin{align*}
(B_n f)'(x) \geq 0
\end{align*}
for every $x \in (0,1)$. Since $B_n f$ is continuous on $[0,1]$ and differentiable on $(0,1)$, the non-negativity of its derivative on $(0,1)$ implies that $B_n f$ is increasing on $[0,1]$.
[/step]
[step:Express the second derivative as a positive combination of second differences]
If $n=1$, then $B_1 f$ is affine and hence convex. Assume now that $n \geq 2$. Differentiating the first-derivative formula gives, for $x \in (0,1)$,
\begin{align*}
(B_n f)''(x) = n(n-1)\sum_{k=0}^{n-1}(a_{k+1}-a_k)\bigl(p_{n-2,k-1}(x)-p_{n-2,k}(x)\bigr).
\end{align*}
Reindexing as before yields
\begin{align*}
(B_n f)''(x) = n(n-1)\sum_{k=0}^{n-2}(a_{k+2}-2a_{k+1}+a_k)p_{n-2,k}(x).
\end{align*}
[/step]
[step:Use non-negative second differences to prove convexity]
Assume that $f$ is convex on $[0,1]$. For each $k \in \{0,\dots,n-2\}$, the point $(k+1)/n$ is the midpoint of $k/n$ and $(k+2)/n$. By convexity,
\begin{align*}
f\left(\frac{k+1}{n}\right) \leq \frac{1}{2}f\left(\frac{k}{n}\right)+\frac{1}{2}f\left(\frac{k+2}{n}\right).
\end{align*}
Equivalently,
\begin{align*}
a_{k+2}-2a_{k+1}+a_k \geq 0.
\end{align*}
Since $p_{n-2,k}(x) \geq 0$ for every $x \in [0,1]$, the second-derivative formula gives
\begin{align*}
(B_n f)''(x) \geq 0
\end{align*}
for every $x \in (0,1)$. Thus $B_n f$ is convex on $(0,1)$, and continuity on $[0,1]$ extends the convexity inequality to all pairs of points in $[0,1]$. Hence $B_n f$ is convex on $[0,1]$.
[/step]
