Simple definition for random polynomial $p:[0,1]\to[0,1]$ This should be a relatively simple problem but I'm looking for a clean definition for a random polynomial, $p$, that passes through the unit square. From suggestions from the comments, the order of the polynomial, $d$, can be whatever is needed and $p$ can be found by any means.
Now I've tried the obvious approach of just taking random points from the square and interpolating from them but sometimes, $p$ tends to jump out of the square.
Half of the time, it's fine and $p$ is within the bounds. The other half of the time, it isn't.
So how can I either mitigate this or redefine $p$ so that it fits these criteria?
That one's good.

That one's not.

Much obliged, Jam.
 A: Depending on how "well-distributed" you want your polynomials to be, here's two simple things that seem worth trying:
1) Every Chebyshev polynomial maps $[-1,1]$ to itself, so taking an arbitrary convex combination of them will always stay within the square $[-1,-1]^2$, which you can easily rescale to fit inside the unit square.  The basic Chebyshev polynomials have the nice property that they fit exactly inside the square, hitting the top and bottom edges a maximal number of times for their degree.
2) Bernstein polynomials should behave better than Lagrange interpolants in terms of boundedness over the entire interval.  They will allow you to find a high-degree polynomial that stays within a bounded envelope of any given continuous curve that you draw in the square.  I believe that by convexity they should stay within the unit square if the original curve is within the unit square; if not you could always rescale by a small amount.
A: Though Erick Wong's provided a great solution, I'll post one I came up with, for posterity.

Let $F_n(x)=\sum_{i=1}^{n}a_i\cos(ix+\phi_i)$ with $a_i$ in $[0,1]$ and $\phi_i$, in $\mathbb{R}$. Then, for all $x$, $\frac{F_n(x)+1}{2\sum a_i}$, will be in $[0,1]$ since $\cos$ is in $[-1,1]$. Lastly, define $p(x)$ as the Taylor Series for $\frac{F_n(x)+1}{2\sum{a_i}}$, at $x=\frac{1}{2}$, with a degree sufficiently high, such that $p$ stays in $[0,1]$.

It's not quite as clean as I'd hoped but works quite well and can probably be simplified. The idea came from a polynomial fitting a random Fourier Series, without terms in $\sin$ since they complicated things. The phase shifts were necessary since each $F_n$ was coming out of the $y$ axis at predictable points.
In the diagram below, I've shown, superimposed, about $20$ to $30$ degree $9$ polynomials, using this method, with $n=4$ and $\phi_i$ in $[-3,3]$. Clearly they all stay within the confines of the unit square, fill up the space quite well and are pretty random so I'd consider this method a success.

When we see what it looks like with about $100$ or so of $p$, it's clear that they cluster around $y=\frac{1}{2}$, which is to be expected, so the method doesn't provide polynomials with an equal chance of being at every point in the unit square but they're reasonably random nonetheless.

