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 (http://www.wolframalpha.com/input/?i=fox+H+function).
You can take a deeper look at this function's theory on:
http://www.wolframalpha.com/input/?i=fox+H+functionhttp://en.wikipedia.org/wiki/Fox_H-functionFor 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
Luan
