Given a certain $n \in \mathbb{N}$, can i construct a discrete positive random variable $X$ that fullfill the following conditions :
$$\forall k \in \{1,...,n\}, \mathbb{E}(X^k) = \mu_k$$
for some $\mu_k$ given (parameters). We suppose that there exists such a random variable (in facts, the $\mu$'s come from a random variable that i want to estimate).
For exemple, for $n=1$, it is enough to set a dirac in $\mu_1$. for $n=2$, take a radom variable with 2 atoms $\mu_1 +x$ and $\mu_1 -x$ with equal probability and choose x such that :
$$\mu_1^2 + x^2 = \mu_2$$
giving $x = \sqrt{\mu_2 - \mu_1^2}$.
What about $n \ge 3$ ?
Edits :
1° Yes, it is possible (prooved to be possible by many diffrent papers, notably Tchakaloff's theorems). 2° Yes, there exists some algorithm to do it, but i did not found a propper one yet.
The Full answer about the existence question is given in the paper of Curto-Fialkow, which we can resume (For the case of Hamburger moment problem) in two different cases, but first we consider the notation : We put $A_k$ is the matrix $(\mu_{i+j})_{i,j\leq k}$ and $v_{k}^l$ is the vector $(\mu_{k+j})_{j\leq l}$.
So the first case is when $n=2k+1$ is odd.
Then we get a positive measure $\iff$ we get a discrete measure $\iff$ $A_k$ is positive semi-definite and $v_{k+1}^k$ is in the range of $A_k$.
The second case is when $n=2k$ is even.
Then we get a positive measure $\iff$ we get a discrete measure $\iff$ $A_k$ is positive semi-definite and $rank(A_k)=\max\{p\leq k\;; \; \det(A_p)>0\}$.
For the fact that every positive measure can be represented as an atomic measures on the space of polynomials of degree $p$, is wildly known in dimension 1 as Gaussian quadrature (Or As Richter-Tchakaloff's theorems if the dimension is bigger then one, for more detail, please check the resent book of K. Schmudgen The moment problem) which state that :
Let $\mu$ be a positive measure supported on $I\subset R$, we assume that $\int_I x^nd\mu(x)<\infty$ for all $n$, then for every fixed $p\in\mathbb{N}$, we can find a set of points of $I$, $(x_i)_{1\leq i\leq p}$ called nodes, and a sequence of positive numbers $(w_i)_{1\leq i\leq p}$ called weights, such that : $$\int_I p(x)d\mu(x)=\sum_{i=1}^p w_i p(x_i)\qquad \qquad \forall p \in \mathbb{R}_{2p-1}[X]$$ Furthemore $(x_i)_i$ are exactely the zeros of the $p^{th}$ orthogonal polynomial $P_p(x)= \det \begin{bmatrix} \mu_0 & \mu_1 & \mu_2 &\cdots & \mu_p \\ \mu_1 & \mu_2 & \mu_3 &\cdots & \mu_{p+1} \\ &&\vdots&& \vdots \\ \mu_{p-1} &\mu_p& \mu_{p+1} &\cdots &\mu_{2p-1}\\ 1 & x & x^2 & \cdots & x^p \end{bmatrix}$
So to give an answer to your second question : If $n=2k+1$, is odd, then the atomic measure will be concentrated on the zero set of $P_{k+1}(x)=\det \begin{bmatrix} \mu_0 & \mu_1 & \mu_2 &\cdots & \mu_{k+1} \\ \mu_1 & \mu_2 & \mu_3 &\cdots & \mu_{k+2} \\ &&\vdots&& \vdots \\ \mu_{k} &\mu_{k+1}& \mu_{k+2} &\cdots &\mu_{2k+1}\\ 1 & x & x^2 & \cdots & x^{k+1} \end{bmatrix}$
and after you can get your weights (for example by resolving a Vandermonde system ...) For the even case $n=2k$ you should extend your sequence to an odd one, and then use the same result.