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.
What now about the construction of such discrete r.v ? Is there some algorithm to do it (i guess not ready-made !) ? Maybe we could find something in there : https://www.ams.org/journals/mcom/1967-21-100/S0025-5718-1967-0222534-4/S0025-5718-1967-0222534-4.pdf , but i'm not shure the adaptation to this case is straightforward.
– lrnv Jul 09 '19 at 07:54