Message Boards Message Boards

17 Replies
15 Total Likes
View groups...
Share this post:

A Preliminary Code for the Fox - H Function

Hey guys, how are you doing? I hope really fine!

While doing my research on pure mathematics / Statistics, I came across this beautiful function called Fox-H function.
This function is quite important to the study of Statistics (Algebra of Random Variables) and Science in general (Fractional Partial Differential Equations, for example).
Unfortunately, Mathematica does not have this function implemented. On the other hand, it has everything you need to build a code to implement the function!
It is interesting to notice that Mathematica has the Meijer-G function implemented, which is a special case of the Fox-H function.
In order to compute some results, I did a quick implementation of the function. In the cases I tested (large scale tests with Alfa-Stable random variables of type 1 and their ratio) the code worked nicely.

Please find the code below:
Needs["NumericalCalculus`"];?FoxH[a_, b_, z_] :=? ? Module[{SPA, SPB, IPA, IPB, T, LeftP, RightP, Poles, RadiusP, ?c, ?c, ?   MaxPossibleResidueIncrementsto?, ?, NRightPolesLessThan?, W, H},?  SPA = Product[?    Gamma[1 - a[[1, j, 1]] - a[[1, j, 2]]*s], {j, 1, Length[a[[1]]]}];?  SPB = Product[?    Gamma[b[[1, j, 1]] + b[[1, j, 2]]*s], {j, 1, Length[b[[1]]]}];?  IPA = Product[?    Gamma[a[[2, j, 1]] + a[[2, j, 2]]*s], {j, 1, Length[a[[2]]]}];?  IPB = Product[?    Gamma[1 - b[[2, j, 1]] - b[[2, j, 2]]*s], {j, 1, Length[b[[2]]]}];?  T := SPA*SPB/(IPA*IPB);?  LeftP[p_] := ?   DeleteDuplicates[?    Flatten[?     Table[-(b[[1, j, 1]] + k)/b[[1, j, 2]], {j, 1, Length[b[[1]]]}, {k, 0, ?       p}]]];?  RightP[p_] := ?   DeleteDuplicates[?    Flatten[?     Table[(1 - a[[1, j, 1]] + k)/a[[1, j, 2]], {j, 1, Length[a[[1]]]}, {k, 0,?        p}]]];?  ?c = Product[a[[1, j, 2]]^(-a[[1, j, 2]]), {j, 1, Length[a[[1]]]}]*?    Product[a[[2, j, 2]]^(-a[[2, j, 2]]), {j, 1, Length[a[[2]]]}]*?    Product[b[[1, j, 2]]^(b[[1, j, 2]]), {j, 1, Length[b[[1]]]}]*?    Product[b[[2, j, 2]]^(b[[2, j, 2]]), {j, 1, Length[b[[2]]]}];?  ?c = Sum[b[[1, j, 2]], {j, 1, Length[b[[1]]]}] + ?    Sum[b[[2, j, 2]], {j, 1, ?      Length[b[[2]]]}] - (Sum[a[[1, j, 2]], {j, 1, Length[a[[1]]]}] + ?      Sum[a[[2, j, 2]], {j, 1, Length[a[[2]]]}]);?  Poles[p_] := Sort[DeleteDuplicates[Flatten[{LeftP[p], RightP[p]}]]];?  RadiusP[p_] := ?   Min[Table[?      Abs[Poles[p][[i + 1]] - Poles[p][[i]]], {i, 1, Length[Poles[p]] - 1}]]/2;?  MaxPossibleResidueIncrementsto? = ?   Ceiling[Re[(Max[LeftP[0]] - Min[RightP[0]])*Max[a[[1, All, 2]]]]];?  If[Max[LeftP[0]] < Min[RightP[0]], ? = (Max[LeftP[0]] + Min[RightP[0]])/2, ?   ? = Max[LeftP[0]] + RadiusP[MaxPossibleResidueIncrementsto?]];?  NRightPolesLessThan? = ?   Catch[Do[?     If[Length[Select[RightP[i], # < ? &]] - ?        Length[Select[RightP[i + 1], # < ? &]] >= 0, Throw[i]], {i, 10, 1000, ?      10}]];?  W = Max[Im[Poles[0]]] + 50;?  If[Abs[z] >= 0.2, Which[?c > 0 ? And[?c == 0, 0 < Abs[z] < ?c],?    H[p1_] := ?     Re[(1/(2*Pi*I))*?        NIntegrate[?         T*z^(-s), {s, ? - I*W, ? + I*W, ? - p1 + I*W, ? - p1 - I*W, ?          ? - I*W}]] - ?      Sum[?       Re[?        NResidue[T*z^(-s), {s, r}, ?         Radius -> Min[0.001, RadiusP[MaxPossibleResidueIncrementsto?]]]], {r,?         Select[RightP[NRightPolesLessThan?], # < ? &]}];?    H[1000],?    ?c < 0 ? And[?c == 0, Abs[z] > ?c],?    H[p1_] := ?     Re[(1/(2*Pi*I))*?        NIntegrate[?         T*z^(-s), {s, ? - I*W, ? + I*W, ? + p1 + I*W, ? + p1 - I*W, ?          ? - I*W}]] - ?      Sum[?       Re[?        NResidue[T*z^(-s), {s, r}, ?         Radius -> Min[0.001, RadiusP[MaxPossibleResidueIncrementsto?]]]], {r,?         Select[RightP[NRightPolesLessThan?], # < ? &]}];?    H[1000]], ?   H[p1_] := ?    Re[(1/(2*Pi*I))*?       NIntegrate[T*z^(-s), {s, ? - I*p1, ? + I*p1}, MaxRecursion -> 40, ?        PrecisionGoal -> 15]] - ?     Sum[Re[?       NResidue[T*z^(-s), {s, r}, ?        Radius -> Min[0.001, RadiusP[MaxPossibleResidueIncrementsto?]]]], {r, ?       Select[RightP[NRightPolesLessThan?], # < ? &]}];?   H[2000]]]

The general idea of the code is presented above. I can send the notebook with the code if anyone wants it.

I basically use 3 possible contours in the complex plane to numerically evaluate the complex integral in Mathematica. Each contour is selected according to existence conditions.

The input insertion is similar to that of Meijer-G function. In the case of the H-Function, for example, each of the elements of the sublists of input a is not a constant, but a list with the values {a_j,alpha_j}, according to the definition (

You can take a deeper look at this function's theory on:

For applications and mathematical definitions, one may check:

Mathai, A.M., Saxena, R.K. and Haubold, H.J. (2010) The H-Function: Theory and Applications, Springer, New York.
M. D. Springer (1979), The Algebra of Random Variables, John Wiley, New York.

I have used the code to write a paper about the analytical obtention and evaluation of the PDF of the ratio of two Alfa-Stable Random Variables. The paper has just been submitted for publication but I can discuss it with[font='times new roman', times, serif] anyone interested in the topic. I also have other papers about the usage of the function itself in pure and applied math (analytically solving special real degree equations, civil engineering applications, etc), which I would be also happy to share =)
Anyway, I guess that if you somehow get into this area, this code would be useful.

That is it guys, please let me know if you have any suggestions on the improvement of the code or any ideas on the subject!!
Best Regards

POSTED BY: Luan Ozelim
17 Replies
Posted 3 years ago

Oh that sounds great! Unfortunately, I haven't been working on the Fox H for a while now, forgot this thread even. We found a way to bypass the Fox H for our purposes... however the availability of it now gives some interesting ideas.

Will definitely take a look and perhaps post back at some point again. In any case, thanks for the good hard work!

Sander Paekivi

POSTED BY: Sander Paekivi

Hey Mathematica fans, how are you, I hope fine!

While doing my research on pure mathematics, I had also the problem, that Wolfram does not implemented the Fox-H function. Therefore I searched for code and found some Preliminary Code. I had to do some changes, know it is working. I would like to share it with you. I add some other crazy functions, too.

POSTED BY: Dirk Lehmkuhl
Posted 5 years ago

Dear Luan Ozelim, I am interested in the FoxH functions. But I don't know how to use it. Can you give me an example of how to use it? Is it Hfunction[{{{}}},{{{0,1}}},0.5,10^(-7)]? But it can't work. Thank you very much.


Dear Sander, how are you doing? I hope really fine!

Long time since I last posted a reply here, sorry for that.

I have been using this code and some variations of it since the time I first posted here. So far, everyhting worked well. In fact, my last implementations of the H function are specific to the problems I am trying to solve.

The parts of the code you mention are related to an automated way to separate the poles. MaxPossibleResidueIncrementsto[Gamma] is used to define gamma, a parameter which positions the contour for the numerical integration. Then, NRightPolesLessThan[Gamma] finds the poles which have been incorrectly included in the straight line gamma contour in order to correct the final result.

The code works, in short, as:

1) mount the products;

2) find a way to separate the poles with a straight line;

3) see if any of the right poles have been mistakenly included in the straight line contour;

4) subtract the residues at the poles found in step 3 from the NIntegrate result.

I think this code can be considerably updated (it is almost 5 years old). But for specific applications, a tailored code would be a better choice imo.

Best regards


POSTED BY: Luan Ozelim
Posted 6 years ago


This thread has been created many years ago now, but the Fox H function is only gaining more interest. I'm studying in a PhD program now and recently encountered this function too and am having trouble interpreting this code and using it. Can anyone verify that this code is a proper implementation of the Fox H function? My confusion with it stems from the part of the script, defining increments: "MaxPossibleResidueIncrementsto[Gamma]" and "NRightPolesLessThan[Gamma]"

Any help would be appreciated greatly!

All the best, Sander Paekivi

POSTED BY: Sander Paekivi
Posted 6 years ago



Dear Sir/Madam, I trust you are doing well. Please I am having an issue of plotting the Fox $H$- Function. Please is there any Mathematica code which could help me to plot the following Fox $H$-Function providing that $ \gamma \in (0,1)$ , $s \in (0,1)$ and $d$ is an integer. Thank you so much for any help

Best regards


enter image description here

Dear luan, I would like to ask you how you can use your code to implement a Fox's H function like the fig.Thanks in advanced enter image description here

POSTED BY: hongbing jiang
Posted 6 years ago


Posted 9 years ago

Dear Luan, i use the function to calculate the fox h function,but some error came out,could you help to look into it? i use it to in the statics of image distribution. the main point is i don't know how to input the parameter when the a is null.

thanks in advance


Posted 10 years ago
Dear luanoz Ozelim,
Would you mind sharing with me the notebook file of Fox-H function?
Thanks in advances.
Dear Mr. Thanh Tu,
Thanks for the interest in the code. You may find an updated code for the H-function below.
 corr[list_] :=
   corr[list] = If[Length[list] > 0, {1 - list[[1]], list[[2]]}, {}];
 SPA[a_] :=
   SPA[a] = Product[
     Gamma[1 - a[[1, j, 1]] - a[[1, j, 2]]*s], {j, 1, Length[a[[1]]]}];
 SPB[b_] :=
   SPB[b] = Product[
    Gamma[b[[1, j, 1]] + b[[1, j, 2]]*s], {j, 1, Length[b[[1]]]}];
IPA[a_] :=
  IPA[a] = Product[
    Gamma[a[[2, j, 1]] + a[[2, j, 2]]*s], {j, 1, Length[a[[2]]]}];
IPB[b_] :=
  IPB[b] = Product[
    Gamma[1 - b[[2, j, 1]] - b[[2, j, 2]]*s], {j, 1, Length[b[[2]]]}];
T[a_, b_, s_] := T[a, b, s] = SPA[a]*SPB[b]/(IPA[a]*IPB[b]);
LeftP[b_, p_] :=
  LeftP[b, p] =
    Flatten[Table[-(b[[1, j, 1]] + k)/b[[1, j, 2]], {j, 1,
       Length[b[[1]]]}, {k, 0, p}]]];
RightP[a_, p_] :=
  RightP[a, p] =
    Flatten[Table[(1 - a[[1, j, 1]] + k)/a[[1, j, 2]], {j, 1,
       Length[a[[1]]]}, {k, 0, p}]]];
