[Scipy-svn] r3181 - in trunk/Lib: sandbox/cdavid special special/tests

scipy-svn@scip... scipy-svn@scip...
Mon Jul 23 02:33:25 CDT 2007


Author: cdavid
Date: 2007-07-23 02:33:16 -0500 (Mon, 23 Jul 2007)
New Revision: 3181

Added:
   trunk/Lib/special/spfun_stats.py
   trunk/Lib/special/tests/test_spfun_stats.py
Modified:
   trunk/Lib/sandbox/cdavid/
   trunk/Lib/special/__init__.py
Log:
add log multivariate gamma + corresponding tests.


Property changes on: trunk/Lib/sandbox/cdavid
___________________________________________________________________
Name: svn:ignore
   + *.pyc
.*.swp


Modified: trunk/Lib/special/__init__.py
===================================================================
--- trunk/Lib/special/__init__.py	2007-07-23 06:57:45 UTC (rev 3180)
+++ trunk/Lib/special/__init__.py	2007-07-23 07:33:16 UTC (rev 3181)
@@ -11,6 +11,7 @@
 from orthogonal import legendre, chebyt, chebyu, chebyc, chebys, \
      jacobi, laguerre, genlaguerre, hermite, hermitenorm, gegenbauer, \
      sh_legendre, sh_chebyt, sh_chebyu, sh_jacobi, poch
+from spfun_stats import multigammaln
 
 __all__ = filter(lambda s:not s.startswith('_'),dir())
 

Added: trunk/Lib/special/spfun_stats.py
===================================================================
--- trunk/Lib/special/spfun_stats.py	2007-07-23 06:57:45 UTC (rev 3180)
+++ trunk/Lib/special/spfun_stats.py	2007-07-23 07:33:16 UTC (rev 3181)
@@ -0,0 +1,87 @@
+#! /usr/bin/env python
+# Last Change: Mon Jul 23 04:00 PM 2007 J
+
+# Copyright (c) 2001, 2002 Enthought, Inc.
+# 
+# All rights reserved.
+# 
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# 
+#   a. Redistributions of source code must retain the above copyright notice,
+#      this list of conditions and the following disclaimer.
+#   b. Redistributions in binary form must reproduce the above copyright
+#      notice, this list of conditions and the following disclaimer in the
+#      documentation and/or other materials provided with the distribution.
+#   c. Neither the name of the Enthought nor the names of its contributors
+#      may be used to endorse or promote products derived from this software
+#      without specific prior written permission.
+# 
+# 
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
+# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
+# DAMAGE.
+
+"""Some more special functions which may be useful for multivariate statistical
+analysis."""
+
+import numpy as N
+from scipy.special import gammaln as loggam
+
+def multigammaln(a, d):
+    """returns the log of multivariate gamma, also sometimes called the
+    generalized gamma. 
+    
+    :Parameters:
+        a : ndarray
+            the multivariate gamma is computed for each item of a
+        d : int
+            the dimension of the space of integration.
+
+    :Returns:
+        res : ndarray 
+            the values of the log multivariate gamma at the given points a.
+
+    Note
+    ----
+
+    The formal definition of the multivariate gamma of dimension d for a real a
+    is :
+
+    \Gamma_d(a) = \int_{A>0}{e^{-tr(A)\cdot{|A|}^{a - (m+1)/2}dA}}
+
+    with the condition a > (d-1)/2, and A>0 being the set of all the positive
+    definite matrices. Note that a is a scalar: the integration is
+    multivariate, the argument is not.
+
+    This can be proven to be equal to the much friendler equation:
+
+    \Gamma_d(a) = \pi^{d(d-1)/4}\prod_{i=1}^{d}{\Gamma(a - (i-1)/2)}.
+    
+    Reference:
+    ----------
+
+    R. J. Muirhead, Aspects of multivariate statistical theory (Wiley Series in
+    probability and mathematical statistics). """
+    a = N.asarray(a)
+    if not N.isscalar(d) or (N.floor(d) != d):
+        raise ValueError("d should be a positive integer (dimension)")
+    if N.any(a <= 0.5 * (d - 1)):
+        raise ValueError("condition a (%f) > 0.5 * (d-1) (%f) not met" \
+                         % (a, 0.5 * (d-1)))
+
+    res = (d * (d-1) * 0.25) * N.log(N.pi)
+    if a.size == 1:
+        axis = -1
+    else:
+        axis = 0
+    res += N.sum(loggam([(a - (j - 1.)/2) for j in range(1, d+1)]), axis)
+    return res

Added: trunk/Lib/special/tests/test_spfun_stats.py
===================================================================
--- trunk/Lib/special/tests/test_spfun_stats.py	2007-07-23 06:57:45 UTC (rev 3180)
+++ trunk/Lib/special/tests/test_spfun_stats.py	2007-07-23 07:33:16 UTC (rev 3181)
@@ -0,0 +1,39 @@
+import numpy as N
+from numpy.testing import *
+
+set_package_path()
+from spfun_stats import multigammaln
+from special import gammaln
+restore_path()
+
+class test_multigammaln(NumpyTestCase):
+    def test1(self):
+        a = N.abs(N.random.randn())
+        assert_array_equal(multigammaln(a, 1), gammaln(a))
+
+    def test_ararg(self):
+        d = 5
+        a = N.abs(N.random.randn(3, 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))
+
+        d = 5
+        a = N.abs(N.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 test_bararg(self):
+        try:
+            multigammaln(0.5, 1.2)
+            raise Exception("Expected this call to fail")
+        except ValueError:
+            pass
+
+if __name__ == '__main__':
+    NumpyTest().run()



More information about the Scipy-svn mailing list