Fourier transform of hypergeometric function only returns where FT vanishes

Question

I have a probability distribution of a random variable $X$ with support between, roughly, $[-3, 3]$. I do not know the probability distribution, but I can simulate instances from it approximately, and I have its moments. My aim is to find the probability distribution from the moments.

Odd moments vanish, and even moments are $\mathbb{E}[X^{2n}] = \frac{\binom{4n}{n}}{3n+1}$. Using the moments, I constructed the Fourier transform of the probability distribution as $\mathbb{E}[e^{i k X}] = \sum_{n=0}^{\infty} \frac{(-1)^n k^{2n}}{(2n)!} \frac{\binom{4n}{n}}{3n+1}$.

Mathematica gives this sum as HypergeometricPFQ[{1/4, 3/4}, {2/3, 1, 4/3}, -((64 k^2)/27)]

Taking the inverse transformation via

  Exp[-I k x] HypergeometricPFQ[{1/4, 3/4}, {2/3, 1, 4/3}, -((64 k^2)/27)], {k, -Infinity, Infinity}]

gives a curious result: ConditionalExpression[0, 1/x^2 < 27/256].

That is, Mathematica is telling me that the distribution vanishes outside $[-\frac{16}{3 \sqrt{3}}, \frac{16}{3 \sqrt{3}}]$, which is in strong agreement with my empirical bounds of roughly $[-3,3]$. However, there's no information about what the resulting function of $x$ looks like within the bounds of $[-\frac{16}{3 \sqrt{3}}, \frac{16}{3 \sqrt{3}}]$!

How can I evaluate this inverse Fourier transform for all real $x$, not just those where the support vanishes? I would prefer a symbolic solution, but I suspect it's likely that this is the best I can extract from Mathematica symbolically, so I would also accept methods that give the transformation numerically.

Answer 1

f=HypergeometricPFQ[{1/4,3/4},{2/3,1,4/3},-64 k^2/27];

(* this returns a symbolic result *)
g=1/Sqrt[2*Pi]*InverseFourierTransform[f,k,x];

(* check normalization, gives 1 *)
NIntegrate[g,{x,-Infinity,Infinity}]

(* plot *)
Plot[g,{x,-4,4},PlotStyle->Thick]

Answer 2

Amplifying on the answer by user293787

Clear["Global`*"]

f[k_] = Sum[
  (-1)^n k^(2 n)/(2 n)! Binomial[4 n, n]/(3 n + 1),
  {n, 0, Infinity}]

(* HypergeometricPFQ[{1/4, 3/4}, {2/3, 1, 4/3}, -((64 k^2)/27)] *)

(pdf1[x_] = Assuming[Element[x, Reals],
  InverseFourierTransform[f[k], k, x,
    FourierParameters -> {1, -1}] // FullSimplify])//
TraditionalForm

The total probability is unity

NIntegrate[pdf1[x], {x, -Infinity, Infinity}]

(* 1. *)

The pdf is an even function

Assuming[x != 0, pdf1[-x] == pdf1[x] // FullSimplify]

(* True *)

The support is on the interval

NIntegrate[pdf1[x], {x, -16 Sqrt[3]/9, 16 Sqrt[3]/9}]

(* 1. *)

Given this constraint,

(pdf2[x_] = Piecewise[{{
    Assuming[-16 Sqrt[3]/9 < x < 16 Sqrt[3]/9,
     pdf1[x] // FullSimplify],
    -16 Sqrt[3]/9 < x < 16 Sqrt[3]/9}}, 0])//
TraditionalForm

NIntegrate[pdf2[x], {x, -Infinity, Infinity}]

(* 1. *)

Further, due to the even symmetry, the pdf can be simplified to

(pdf3[x_] = Piecewise[{{
    Assuming[0 < x < 16 Sqrt[3]/9,
      pdf1[x] // FullSimplify] /. x -> Sqrt[x^2],
    -16 Sqrt[3]/9 < x < 16 Sqrt[3]/9}}, 0])//
TraditionalForm

NIntegrate[pdf3[x], {x, -Infinity, Infinity}]

(* 1. *)

EDIT: Following through on the comment by user196574

(ratio = (8*3^(1/4)*Gamma[-(1/12)]*Gamma[1/4]^2*Gamma[7/12])/
       (3*3^(3/4)*Gamma[-(1/4)]*Gamma[5/12]*Gamma[3/4]*
          Gamma[13/12])) // N

(* 32. *)

Off[N::meprec]

Block[{$MaxExtraPrecision = 500}, ratio - 32 // N[#, 200] &]

(* 0.*10^-698 *)

(repl = Solve[
      (ratio /. Gamma[7/12] -> Inactive[Gamma][7/12]) == 32, 
      Inactive[Gamma][7/12]][[1]] // Activate) //
 TraditionalForm

pdf4[x_] = pdf3[x] /. repl // FullSimplify

NIntegrate[pdf4[x], {x, -Infinity, Infinity}]

(* 1. *)

The simplification results

LeafCount /@ {pdf1[x], pdf2[x], pdf3[x], pdf4[x]}

(* {300, 218, 172, 146} *)