Poles[a_, b_, p_] :=
  Poles[a, b, p] =
   Sort[DeleteDuplicates[Flatten[{LeftP[b, p], RightP[a, p]}]]];
RadiusP[a_, b_, p_] :=
  RadiusP[a, b, p] =
      Abs[Poles[a, b, p][[i + 1]] - Poles[a, b, p][[i]]], {i, 1,
       Length[Poles[a, b, p]] - 1}]]/2;
MaxPossibleResidueIncrementsto\[Gamma][a_, b_] :=
  MaxPossibleResidueIncrementsto\[Gamma][a, b] =
   Ceiling[Re[(Max[LeftP[b, 0]] - Min[RightP[a, 0]])*
      Max[a[[1, All, 2]]]]];
\[Gamma][a_, b_] := \[Gamma][a, b] =
   If[Max[LeftP[b, 0]] <
     Min[RightP[a, 0]], (Max[LeftP[b, 0]] + Min[RightP[a, 0]])/2,
    Max[LeftP[b, 0]] +
     RadiusP[a, b, MaxPossibleResidueIncrementsto\[Gamma][a, b]]];
NRightPolesLessThan\[Gamma][a_, b_] :=
  NRightPolesLessThan\[Gamma][a, b] =
     If[Length[Select[RightP[a, i], # < \[Gamma][a, b] &]] -
        Length[Select[RightP[a, i + 1], # < \[Gamma][a, b] &]] >= 0,
      Throw[i]], {i, 10, 1000, 10}]];
H[a_, b_, z_, p1_, pg_] :=
     T[a, b, s]*
      z^(-s), {s, \[Gamma][a, b] - I*p1, \[Gamma][a, b] + I*p1},
     MaxRecursion -> 20, PrecisionGoal -> pg]] -
  Sum[Re[NResidue[T[a, b, s]*z^(-s), {s, r},
     Radius ->
       RadiusP[a, b,
        MaxPossibleResidueIncrementsto\[Gamma][a, b]]]]], {r,
      NRightPolesLessThan\[Gamma][a, b]], # < \[Gamma][a, b] &]}]

