Let $f\colon [0,1]^n\to [0,\infty)$ be a symmetric (continuous) function, i.e., $f(x_1,...,x_n)=f\left(x_{\sigma(1)},...,x_{\sigma(n)}\right)$ for every permutation $\sigma$ of $\{1,...,n\}$ such that the following hold:
For any $i\in \{1,...,n\}$ and any $(a_1,...,a_{i-1}, a_{i+1},...,a_n)\in [0,1]^{n-1}$, the function $$[0,1]\ni t\longmapsto f(a_1,...,a_{i-1},t,a_{i+1},...,a_n)\in [0,\infty)$$ is monotone increasing.
Question: Given any $\varepsilon>0$, is it possible to find a symmetric polynomial $p$ with non-negative coefficients so that $$\big|p(\mathbf x)-f(\mathbf x)\big|<\varepsilon\text{ for all }\mathbf x\in [0,1]^n?$$
Thought: If $q\colon [0,1]^n\to \Bbb R$ is a polynomial approximation of $f$ (Weierstrass's theorem), then the following symmetric polynomial $$p(x_1,x_2,..,x_n):=\frac 1 {n!} \sum_{\sigma\in \text{sym}(n)} q(x_{\sigma(1)},x_{\sigma(2)},...,x_{\sigma(n)})$$ also approximates $f$. So all we need to know is whether $q$ has non-negative coefficients or not.