Perhaps a partial answer will help to focus the issue of defining a probability distribution governing "random" choice of polynomials.
The region of feasible polynomials, increasing on $[0,1]$ of degree (at most) $p \gt 0$ is unbounded in the coefficient space. For one thing any constant can be added to such a polynomial and obtain another; similarly any positive constant can be multiplied to give another feasible polynomial.
If the feasible region is unbounded, we cannot glibly appeal to "uniform" distribution of probability over the region. This motivates adding a couple of constraints that make the region bounded. If strongly monotone increasing is the intended meaning, then consider adding the normalizing requirements $y(0) = 0$ and $y(1) = 1$, i.e. that $y$ maps $[0,1]$ onto $[0,1]$.
Now the case $p = 1$ becomes trivial, since $y = x$ is the only feasible polynomial.
Moreover the tests proposed by the OP, showing monotonicity by sampling the polynomial at intermediate points $\{x_1,...,x_n\}$, can be carried out in a deterministic way with $n = p-1$.
Corrected: Consider as an illustration the case $p = 2$, quadratic polynomials. Sampling these at $x_1 = \frac{1}{2}$ gives a parameterization of the solutions with $y(x_1) \in [\frac{1}{4},\frac{3}{4}]$. If we wish, the probability distribution can be taken to make $x_1$ be uniform on that interval, or to accord with some other distribution. The feasible quadratics are convex combinations of $y_0(x) = x^2$ and $y_1(x) = 2x - x^2$.
Conjecture: Looking at the lower degree examples, the general degree (at most) $p$ case appears to be this:
$ y(x) = x^p + \sum_{k=1}^{p-1} c_k \binom{p}{k} x^{p-k} (1-x)^k $
where $c_k \in [0,1]$ for $k=1,\ldots,p-1$, so that the normalized feasible polynomials can be parameterized by $[0,1]^{p-1}$, using perhaps a uniform distribution over that hypercube. Among other evidence this set of functions is conserved by the monotone-preserving symmetry $y(x) \mapsto 1 - y(1-x)$.