Is there a name for an extension of Bernstein basis polynomial based on multinomial distribution? I wonder if we can replace the binomial terms in Bernstein polynomial with multinomial terms. To be more specific, given a vector $f=(f_1, f_2,\dotsc, f_M)\in \mathbb{Z}^M$ such that $\forall m, f_m\ge0, \sum_mf_m=N\in\mathbb{Z}$, define
$F_{N,f}(\sigma)=\frac{N!}{\prod_{m} f_{m}!}\prod_{m:\sigma_m>0}(\sigma_{m})^{f_{m}}$, where $\sigma\in\Delta^{M-1}$.
I think when $M=2$, this is Bernstein polynomial. Is there a name for $F_{N, f}$ for cases $M>2$ and what are its properties?
 A: $\newcommand{\de}{\delta}\newcommand{\De}{\Delta}\newcommand\R{\mathbb R}\newcommand{\om}{\omega}\newcommand{\si}{\sigma}$Such polynomials are called multivariate Bernstein polynomials.
Many of the properties of Bernstein polynomials should be extendible to these multivariate versions.
Perhaps the most interesting of the properties of Bernstein polynomials is that they provide an explicit bound on the uniform approximation error of continuous functions by polynomials.
Let us show that such a property holds in the multivariate case as well.
Let $g\colon\De^{M-1}\to\R$ be any continuous function with the modulus of continuity $\om$ given by the formula
\begin{equation*}
    \om(\de):=\sup\{|g(x)-g(y)|\colon\|x-y\|_2\le\de\}
\end{equation*}
for real $\de>0$, where $\|\cdot\|_2$ is the Euclidean norm.
For any $\si=(\si_1,\dots,\si_M)\in\De^{M-1}$, let a random $M$-tuple $X=(X_1,\dots,X_M)$ have the multinomial distribution with parameters $N;\si_1,\dots,\si_M$. Then
\begin{equation*}
    F_{N,f}(\si)=P(X_1=f_1,\dots,X_M=f_M)
\end{equation*}
for $f=(f_1,\dots,f_M)\in\{0,1,\dots\}^M$ with $\sum_{j=1}^M f_j=N$. So,
\begin{equation*}
    p_{N,f}(\si):=Eg\Big(\frac XN\Big)=\sum_f g\Big(\frac fN\Big)F_{N,f}(\si),
\end{equation*}
and $p_{N,f}$ is a polynomial; of course, here $\sum_f$ denotes the summation over all $f=(f_1,\dots,f_M)\in\{0,1,\dots\}^M$ with $\sum_{j=1}^M f_j=N$.
Moreover, for any real $\de>0$ and all $\si\in\De^{M-1}$,
\begin{equation*}
\begin{aligned}
|p_{N,f}(\si)-g(\si)|&\le E\Big|g\Big(\frac XN\Big)-g(\si)\Big| \\ 
&=E\Big|g\Big(\frac XN\Big)-g(\si)\Big|\,1\Big(\Big\|\frac XN-\si\Big\|_2\le\de\Big) \\ 
&+E\Big|g\Big(\frac XN\Big)-g(\si)\Big|\,1\Big(\Big\|\frac XN-\si\Big\|_2>\de\Big) \\ 
&\le\om(\de)+2\|g\|_\infty\,P\Big(\Big\|\frac XN-\si\Big\|_2>\de\Big) \\ 
&\le\om(\de)+\frac{2\|g\|_\infty}{\de^2}\,E\Big\|\frac XN-\si\Big\|_2^2 \\ 
&\le\om(\de)+\frac{2\|g\|_\infty}{\de^2}\,\sum_{j=1}^M E\Big(\frac{X_j}N-\si_j\Big)^2 \\ 
&=\om(\de)+\frac{2\|g\|_\infty}{\de^2\,N}\,\sum_{j=1}^M \si_j(1-\si_j) \\ 
&\le\om(\de)+\frac{2\|g\|_\infty}{\de^2\,N}, 
\end{aligned}
\tag{1}\label{1}    
\end{equation*}
where $\|g\|_\infty:=\max_{\si\in\De^{M-1}}|g(\si)|$.
The penultimate inequality and the last equality in \eqref{1} follow because the $X_j$'s are negatively correlated random variables with respective binomial distributions with parameters $N,\si_j$, and the last equality in \eqref{1} follows because $\si\in\De^{M-1}$.
So, $\limsup_{N\to\infty}\|p_{N,f}(\si)-g(\si)\|_\infty\le\om(\de)$, for any real $\de>0$. Since $g$ is continuous, it is uniformly continuous on $\De^{M-1}$, so that $\om(\de)\to0$ as $\de\downarrow0$. We conclude that
\begin{equation}
    \lim_{N\to\infty}\|p_{N,f}(\si)-g(\si)\|_\infty=0,
\end{equation}
as claimed. $\quad\Box$
