Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

BUG: special: multigammaln returned wrong size array (closes ticket #1849) #455

Merged
merged 1 commit into from

2 participants

@WarrenWeckesser
Collaborator

The only change to the function is to remove a few lines that set axis=-1 if the size of the input array is 1. I'm guessing there was once a good reason for this, but later changes made this incorrect.

I added a few more tests, and refactored the array result tests to use a test generator.

The ticket on trac is http://projects.scipy.org/scipy/ticket/1849

@pv pv merged commit e5b5de4 into scipy:master
@pv
Owner
pv commented

Thanks, it's correct.

@ClemensFMN ClemensFMN referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
@ClemensFMN ClemensFMN referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
7 scipy/special/spfun_stats.py
@@ -88,10 +88,5 @@ def multigammaln(a, d):
% (a, 0.5 * (d-1)))
res = (d * (d-1) * 0.25) * np.log(np.pi)
- if a.size == 1:
- axis = -1
- else:
- axis = 0
-
- res += np.sum(loggam([(a - (j - 1.)/2) for j in range(1, d+1)]), axis)
+ res += np.sum(loggam([(a - (j - 1.)/2) for j in range(1, d+1)]), axis=0)
return res
View
65 scipy/special/tests/test_spfun_stats.py
@@ -2,42 +2,63 @@
import numpy as np
from numpy.testing import assert_array_equal, TestCase, run_module_suite, \
- assert_array_almost_equal_nulp
+ assert_array_almost_equal_nulp, assert_raises, assert_almost_equal
from scipy.special import gammaln, multigammaln
class TestMultiGammaLn(TestCase):
+
def test1(self):
+ # A test of the identity
+ # Gamma_1(a) = Gamma(a)
np.random.seed(1234)
a = np.abs(np.random.randn())
assert_array_equal(multigammaln(a, 1), gammaln(a))
- def test_ararg(self):
- d = 5
- np.random.seed(1234)
- a = np.abs(np.random.randn(3, 2)) + d
+ def test2(self):
+ # A test of the identity
+ # Gamma_2(a) = sqrt(pi) * Gamma(a) * Gamma(a - 0.5)
+ a = np.array([2.5, 10.0])
+ result = multigammaln(a, 2)
+ expected = np.log(np.sqrt(np.pi)) + gammaln(a) + gammaln(a - 0.5)
+ assert_almost_equal(result, expected)
- tr = multigammaln(a, d)
- assert_array_equal(tr.shape, a.shape)
- for i in range(a.size):
- assert_array_almost_equal_nulp(tr.ravel()[i],
- multigammaln(a.ravel()[i], d))
+ def test_bararg(self):
+ assert_raises(ValueError, multigammaln, 0.5, 1.2)
- d = 5
- a = np.abs(np.random.randn(1, 2)) + d
- tr = multigammaln(a, d)
- assert_array_equal(tr.shape, a.shape)
- for i in range(a.size):
- assert_array_equal(tr.ravel()[i], multigammaln(a.ravel()[i], d))
+def _check_multigammaln_array_result(a, d):
+ # Test that the shape of the array returned by multigammaln
+ # matches the input shape, and that all the values match
+ # the value computed when multigammaln is called with a scalar.
+ result = multigammaln(a, d)
+ assert_array_equal(a.shape, result.shape)
+ a1 = a.ravel()
+ result1 = result.ravel()
+ for i in range(a.size):
+ assert_array_almost_equal_nulp(result1[i], multigammaln(a1[i], d))
- def test_bararg(self):
- try:
- multigammaln(0.5, 1.2)
- raise Exception("Expected this call to fail")
- except ValueError:
- pass
+
+def test_multigammaln_array_arg():
+ # Check that the array returned by multigammaln has the correct
+ # shape and contains the correct values. The cases have arrays
+ # with several differnent shapes.
+ # The cases include a regression test for ticket #1849
+ # (a = np.array([2.0]), an array with a single element).
+ np.random.seed(1234)
+
+ cases = [
+ # a, d
+ (np.abs(np.random.randn(3, 2)) + 5, 5),
+ (np.abs(np.random.randn(1, 2)) + 5, 5),
+ (np.arange(10.0, 18.0).reshape(2, 2, 2), 3),
+ (np.array([2.0]), 3),
+ (np.float64(2.0), 3),
+ ]
+
+ for a, d in cases:
+ yield _check_multigammaln_array_result, a, d
if __name__ == '__main__':
Something went wrong with that request. Please try again.