Inductive definition of Bernstein polynomials For $n\in \mathbb{N}$ let $B_n$ be the linear operator taking a function $f$ on the unit interval $I=[0,1]$ to its $n$-th Bernstein polynomial $B_nf$, 
$$   B_nf(x):=\sum_{k=0}^n\binom{n}{k} f\Big(\frac{k}{n}\Big)x^k(1-x)^{n-k}\label{1}\tag{1}$$ 
The polynomial $B_nf(x)$ has a natural probabilistic interpretation, namely, it is the expected value of $f(\xi)$, where   $\xi=\frac{1}{n}\sum_{j=1}^n \omega_j$ is the average  of $n$ independent random variables with identical Bernoulli distribution of parameter $x$, that is, $\mathbb{P}(\omega_j=1)=x$.  In fact, this is the starting point in the beautiful Bernstein's proof of the Weierstrass' density theorem via the WLLN. However, this question is about an alternative definition of the sequence $(B_n)_{n\ge0}$. 
Let $D:C^1(I)\to C^0(I)$ be the derivative operator, and for all $n\ge1$, let  $D_n:C^0(I)\to C^0(I)$ be the approximate discrete derivative given by the incremental ratio
$$D_nf(x):=\frac{f\big( \frac{n-1}{n} x+\frac{1}{n}\big)-f\big( \frac{n-1}{n} x\big)}{\frac{1}{n}}, $$ 
(which is well-defined for  $f\in C^0(I)$ and $x\in I$). 
It is easy to check that definition \eqref{1} implies 
$$DB_n=B_{n-1}D_n\label{2}\tag{2}$$
together with:
$$B_0f(x)=B_nf(0)=f(0)\label{3}\tag{3}$$
Conversely these two imply formula \eqref{1}, as it follows immediately by induction, at least, if we already have it (quite a common situation of formulas proven by induction). Thus, since \eqref{2} and \eqref{3} characterize $(B_n)_n$, we may even take them as an inductive definition of $(B_n)_n$. Note that replacing $D_n$ with $D$ in \eqref{2} gives the analogous inductive definition for the Taylor polynomials in $0$. (Incidentally, formula \eqref{2} is relevant in the approximation theory, in that it implies that for $f\in C^k(I)$ one has $B_nf\to f$ in $C^k$: this by induction from the case $k=0$, since $D_n$ converges strongly to $D$. Also, it says that if some derivative $f^{(k)}$ is non-negative on $I$, so is $(B_nf)^{(k)}$.)
Question: How can we deduce  naturally formula \eqref{1} (i.e., assuming we don't know it, and we do not have a crystal ball to guess it) from \eqref{2} and \eqref{3}?
 A: $\newcommand{\De}{\Delta}$
Iterating your condition \eqref{2}, for $k=0,\dots,n$ we have 
\begin{equation*}
 D^kB_n=\frac{n!}{(n-k)!}\,B_{n-k}P_{n,k},\label{a}\tag{a}
\end{equation*}
where 
\begin{equation*}
 P_{n,k}:=\De_{n-k+1}\cdots\De_n,\quad \De_j:=\tfrac1j\,D_j. 
\end{equation*}
By induction on $k=0,\dots,n$, 
\begin{equation*}
 (P_{n,k}f)(x)=\sum_{i=0}^k(-1)^{k-i}\binom ki f\Big(\frac{n-k}n\,x+\frac in\Big),\label{b}\tag{b}
\end{equation*}
whence, using \eqref{a} and taking your condition \eqref{3} into account, we have 
\begin{equation*}
 \frac{(n-k)!}{n!}\,(D^kB_n f)(0)=(B_{n-k}P_{n,k}f)(0)=(P_{n,k}f)(0)
 =\sum_{i=0}^k(-1)^{k-i}\binom ki f\Big(\frac in\Big). 
\end{equation*}
Also, using again \eqref{a} and \eqref{b}, and again taking your condition \eqref{3} into account, we have 
\begin{equation*}
 \frac1{n!}\,(D^nB_n f)(x)=(B_0P_{n,n} f)(x)=(P_{n,n} f)(0) 
 =\sum_{i=0}^n(-1)^{n-i}\binom ni f\Big(\frac in\Big), 
\end{equation*}
a constant. So, $B_n f$ is a polynomial of degree $\le n$, and hence 
\begin{align*}
 (B_n f)(x)&=\sum_{k=0}^n \frac{(D^kB_n f)(0)}{k!}\,x^k \\ 
 &=\sum_{k=0}^n\binom nk x^k 
 \sum_{i=0}^k(-1)^{k-i}\binom ki f\Big(\frac in\Big) \\ 
 &=\sum_{i=0}^n f\Big(\frac in\Big)\sum_{k=i}^n (-1)^{k-i}\binom nk \binom ki x^k  \\  
 &=\sum_{i=0}^n f\Big(\frac in\Big)\binom ni x^i(1-x)^{n-i}, 
 \end{align*}
as desired. 
A: A comment  on Josif Pinelis' formula $(b)$ for $\Delta_{n-k+1} \dots\Delta_{n-1}\Delta_{n}$, which is a main point of the computation. Let $\{\tau_{a}\}_{a\in\mathbb{R}}$ and $\{\delta_{b}\}_{a\in\mathbb{R}_+}$ denote respectively the linear group of translations on functions (that we may think defined on the whole real line w.l.o.g.), $f(\cdot)\mapsto f(\cdot+a)$, and the linear group of dilations,  $f(\cdot)\mapsto f(\cdot b)$. So $$\tau_{a+b}=\tau_a\tau_b,$$ $$\delta_{ab}=\delta_a\delta_b,$$ $$\tau_{ab}=\delta_b^{-1}\tau_a\delta_b$$
Since $\Delta_n:=\delta_{\frac{n-1}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)$, moving all dilations on the left by the above relations imply nicely 
$$\Delta_{n-k+1} \dots\Delta_{n-1}\Delta_{n}=\delta_{\frac{n-k}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k},$$
whence
$$\frac{1}{k!} D^kB_n=\frac{1}{k!}B_{n-k} D _{n-k+1} \dots D _{n-1} D _{n}=\Big({n\atop k}\Big)B_{n-k}\delta_{\frac{n-k}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k},$$
which we can expand to  formula $(b)$.
edit. In fact we may skip the last expansion too, keeping all Josif's formulas on the level of operators.  Since the $D_k$'s lower the degree of polynomials, $(2)$ and $(3)$ imply that $B_n$  takes values on polynomials of degree less than or equal to $n$, as said. So, for any $x$, denoting $e_x$ the evaluation form,
$$ e_xB_n=e_0\bigg[\sum_{k=0}^n \frac{x^k}{k!}D^kB_n\bigg]=e_0\bigg[\sum_{k=0}^n  x^k \Big({n\atop k}\Big)B_{n-k}\delta_{\frac{n-k}{n}}\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k}\bigg]=$$
$$=e_0\bigg[\sum_{k=0}^n  \Big({n\atop k}\Big)x^k\big(\tau_{\frac{1}{n}}-\mathbb{1}\big)^{k}\bigg]=e_0\bigg( \mathbb{1} + x \big(\tau_{\frac{1}{n}}-\mathbb{1}\big) \bigg)^n =e_0\bigg( x  \tau_{\frac{1}{n}} + (1-x)\mathbb{1}  \bigg)^n$$
$$=e_0\bigg(\sum_{k=0}^n  \Big({n\atop k}\Big)x^k(1-x)^{n-k}\tau_{\frac{k}{n}} \bigg) $$
which indeed takes $f$ to the original $(B_nf)(x)$ given by $(1)$.
