Added multigamma function #16843

ritesh99rakesh · 2019-05-16T13:22:59Z

References to other Issues or PRs

This is a new PR which was created because
I made a blunder with old PR: #16822

Brief description of what is fixed or changed

Multivariate Gamma function is added in sympy/functions/special/gamma_functions.py
Need to discuss with mentors if this is needed or not. Will add test after finalization.

Other comments

Multivariate Gamma function is a generalization of Gamma function and is useful in
multivariate statistics appearing in the probability density function of the Wishart and inverse Wishart distributions, and the matrix variate beta distribution.

Release Notes

functions
- Added multigamma function in sympy/functions/special/gamma_functions.py

sympy-bot · 2019-05-16T13:23:00Z

✅

Hi, I am the SymPy bot (v147). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

functions
- Added multigamma function in sympy/functions/special/gamma_functions.py (#16843 by @ritesh99rakesh)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.5.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234". See
https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
write a comment on that issue linking back to this pull request once it is
open. -->
This is a new PR which was created because
I made a blunder with old PR: https://github.com/sympy/sympy/pull/16822

#### Brief description of what is fixed or changed
Multivariate Gamma function is added in `sympy/functions/special/gamma_functions.py`
Need to discuss with mentors if this is needed or not. Will add test after finalization.

#### Other comments
`Multivariate Gamma function` is a generalization of `Gamma function` and is useful in 
multivariate statistics appearing in the probability density function of the `Wishart` and `inverse Wishart` distributions, and the `matrix variate beta` distribution.

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- functions
    - Added `multigamma` function in `sympy/functions/special/gamma_functions.py`
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

ritesh99rakesh · 2019-05-16T13:25:37Z

@sidhantnagpal this is a new PR since I did a blunder with older PR. Please review this.

sidhantnagpal · 2019-05-16T14:34:49Z

sympy/functions/special/gamma_functions.py

+        raise TypeError('p must be a positive integer')
+    from sympy import Product
+    from sympy.abc import k
+    return pi**(psym*(psym-1)/4)*Product(gamma(x + (1 - k)/2), (k, 1, psym))


Lower case methods return evaluated results when possible, that is they shouldn't require doit.
That said, unevaluated expressions can still be obtained by passing evaluate=False.

>>> gamma(5) 24 >>> gamma(S(5)/2) 3*sqrt(pi)/4 >>> gamma(7, evaluate=False) gamma(7)

What do think about the following diff:

--- a/sympy/functions/special/gamma_functions.py +++ b/sympy/functions/special/gamma_functions.py @@ -1010,7 +1010,7 @@ def trigamma(x): """ return polygamma(1, x) -def multivariate_gamma(x, p): +def multivariate_gamma(x, p, evaluate=True): r""" The multivariate gamma function is a generalization of the gamma function i.e, @@ -1065,8 +1065,11 @@ def multivariate_gamma(x, p): .. [1] https://en.wikipedia.org/wiki/Multivariate_gamma_function """ psym = sympify(p) - if fuzzy_not(fuzzy_and([psym.is_integer, psym.is_positive])): + if psym.is_integer is False or psym.is_positive is False: raise TypeError('p must be a positive integer') from sympy import Product from sympy.abc import k + if evaluate: + return (pi**(psym*(psym-1)/4)*Product(gamma(x + (1 - k)/2), + (k, 1, psym))).doit() return pi**(psym*(psym-1)/4)*Product(gamma(x + (1 - k)/2), (k, 1, psym))

Example usage:

>>> from sympy import * >>> x, p = symbols('x p') >>> mg = multivariate_gamma >>> mg(1, 1) 1 >>> mg(1, 1, evaluate=False) Product(gamma(3/2 - k/2), (k, 1, 1)) >>> mg(x, 1) gamma(x) >>> mg(x, 1, evaluate=False) Product(gamma(-k/2 + x + 1/2), (k, 1, 1)) >>> mg(x, p) pi**(p*(p - 1)/4)*Product(gamma(-k/2 + x + 1/2), (k, 1, p)) >>> mg(1, p)

I think multivariate_gamma could be implemented along the lines of gamma inheriting from Function (and not digamma or trigamma which are just Python functions based off polygamma).

@sidhantnagpal I am trying to implement multivariate gamma function along the lines of gamma class. I found the following bug with gamma function.

>>> gamma(-1) # -1 is pole for gamma function zoo >>> gamma(-1.0) Traceback (most recent call last): File "/home/ritesh/Desktop/sympy/sympy/core/evalf.py", line 1308, in evalf rf = evalf_table[x.func] KeyError: gamma During handling of the above exception, another exception occurred: Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/home/ritesh/Desktop/sympy/sympy/core/cache.py", line 94, in wrapper retval = cfunc(*args, **kwargs) File "/home/ritesh/Desktop/sympy/sympy/core/function.py", line 463, in __new__ result = result.evalf(mlib.libmpf.prec_to_dps(pr)) File "/home/ritesh/Desktop/sympy/sympy/core/evalf.py", line 1442, in evalf result = evalf(self, prec + 4, options) File "/home/ritesh/Desktop/sympy/sympy/core/evalf.py", line 1314, in evalf xe = x._eval_evalf(prec) File "/home/ritesh/Desktop/sympy/sympy/core/function.py", line 586, in _eval_evalf v = func(*args) File "/home/ritesh/venv/lib/python3.6/site-packages/mpmath/ctx_mp_python.py", line 1012, in f return ctx.make_mpf(mpf_f(x._mpf_, prec, rounding)) File "/home/ritesh/venv/lib/python3.6/site-packages/mpmath/libmp/gammazeta.py", line 1954, in mpf_gamma raise ValueError("gamma function pole") ValueError: gamma function pole

I have opened this issue for this. After this is sorted, we can multivariate_gamma will work even for negative numbers. Meanwhile, I will share my code here, for you to review.

sympy/functions/special/gamma_functions.py

ritesh99rakesh · 2019-05-17T07:33:18Z

@oscargus this PR is related to functions module. Please remove stats and add functions label.

czgdp1807 · 2019-05-18T03:24:37Z

@ritesh99rakesh Please explain the issues which you faced and compelled you to close the old PR.

czgdp1807 · 2019-05-18T03:27:19Z

sympy/functions/special/gamma_functions.py

+    if fuzzy_not(fuzzy_and([psym.is_integer, psym.is_positive])):
+        raise TypeError('p must be a positive integer')
+    from sympy import Product
+    from sympy.abc import k


Please shift the above two imports to the top. It would be great if you can define your symbol k under appropriate assumptions rather than importing it from sympy.abc because it is slower.

czgdp1807 · 2019-05-18T03:28:41Z

sympy/functions/special/gamma_functions.py

+
+    .. [1] https://en.wikipedia.org/wiki/Multivariate_gamma_function
+    """
+    psym = sympify(p)


Would it be possible to use p = sympify(p)?

czgdp1807 · 2019-05-18T03:53:07Z

sympy/functions/special/tests/test_gamma_functions.py

+    assert multivariate_gamma(n, 1).rewrite(factorial).doit() == factorial(n - 1)
+    assert multivariate_gamma(n, 2).rewrite(factorial).doit() == sqrt(pi)*\
+        factorial(n - S(3)/2)*factorial(n - 1)
+    assert multivariate_gamma(n, 3).rewrite(factorial).doit() == pi**(S(3)/2)*\


So many doit() s? Please don't do it. See how, integrate function has been defined and go according to the same lines. @sidhantnagpal What do you say?

#16843 (comment)

ritesh99rakesh · 2019-05-30T18:10:42Z

@sidhantnagpal I have replaced the multivariate_gamma function with multivariate_gamma class. Please review it. Also, I am having trouble with the conjugate part. Since p is always real and positive,

>>> from sympy import symbols, conjugate, multivariate_gamma as mg
>>> x, p = symbols('x p')

>>> conjugate(mg(x, p)) # expected result
pi**((p - 1)*p/4)*Product(gamma(-_k/2 + conjugate(x) + 1/2), (_k, 1, p))

>>> conjugate(mg(x, p)) # present result
pi**((conjugate(p) - 1)*conjugate(p)/4)*Product(gamma(-_k/2 + conjugate(x) + 1/2), (_k, 1, p))

How can I explicitly impose the assumption that p is real and positive.

oscarbenjamin · 2019-05-30T19:14:02Z

How can I explicitly impose the assumption that p is real and positive.

p = Symbol('p', positive=True) # positive implies real

sympy/functions/special/gamma_functions.py

ritesh99rakesh · 2019-05-31T17:19:32Z

@sidhantnagpal I will make the suggested changes and also correct the tests. I wanted to know your view on the class implementation. What additions should I do.

sidhantnagpal · 2019-06-01T04:48:02Z

The class implementation looks OK to me, consistent with other classes for gamma functions. There could be another method _eval_evalf (not necessarily to be added in this PR) to allow the user to work with specified precision for multivariate_gamma.

ritesh99rakesh · 2019-06-01T14:19:08Z

@sidhantnagpal Please review.

sympy/functions/special/gamma_functions.py

sympy/functions/special/tests/test_gamma_functions.py

sympy/functions/special/gamma_functions.py

ritesh99rakesh · 2019-06-02T16:20:42Z

@sidhantnagpal Please review this and if this looks good to you, please consider merging. After this PR gets merged, I will be able to implement various multivariate distributions.

sympy/functions/special/gamma_functions.py

ritesh99rakesh · 2019-06-07T06:07:17Z

@sidhantnagpal Please review this PR.

ritesh99rakesh · 2019-06-13T16:37:45Z

@sidhantnagpal I have made the changes from the above suggestions. Please review.

I can increase the test coverage in a new PR.

ritesh99rakesh · 2019-06-15T09:05:28Z

ping @sidhantnagpal

sidhantnagpal · 2019-06-15T13:21:09Z

sympy/functions/special/gamma_functions.py

+    def _eval_is_real(self):
+        x, p = self.args
+        y = 2*x
+        if isinstance(x, Symbol) or isinstance(p, Symbol):


Why restrict to symbols? It should generalise to expressions as well.

@sidhantnagpal what do you think about the following:

if x.has(Symbol) or p.has(Symbol)

If there is another better method please let me know.

Is it possible to just drop these checks? (the code could be made generic to symbolic/numeric expressions)

ritesh99rakesh · 2019-06-16T12:57:28Z

@sidhantnagpal please review.

sidhantnagpal

Looks good. Minor comments.

sidhantnagpal · 2019-06-18T05:47:27Z

sympy/functions/special/tests/test_gamma_functions.py

+    from sympy import Product
+    p = Symbol('p')
+
+    assert str(multigamma(x, p)) == ('pi**(p*(p - 1)/4)*Product(gamma(-_k/2 + x + 1/2), '


You can avoid string comparison by making use of dummy_eq (here and in other places).

sympy/functions/special/gamma_functions.py

Upabjojr · 2019-06-18T08:38:14Z

sympy/functions/special/gamma_functions.py

@@ -137,8 +136,7 @@ def eval(cls, arg):

                    if arg.is_positive:
                        return coeff*sqrt(S.Pi) / 2**n
-                    else:
-                        return 2**n*sqrt(S.Pi) / coeff
+                    return 2**n*sqrt(S.Pi) / coeff


unnecessary change

sympy/functions/special/gamma_functions.py

Upabjojr · 2019-06-18T08:39:17Z

sympy/functions/special/gamma_functions.py

@@ -659,8 +659,7 @@ def eval(cls, n, z):
                        if n.is_Number:
                            if n is S.Zero:
                                return S.Infinity
-                            else:
-                                return S.Zero
+                            return S.Zero


unnecessary change

Upabjojr · 2019-06-18T08:39:34Z

sympy/functions/special/gamma_functions.py

@@ -739,24 +737,21 @@ def _eval_expand_func(self, **hints):
    def _eval_rewrite_as_zeta(self, n, z, **kwargs):
        if n >= S.One:
            return (-1)**(n + 1)*factorial(n)*zeta(n + 1, z)
-        else:
-            return self
+        return self


unnecessary change

sympy/functions/special/gamma_functions.py

Upabjojr · 2019-06-18T08:39:46Z

sympy/functions/special/gamma_functions.py


    def _eval_as_leading_term(self, x):
        from sympy import Order
        n, z = [a.as_leading_term(x) for a in self.args]
        o = Order(z, x)
        if n == 0 and o.contains(1/x):
            return o.getn() * log(x)
-        else:
-            return self.func(n, z)
+        return self.func(n, z)


unnecessary change

sympy/functions/special/gamma_functions.py

sympy/functions/special/tests/test_gamma_functions.py

…or .is_real in multigamma and gamma

ritesh99rakesh · 2019-06-28T16:26:35Z

@sidhantnagpal I have made the changes. Please review.

sidhantnagpal · 2019-06-28T16:51:01Z

Not all reviews are addressed. See #16843 (comment), #16843 (comment).

Upabjojr · 2019-06-28T19:11:37Z

Yes, there are still some open reviews.

codecov · 2019-07-05T15:42:43Z

Codecov Report

Merging #16843 into master will decrease coverage by 0.007%.
The diff coverage is 64.864%.

@@             Coverage Diff              @@
##            master   #16843       +/-   ##
============================================
- Coverage   74.507%    74.5%   -0.008%     
============================================
  Files          623      623               
  Lines       161379   161412       +33     
  Branches     37875    37882        +7     
============================================
+ Hits        120240   120252       +12     
- Misses       35809    35829       +20     
- Partials      5330     5331        +1

ritesh99rakesh · 2019-07-05T15:45:24Z

@sidhantnagpal done. Please review.

sidhantnagpal · 2019-07-06T04:48:36Z

sympy/functions/special/tests/test_gamma_functions.py

+    p = Symbol('p')
+    _k = Dummy('_k')
+
+    assert multigamma(x, p).dummy_eq(pi**(p*(p - 1)/4)*\


Continued line markers \ are not required when enclosed in braces.

sidhantnagpal · 2019-07-06T04:52:16Z

I left a comment for future PRs. This PR looks ready to me.
I updated the title and release notes entry to reflect the function name multigamma.

sidhantnagpal reviewed May 16, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

oscargus added GSoC stats labels May 16, 2019

czgdp1807 reviewed May 18, 2019

View reviewed changes

sidhantnagpal reviewed May 31, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

sidhantnagpal reviewed May 31, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

oscargus added functions and removed stats labels May 31, 2019

czgdp1807 reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Show resolved Hide resolved

czgdp1807 reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

czgdp1807 reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Show resolved Hide resolved

czgdp1807 reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Show resolved Hide resolved

czgdp1807 reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/tests/test_gamma_functions.py Outdated Show resolved Hide resolved

sidhantnagpal reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

sidhantnagpal reviewed Jun 1, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

oscarbenjamin reviewed Jun 2, 2019

View reviewed changes

sympy/functions/special/gamma_functions.py Outdated Show resolved Hide resolved

sidhantnagpal reviewed Jun 15, 2019

View reviewed changes

sidhantnagpal reviewed Jun 18, 2019

View reviewed changes

Upabjojr reviewed Jun 18, 2019

View reviewed changes

ritesh99rakesh added 11 commits June 28, 2019 21:36

Added multivariate_gamma function

f55eb80

Added evaluate parameter

87bb67e

Replaced multivariate_gamma func with multivariate_gamma class

8c31631

Added check conditions

f96d9d0

Minor changes

f2e477e

Removed replacing of symbol

d66f46e

Changed name from multivariate_gamma to multigamma, corrected logic f…

4e6c9c6

…or .is_real in multigamma and gamma

Corrected tests

a233708

Corrected is_real function

c7da9fd

Corrected is_real function

c671655

Made tests generic

d877880

ritesh99rakesh force-pushed the multivariateGamma branch from 17202af to d877880 Compare June 28, 2019 16:25

Added multivariate gamma function

b9061ba

Removed unecessary change, used dummy_eq

2eeb7a7

ritesh99rakesh force-pushed the multivariateGamma branch from 398fcd0 to 2eeb7a7 Compare July 5, 2019 15:44

sidhantnagpal reviewed Jul 6, 2019

View reviewed changes

sidhantnagpal changed the title ~~Added multivariate_gamma function~~ Added multigamma function Jul 6, 2019

sidhantnagpal merged commit ff7492f into sympy:master Jul 6, 2019

Added multigamma function #16843

Added multigamma function #16843

Conversation

ritesh99rakesh commented May 16, 2019 • edited by sidhantnagpal

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

sympy-bot commented May 16, 2019 • edited

ritesh99rakesh commented May 16, 2019

Choose a reason for hiding this comment

ritesh99rakesh May 16, 2019 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ritesh99rakesh commented May 17, 2019

czgdp1807 commented May 18, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ritesh99rakesh commented May 30, 2019

oscarbenjamin commented May 30, 2019

ritesh99rakesh commented May 31, 2019

sidhantnagpal commented Jun 1, 2019

ritesh99rakesh commented Jun 1, 2019

ritesh99rakesh commented Jun 2, 2019

ritesh99rakesh commented Jun 7, 2019

ritesh99rakesh commented Jun 13, 2019

ritesh99rakesh commented Jun 15, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ritesh99rakesh commented Jun 16, 2019

sidhantnagpal left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ritesh99rakesh commented Jun 28, 2019

sidhantnagpal commented Jun 28, 2019

Upabjojr commented Jun 28, 2019

codecov bot commented Jul 5, 2019 • edited

Codecov Report

ritesh99rakesh commented Jul 5, 2019

Choose a reason for hiding this comment

sidhantnagpal commented Jul 6, 2019

ritesh99rakesh commented May 16, 2019 •

edited by sidhantnagpal

sympy-bot commented May 16, 2019 •

edited

ritesh99rakesh May 16, 2019 •

edited

codecov bot commented Jul 5, 2019 •

edited