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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
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

Improve this answer
$\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$ Commented 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} *)
Improve this answer
$\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$ Commented 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$ Commented Oct 15, 2022 at 0:05

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.