Numerical values of Fox H function

I am playing around with the new function from 12.3, FoxH and found the result is not stable. It is interesting that the results are pretty different with the introduction of higher accuracy, and the one under higher accuracy requirement is much more time saving. The first two graphs are under machine precision and the later two are under 30 precision.

Two functions are chosen from a paper "2018On densities of the product, quotient and power of independent subordinators", Theorem 4.1 and Proposition 4.4, the probability density function and the probability distribution function of the inverse $$\beta$$ stable subordinator. In that case, the area under the curve $$f_X(x)$$ with x-axis shall be 1 and and $$F_X(x)$$ shall be never above 1, while the numerical results contradict that. Simply based on the graphs, I feel the one with lower accuracy makes more sense.

The code is pasted below:

Block[{t = 5},
Plot[ParallelTable[
1/t^\[Beta] FoxH[{{}, {{1 - \[Beta], \[Beta]}}}, {{{0, 1}}, {}},
x/t^\[Beta]], {\[Beta], {0.7, 0.8, 0.9}}] // Evaluate, {x, 0,
10}, Filling -> Axis,
AxesLabel -> {"x", "\!$$\*SubscriptBox[\(f$$, $$X$$]\)(x)"}]]
Block[{t = 5},
Plot[ParallelTable[
FoxH[{{{1, 1}}, {{1, \[Beta]}}}, {{{1, 1}}, {{0, 1}}}, x/
t^\[Beta]], {\[Beta], {0.7, 0.8, 0.9}}] // Evaluate, {x, 0, 10},
Filling -> Axis,
AxesLabel -> {"x", "\!$$\*SubscriptBox[\(F$$, $$X$$]\)(x)"}]]
$PreRead = (# /. s_String /; StringMatchQ[s, NumberString] && Precision@ToExpression@s == MachinePrecision :> s <> "50." &); Block[{t = 5}, Plot[ParallelTable[ 1/t^\[Beta] FoxH[{{}, {{1 - \[Beta], \[Beta]}}}, {{{0, 1}}, {}}, x/t^\[Beta]], {\[Beta], {0.7, 0.8, 0.9}}] // Evaluate, {x, 0, 10}, Filling -> Axis, AxesLabel -> {"x", "\!$$\*SubscriptBox[\(f$$, $$X$$]\)(x)"}, WorkingPrecision -> 30]] Block[{t = 5}, Plot[ParallelTable[ FoxH[{{{1, 1}}, {{1, \[Beta]}}}, {{{1, 1}}, {{0, 1}}}, x/ t^\[Beta]], {\[Beta], {0.7, 0.8, 0.9}}] // Evaluate, {x, 0, 10}, Filling -> Axis, AxesLabel -> {"x", "\!$$\*SubscriptBox[\(F$$, $$X$$]\)(x)"}, WorkingPrecision -> 30]]  Looking for your help. • ParallelTable here only parallelizes the construction of the list of function, but not their numerical evaluation, and so is pointless. Jan 6, 2022 at 8:11 • Try: Plot[{FoxH[{{}, {{0.210, 0.810}}}, {{{0, 1}}, {}}, x], FoxH[{{}, {{1/5, 4/5}}}, {{{0, 1}}, {}}, x]}, {x, 0, 3}] Jan 6, 2022 at 8:52 • @Roman Thanks for that. What is the smart way to plot a graph parallelly? Jan 7, 2022 at 5:32 2 Answers It looks like there are two problems: (1) the exact expressions are wrong and (2) the numerical expressions are unstable. Maybe file a bug report with Wolfram? Let's define four functions that should in principle all be the same. First, the exact Fox $$H$$ function (one of the functions you use): f1[x_] = FoxH[{{}, {{1/5, 4/5}}}, {{{0, 1}}, {}}, x];  Next, the numerical version of the same: f2[x_] = FoxH[{{}, {{0.2, 0.8}}}, {{{0, 1}}, {}}, x];  Then the corresponding Mellin–Barnes numerical integral taken at $$\gamma=\frac12$$ (to thread the singularities at $$s=0$$ and $$s=1$$ correctly): J[z_, s_] = z^-s Gamma[s]/Gamma[1/5 + 4s/5]; f3[x_?NumericQ] := 1/(2π) Re[NIntegrate[J[x, 1/2 + I*w], {w, -∞, ∞}]]  Finally, the Fox $$H$$ function expressed in terms of hypergeometric functions: f4[x_] = HypergeometricPFQ[{9/20, 7/10, 19/20}, {2/5, 3/5, 4/5}, -256 x^5/3125]/Gamma[1/5] - x HypergeometricPFQ[{13/20, 9/10, 23/20}, {3/5, 4/5, 6/5}, -256 x^5/3125]/Gamma[-3/5] + x^2 HypergeometricPFQ[{17/20, 11/10, 27/20}, {4/5, 6/5, 7/5}, -256 x^5/3125]/(2 Gamma[-7/5]) - x^3 HypergeometricPFQ[{21/20, 13/10, 31/20}, {6/5, 7/5, 8/5}, -256 x^5/3125]/(6 Gamma[-11/5]);  Note that $$f_3(x)$$ and $$f_4(x)$$ are the "truth" that the other two functions should also give. Several observations: 1. $$f_3(x)$$ and $$f_4(x)$$ are indeed indistinguishable within numerical accuracy. Try Plot[f3[x]-f4[x], {x, 0, 3}] 2. The exact Fox function $$f_1(x)$$ is not evaluated correctly: it diverges from the correct result for $$x>0.7$$. This looks like a bug. 3. The numerical Fox function $$f_2(x)$$ oscillates between $$f_1(x)$$ and $$f_{3,4}(x)$$ and is unusable. 4. The series expansions of $$f_1(x)$$ and $$f_4(x)$$ match at $$x=0$$, which is strange: Series[f1[x]-f4[x], {x, 0, 100}] gives $$O(x)^{101}$$. We see these points in a plot: Plot[{f1[x], f2[x], f3[x], f4[x]}, {x, 0, 3}, PlotLegends -> {"f1", "f2", "f3", "f4"}]  • Might be enlighting to manually$\texttt{NIntegrate}$the$\texttt{FoxH}$integrand over$\mathscr{L}_{-\infty}$,$\mathscr{L}_{\infty}$, and$\mathscr{L}_{\gamma+i\infty}\$ to identify the problem. Designing the contours would be a challenge though.
– josh
Jan 6, 2022 at 13:12
• @josh I've added the contour integral. It's particularly easy in this case because there are only two singularities and we can pass straight between them. Jan 6, 2022 at 18:41
• May I ask what is the difference between adding "?NumericQ" and not adding in the definition of f3? Jan 7, 2022 at 6:11
• It means that f3 can only be called with a numerical argument. Try f3[u] and you'll see that it won't even attempt to integrate numerically. Jan 7, 2022 at 10:03

Many thanks for the help.

I have not identified any errors among these theoretical results. If we apply the Mellin Barnes integral directly, which is f3 in @Roman's answer, the numerical fits the expectation.

I first identified the contradiction in Matheamtica 12.3 for $$\beta=\frac{1}{2}$$. Following the notation in @Roman's answer, consider one example

f1[x_] = 1/(\[Beta] t^(1/\[Beta])) FoxH[{{{1 - 1/\[Beta], 1/\[Beta]}}, {}}, {{}, {{0, 1}}}, x/t^(1/\[Beta])];
p[z_, s_] = Gamma[1/\[Beta] - s/\[Beta]]/Gamma[1 - s] z^-s ;
f3[x_] := 1/(2 \[Pi]) 1/(\[Beta] t^(1/\[Beta]))
Re[NIntegrate[p[x/t^(1/\[Beta]), 1/2 + I*w], {w, -\[Infinity], \[Infinity]}]];


When $$\beta=\frac{1}{2}$$, we know it could be simplified to (Mathematica 13 does not recognise this)

t/(2 Sqrt[\[Pi]]) x^(-(3/2)) E^(-(t^2/(4 x)))


and we can see the f3 fits the theoretical result while f1 does not, see below:

This example somehow verifies Theorem 4.1 in the paper. Proposition 4.4 is easier to be verified since the function shall be non-decreasing, non-negative, bounded below by 1.

Here is the code:

p[z_, s_] = Gamma[s]/Gamma[1/5 + 4 s/5] z^-s;
f[x_] := 1/(2 \[Pi])
Re[NIntegrate[p[x, 1/2 + I*w], {w, -\[Infinity], \[Infinity]}]];
P[z_, s_] = (Gamma[-s] Gamma[1 + s])/(Gamma[1 + 4 s/5] Gamma[1 - s])
z^-s ;
F[x_] := 1/(2 \[Pi])
Re[NIntegrate[P[x, -(1/2) + I*w], {w, -\[Infinity], \[Infinity]}]];
Plot[{f[x], F[x]}, {x, 0, 3}, PlotLegends -> {"PDF", "CDF"}]
`

My current understanding is there are bugs in Fox H function in Mathematica 13.