3
$\begingroup$

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.

$\endgroup$

2 Answers 2

7
$\begingroup$
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]

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Really nice to see! Curious that InverseFourierTransform fares better than direct integration, but it likely makes sense, since the former can be more specialized. This solution also helped me realize a curious divergence structure at small arguments; I'd been expecting $x^{-2/3}$, but this diverges as $x^{-1/2}$, confirmed in my approximate sampling histograms. $\endgroup$
    – user196574
    Oct 14, 2022 at 18:02
3
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

(* 1. *)

The simplification results

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

(* {300, 218, 172, 146} *)
$\endgroup$
2
  • $\begingroup$ Smart simplifications! I was working out removing the Heaviside thetas, but you've gotten a much crisper result than I was getting! Nice trick also with writing $x$ as $\sqrt{x^2}$ to make the evenness manifest. $\endgroup$
    – user196574
    Oct 14, 2022 at 23:41
  • $\begingroup$ I think there's even one more (small) simplification, which is that the numerical prefactors of the first two terms are the same up to a factor of $32$, that is, N[(8 3^(1/4) Gamma[-1/12] Gamma[1/4]^2 Gamma[7/12])/( 3^(7/4) Gamma[-1/4] Gamma[5/12] Gamma[3/4] Gamma[13/12])] is 32. $\endgroup$
    – user196574
    Oct 15, 2022 at 0:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.