[Scipy-svn] r3181 - in trunk/Lib: sandbox/cdavid special special/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Jul 23 03:33:25 EDT 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