Hfunction[as_, bs_, z1_, pg_] :=
Module[{W, a, b, z},
  If[Abs[z1] > 1, z = z1; a = as; b = bs, z = 1/z1;
   a = Map[corr, bs, {2}];
   b = Map[corr, as, {2}]];
  W = Max[Im[Poles[a, b, 0]]] + 50;
  H[a, b, z, 1000, pg]]
You may also find the notebook attached to this post. The function has been now named Hfunction instead of FoxH. Also, now you may set the precision goal by changing the last parameter of the function.
Since in my case I used the same parameters a lot, the f := f = structure has been used. You may change that if you wish.
Best Regards

POSTED BY: Luan Ozelim

Dear luanoz Ozelim, When I tried to open the attached file in the previous reply, it give an error. would you please email me the code that generate the plots. My email is Thanks ~Moh

POSTED BY: mohammed aloqlah
Posted 11 years ago
Dear Mr.  luanoz Ozelim, 
I would like to ask you how you can use your code to implement any-function, (e.g. H_{1,1}^{1,1}=[5|{1,2};{2,5}]), because I have a paper include the a code to implement H function on Mathematica, but unfortunatelly I can not use it in Mathematica
Thanks in advanced
POSTED BY: Hosam rh
Hello Mrs. Hosam rh.
The structure you are using is not correct for this implementation of the H function. Which equation in the paper are you trying to implement?
Let me know so that I can suggest a proper structure!
Best Regards,

