Finding probability distribution from given moments Some quantum mechanics derivation gave me the expectation values of some operator. I got some clean analytical result which makes me think that there is an analytical answer for the probability distribution.The odd moments are all $0$ and the even moments are given by:
$$\langle p^{2n}\rangle=2^{-2n}\binom{2n}{n}(k)^{2n}$$
$ghk$ is just some parameter. The most significant result for me will happen if this will result in some probability distribution based on Bessel $J$ functions but anything is good.
 A: Let us denote $q := ghk$ and your variable by $X.$ You can consider the moment-generating function of the variable $X,$ $M_X(t) = \mathbb{E}[e^{tX}] = \sum_{0 \leq k} \mathbb{E}[X^k] \frac{t^k}{k!},$ which in your case would be $$M_X(t) = \sum_{0 \leq k} 2^{-2k} \binom{2k}{k} \frac{(tq)^{2k}}{(2 k)!} = \sum_{k \geq 0} 2^{-2 k} \frac{(2 k)!}{(k!)^2} \frac{(t q)^{2k}}{(2k)!} = \sum_{k \geq 0} (\frac{t^k q^k}{2^k k!})^2.$$ Note that we already get that $M_X(t) \leq e^{tq},$ so this series is actually well defined in this case (it is not always the case with moment generating functions),. In fact, this is just a modified Bessel function of the first kind. More precisely, $M_X(t) = I_0(t q),$ so your intuition was right that we would get some kind of Bessel function somehow.
P.S.: We can determine the PDF of $X$ by using the inversion formula. Note that by the same reasoning as above, we deduce that the characteristic function of $X$ is given by $\varphi_X(t) = J_0(q t),$ a Bessel function of the first kind. It follows from the Inversion Theorem that $X$ admits a probability density function, which is given simply as the Fourier Transform of the characteristic function of $X.$ This has been computed before on this very site (see for example here or here). We have that $$\mathcal{F}(J_0)(\xi) = \chi_{(-1, 1)}(\xi)\frac1\pi \frac{1}{\sqrt{1 - \xi^2}}.$$ However, we are interested in a scaling of this function, so performing a simple substitution in the Fourier Transform integral, we get the expression of the PDF of $X:$ $$f_X(\xi) = \chi_{(-q, q)}(\xi) \frac{1}{\pi} \frac{1}{\sqrt{q^2 - \xi^2}}.$$
A: The arcsine distribution stretched to have support on $[-ghk, ghk]$ has the moments you are looking for.   So, with $-ghk < x < ghk$, the probability density function is $$f(x)=\frac{1}{\pi  \sqrt{{{(ghk)}^{2}}-{{x}^{2}}}}$$ and the cumulative distribution function is $$F(x)=\mathbb P(X\le x) = \frac1\pi \cos^{-1}\left(\frac{-x}{ghk}\right)$$
As a quick check in R of the first twelve moments, setting $ghk=2$ gives the central binomial coefficients:
qstretchedarcsine <- function(prob, ghk){ - ghk * cos(prob * pi) }
thisdata <- qstretchedarcsine(ppoints(10^6), ghk=2)
round(colMeans(outer(thisdata, 1:12, "^")), 12)
#    0   2   0   6   0  20   0  70   0 252   0 924

See https://math.stackexchange.com/a/2629048/6460 for the relationship between the arcsine distribution and Bessel functions
A: Let $f(x)$ be the (unknown) density function. Then the characteristic function is $h(t)=\int_R e^{itx}f(x)dx=\sum_{n=0}^\infty \frac{(it)^nM_n}{n!}$  Then the density function is given by $f(x)=\frac{1}{2\pi}\int_Re^{-itx}h(t)dt$.
Knowing the moments, $h(t)$ can be computed and then $f(x)$ using the inverse Fourier transform.
