Two random variables have the same $n^{th}$ moment for all $n\geq 0$. Then they have the same distribution

Question

I am working this question:

Let $X$ and $Y$ be $[0,1]-$valued random variables such that $E[X^{n}]=E[Y^{n}]$ for every integer $n\geq 0$. Show that $E[f(X)]=E[f(Y)]$ for every continuous function $f:[0,1]\longrightarrow\mathbb{R}$ and conclude that $X=_{d}Y$. (Hint: use the Weierstrass approximation theorem)

Similar questions have been posted here:Show two random variables have same distribution

and here: Proof of if two random variables have the same distribution then they have the same moment generating function.

However, it seems that this question is asking me firstly to prove $E[f(X)]=E[f(Y)]$ for every continuous function, and then use this fact to deduce they have the same distribution.

The Weierstrass approximation theorem is here:

If $f$ is a real-valued function on $[a,b]$ and if any $\epsilon>0$ is given, then there exists a polynomial $P$ on $[a,b]$ such that $|f(x)-P(x)|<\epsilon$ for all $x\in [a,b]$.

So firstly I show that $n^{th}$ moment of random variable is always a polynomial? Even if I showed this, it seems that the argument is backward.

Answer with some more details will be really appreciated since I got really lost...

Thank you!

Edit 1:

As what the comment suggested, I should show $E[p(X)]=E[p(Y)]$ for all polynomial. However, I am reading Durrett, and there is nothing related to the relation between $n^{th}$ moment of a random variable and the expectation of a polynomial. I don't really know what to do here.

Also, as I pointed out in the comment of the first answer, if I know they have the same $n^{th}$ moment for each $n$, then using series expansion to the moment generating function immediately yields me that they have the same moment generating function and thus they have the same distribution. So why would I need to deduce this fact from $E[f(X)]=E[f(Y)]$?

Answer 1

Take an arbitrary polynomial $p(x) = c_n x^n + c_{n-1}x^{n-1}+ \dots + c$. Then

$$ \mathbb{E}[p(X)] = c_n\mathbb{E}[X^n] + c_{n-1} \mathbb{E}[X^{n-1}] + \dots + c = \\ c_n\mathbb{E}[Y^n] + c_{n-1} \mathbb{E}[Y^{n-1}] + \dots + c = \mathbb{E}[p(Y)] $$

Now take any function $f : [0,1] \to \mathbb{R}$ and corresponding polynomial $p$ from Weierstrass then

$$ \left|\mathbb{E}[f(X)] - \mathbb{E}[f(Y)]\right| = \\ \left|\mathbb{E}[f(X) - p(X) + p(X)] - \mathbb{E}[f(Y) - p(Y) + p(Y)]\right| \leq \\ 2 \epsilon + \underbrace{\left|\mathbb{E}[p(X) - p(Y)] \right|}_{equal} \leq 2 \epsilon $$

Now since $X,Y$ are bounded the moment generating functions exists and since $\mathbb{E}[f(X)] = \mathbb{E}[f(Y)]$. The moment generating functions are equal and hence $X =^{d} Y$