POSTED BY: Luan Ozelim
Hi Vitaliy,

Please find below a brief comparison between the built-in PDF of StableDistributions of type 1 and the correspondend expression in H-function. This shows that the algorithm works nicely in these cases. Any other thing you need, just let me know!

Assuming the function from my original post is defined:
Levy?StableDistribution[x_, ?_, ?_, ?_, ?_] :=?  Module[{?1, ?1, K, F},?   ?1 = ?*(1 + (?^2)*(Tan[?*Pi/2]^2))^(1/(2*?));?   K = Piecewise[{{?, ? < 1}, {? - 2, ? > 1}}];?   ?1 = Piecewise[{{(2/(Pi*?))*ArcTan[?*Tan[?*Pi/2]], ?       0 < ? < 1}, {(2/(Pi*(? - 2)))*ArcTan[?*Tan[(? - 2)*Pi/2]], 1 < ? < 2}}];?   F = Piecewise[{{(1/(?*?1))*?        FoxH[{{{1 - 1/?, 1/?}}, {{1/2 - ?1*K/(2*?), 1/2 + ?1*K/(2*?)}}}, {{{0,?             1}}, {{1/2 - ?1*K/(2*?), 1/2 + ?1*K/(2*?)}}}, (x - ?)/?1], ?       x >= ?}, {(1/(?*?1))*?        FoxH[{{{1 - 1/?, 1/?}}, {{1/2 + ?1*K/(2*?), 1/2 - ?1*K/(2*?)}}}, {{{0,?             1}}, {{1/2 + ?1*K/(2*?), 1/2 - ?1*K/(2*?)}}}, (? - x)/?1], ?       x < ?}}];?   F];
??Tabe = {{0.6, 0.4, 0.5, 2}, {1.5, -0.7, -1, 5}};??
dadosdst3 = ?  Quiet[Table[{x, Levy?StableDistribution[x, Sequence @@ Tabe[[i]]]}, {i, 1, ?     2}, {x, -20.11, 20.11, 0.2}]];??
GraphicsRow[? Table[Show[{ListPlot[dadosdst3[[i]], PlotRange -> Full, PlotStyle -> Red, ?     ImageSize -> 300, PlotLegends -> {"H-function Evaluation"}], ?    Plot[PDF[StableDistribution[1, Sequence @@ Tabe[[i]]], x], {x, -20.2, ?      20.2}, PlotRange -> Full, PlotStyle -> Black, ImageSize -> 300, ?     PlotLegends -> {"Built-in StableDistribution PDF"}]}], {i, 1, 2}], ? ImageSize -> {1000}]
POSTED BY: Luan Ozelim
Do you have example of usage for your function? Something simple - like a plot ?
POSTED BY: Vitaliy Kaurov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract