[Scipy-svn] r2280 - in trunk/Lib/sandbox/pyem: . pyem pyem/tests

scipy-svn at scipy.org scipy-svn at scipy.org
Thu Oct 12 08:48:05 CDT 2006


Author: cdavid
Date: 2006-10-12 08:47:57 -0500 (Thu, 12 Oct 2006)
New Revision: 2280

Added:
   trunk/Lib/sandbox/pyem/pyem/tests/
   trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py
Modified:
   trunk/Lib/sandbox/pyem/.bzrignore
   trunk/Lib/sandbox/pyem/pyem/densities.py
   trunk/Lib/sandbox/pyem/pyem/gauss_mix.py
   trunk/Lib/sandbox/pyem/pyem/gmm_em.py
   trunk/Lib/sandbox/pyem/pyem/info.py
Log:
[pyem @ david at ar.media.kyoto-u.ac.jp-20060911092105-c3e0d6a0da79b0ca]
- Improve plot1d method to be more useful

- create tests dir (do not do anything yet)

- update doc
David Cournapeau <david at ar.media.kyoto-u.ac.jp> | 2006-09-11 18:21:05 +0900 (Mon, 11 Sep 2006)

Modified: trunk/Lib/sandbox/pyem/.bzrignore
===================================================================
--- trunk/Lib/sandbox/pyem/.bzrignore	2006-10-12 13:47:46 UTC (rev 2279)
+++ trunk/Lib/sandbox/pyem/.bzrignore	2006-10-12 13:47:57 UTC (rev 2280)
@@ -22,3 +22,4 @@
 matcode
 ../MSG
 MSG
+exinfo.py

Modified: trunk/Lib/sandbox/pyem/pyem/densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/densities.py	2006-10-12 13:47:46 UTC (rev 2279)
+++ trunk/Lib/sandbox/pyem/pyem/densities.py	2006-10-12 13:47:57 UTC (rev 2280)
@@ -1,7 +1,7 @@
 #! /usr/bin/python
 #
 # Copyrighted David Cournapeau
-# Last Change: Tue Aug 29 04:00 PM 2006 J
+# Last Change: Wed Aug 30 12:00 PM 2006 J
 
 import numpy as N
 import numpy.linalg as lin
@@ -88,8 +88,6 @@
     of gaussian pdf in scalar case. It assumes all args
     are conformant, so it should not be used directly
     
-    ** Expect centered data (ie with mean removed) **
-
     Call gauss_den instead"""
     d       = mu.size
     inva    = 1/va
@@ -113,8 +111,6 @@
     of gaussian pdf in scalar case. It assumes all args
     are conformant, so it should not be used directly
     
-    ** Expect centered data (ie with mean removed) **
-
     Call gauss_den instead"""
     # Diagonal matrix case
     d   = mu.size
@@ -151,8 +147,6 @@
     It assumes all args are conformant, so it should 
     not be used directly Call gauss_den instead
     
-    ** Expect centered data (ie with mean removed) **
-
     Does not check if va is definite positive (on inversible 
     for that matter), so the inverse computation and/or determinant
     would throw an exception."""
@@ -234,153 +228,6 @@
 
     return elps[0, :], elps[1, :]
 
-def generate_test_data(n, d, mode = 'diag', file='test.dat'):
-    """Generate a set of data of dimension d, with n frames,
-    that is input data, mean, var and output of gden, so that
-    other implementations can be tested against"""
-    mu  = randn(1, d)
-    if mode == 'diag':
-        va  = abs(randn(1, d))
-    elif mode == 'full':
-        va  = randn(d, d)
-        va  = dot(va, va.transpose())
-
-    input   = randn(n, d)
-    output  = gauss_den(input, mu, va)
-
-    import tables
-    h5file  = tables.openFile(file, "w")
-
-    h5file.createArray(h5file.root, 'input', input)
-    h5file.createArray(h5file.root, 'mu', mu)
-    h5file.createArray(h5file.root, 'va', va)
-    h5file.createArray(h5file.root, 'output', output)
-
-    h5file.close()
-
-def test_gauss_den():
-    """"""
-    # import tables
-    # import numpy as N
-    # 
-    # filename    = 'dendata.h5'
-
-    # # # Dimension 1
-    # # d   = 1
-    # # mu  = 1.0
-    # # va  = 2.0
-
-    # # X   = randn(1e3, 1)
-
-    # # Y   = gauss_den(X, mu, va)
-
-    # # h5file      = tables.openFile(filename, "w")
-
-    # # h5file.createArray(h5file.root, 'X', X)
-    # # h5file.createArray(h5file.root, 'mu', mu)
-    # # h5file.createArray(h5file.root, 'va', va)
-    # # h5file.createArray(h5file.root, 'Y', Y)
-
-    # # h5file.close()
-
-    # # # Dimension 2, diag
-    # # d   = 2
-    # # mu  = N.array([1.0, -2.0])
-    # # va  = N.array([1.0, 2.0])
-
-    # # X   = randn(1e3, 2)
-
-    # # Y   = gauss_den(X, mu, va)
-
-    # # h5file      = tables.openFile(filename, "w")
-
-    # # h5file.createArray(h5file.root, 'X', X)
-    # # h5file.createArray(h5file.root, 'mu', mu)
-    # # h5file.createArray(h5file.root, 'va', va)
-    # # h5file.createArray(h5file.root, 'Y', Y)
-
-    # # Dimension 2, full
-    # d   = 2
-    # mu  = N.array([[0.2, -1.0]])
-    # va  = N.array([[1.2, 0.1], [0.1, 0.5]])
-
-    # X   = randn(1e3, 2)
-
-    # Y   = gauss_den(X, mu, va)
-
-    # h5file      = tables.openFile(filename, "w")
-
-    # h5file.createArray(h5file.root, 'X', X)
-    # h5file.createArray(h5file.root, 'mu', mu)
-    # h5file.createArray(h5file.root, 'va', va)
-    # h5file.createArray(h5file.root, 'Y', Y)
-
-    # h5file.close()
-
-    import numpy.testing as testing
-    #=================
-    # Small test in 1d
-    #=================
-    va  = 2.0
-    mu  = 1.0
-    X   = N.linspace(-2, 2, 10)[:, N.NewAxis]
-
-    Yt  = N.array([0.02973257230591, 0.05512079811082, 0.09257745306945, 
-            0.14086453882683,
-            0.19418015562214, 0.24250166773127, 0.27436665745048, 0.28122547107069,
-            0.26114678964743, 0.21969564473386])
-
-    Y   = gauss_den(X, mu, va)
-    try:
-        testing.assert_array_almost_equal(Y, Yt)
-        print "1d test succeded"
-    except AssertionError:
-        print "test fails in 1d"
-
-    #============================
-    # Small test in 2d (diagonal)
-    #============================
-    mu  = N.atleast_2d([-1.0, 2.0])
-    va  = N.atleast_2d([2.0, 3.0])
-    X1  = N.linspace(-2, 2, 10)[:, N.NewAxis]
-    X2  = N.linspace(-1, 3, 10)[:, N.NewAxis]
-    X   = N.concatenate(([X1, X2]), 1)
-    
-    Yt  = N.array([0.01129091565384, 0.02025416837152, 0.03081845516786, 
-            0.03977576221540, 0.04354490552910, 0.04043592581117, 
-            0.03184994053539, 0.02127948225225, 0.01205937178755, 
-            0.00579694938623 ])
-
-    Y   = gauss_den(X, mu, va)
-    try:
-        testing.assert_array_almost_equal(Y, Yt)
-        print "2d diag test succeded"
-    except AssertionError:
-        print "test fails in 2d diag"
-
-    #============================
-    # Small test in 2d (full mat)
-    #============================
-    mu  = N.array([[0.2, -1.0]])
-    va  = N.array([[1.2, 0.1], [0.1, 0.5]])
-    X1  = N.linspace(-2, 2, 10)[:, N.NewAxis]
-    X2  = N.linspace(-3, 3, 10)[:, N.NewAxis]
-    X   = N.concatenate(([X1, X2]), 1)
-    
-    Yt  = N.array([0.00096157109751, 0.01368908714856,
-        0.07380823191162, 0.15072050533842, 
-        0.11656739937861, 0.03414436965525,
-        0.00378789836599, 0.00015915297541, 
-        0.00000253261067, 0.00000001526368])
-
-    Y   = gauss_den(X, mu, va)
-    try:
-        testing.assert_array_almost_equal(Y, Yt)
-        print "2d full test succeded"
-    except AssertionError:
-        print "test fails in 2d full"
-         
-
 if __name__ == "__main__":
     import pylab
 

Modified: trunk/Lib/sandbox/pyem/pyem/gauss_mix.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/gauss_mix.py	2006-10-12 13:47:46 UTC (rev 2279)
+++ trunk/Lib/sandbox/pyem/pyem/gauss_mix.py	2006-10-12 13:47:57 UTC (rev 2280)
@@ -1,5 +1,5 @@
 # /usr/bin/python
-# Last Change: Tue Aug 29 08:00 PM 2006 J
+# Last Change: Mon Sep 11 05:00 PM 2006 J
 
 # Module to implement GaussianMixture class.
 
@@ -28,6 +28,9 @@
 #   methods to be used, ie for implementing random index.
 #   - there is no check on internal state of the GM, that is does w, mu and va values
 #   make sense (eg singular values)
+#   - plot1d is still very rhough. There should be a sensible way to 
+#   modify the result plot (maybe returns a dic with global pdf, component pdf and
+#   fill matplotlib handles). Should be coherent with plot
 class GmParamError:
     """Exception raised for errors in gmm params
 
@@ -84,8 +87,8 @@
         if not d == self.d:
             raise GmParamError("Dimension of the given model is %d, expected %d" 
                     % (shape(d), shape(self.d)))
-        if not mode == self.mode:
-            raise GmParamError("Given covariance mode is %s, expected %d"
+        if not mode == self.mode and not d == 1:
+            raise GmParamError("Given covariance mode is %s, expected %s"
                     % (mode, self.mode))
         self.w  = weights
         self.mu = mu
@@ -260,8 +263,17 @@
         except ImportError:
             raise GmParamError("matplotlib not found, cannot plot...")
 
-    def plot1d(self, level = 0.5):
-        """TODO: this is not documented"""
+    def plot1d(self, level = 0.5, fill = 0, gpdf = 0):
+        """This function plots the pdfs of each component of the model. 
+        If gpdf is 1, also plots the global pdf. If fill is 1, fill confidence
+        areas using level argument as a level value
+        
+        Returns a dictionary h of plot handles so that their properties can
+        be modified (eg color, label, etc...):
+            - h['pdf'] is a list of lines, one line per component pdf
+            - h['gpdf'] is the line for the global pdf
+            - h['conf'] is a list of filling area
+        """
         # This is not optimized at all, may be slow. Should not be
         # difficult to make much faster, but it is late, and I am lazy
         if not self.d == 1:
@@ -271,36 +283,47 @@
         pval    = N.sqrt(self.va[:,0]) * nrm.ppf((1+level)/2)
 
         # Compute reasonable min/max for the normal pdf
-        mc  = 4
+        mc  = 3
         std = N.sqrt(self.va[:,0])
-        m   = N.amin(self.mu[:, 0] - 5 * std)
-        M   = N.amax(self.mu[:, 0] + 5 * std)
+        m   = N.amin(self.mu[:, 0] - mc * std)
+        M   = N.amax(self.mu[:, 0] + mc * std)
 
         np  = 500
         x   = N.linspace(m, M, np)
         Yf  = N.zeros(np)
         Yt  = N.zeros(np)
+
+        # Prepare the dic of plot handles to return
+        ks  = ['pdf', 'conf', 'gpdf']
+        hp  = dict((i,[]) for i in ks)
         try:
             import pylab as P
             for c in range(self.k):
                 y   = self.w[c]/(N.sqrt(2*N.pi) * std[c]) * \
                         N.exp(-(x-self.mu[c][0])**2/(2*std[c]**2))
                 Yt  += y
-                P.plot(x, y, 'r:')
-                #P.axvspan(-pval[c] + self.mu[c][0], pval[c] + self.mu[c][0], 
-                #        facecolor = 'b', alpha = 0.2)
-                id1 = -pval[c] + self.mu[c]
-                id2 = pval[c] + self.mu[c]
-                xc  = x[:, N.where(x>id1)[0]]
-                xc  = xc[:, N.where(xc<id2)[0]]
-                Yf  = self.w[c]/(N.sqrt(2*N.pi) * std[c]) * \
-                        N.exp(-(xc-self.mu[c][0])**2/(2*std[c]**2))
-                xc  = N.concatenate(([xc[0]], xc, [xc[-1]]))
-                Yf  = N.concatenate(([0], Yf, [0]))
-                P.fill(xc, Yf, facecolor = 'b', alpha = 0.2)
-                #P.fill([xc[0], xc[0], xc[-1], xc[-1]], 
-                #        [0, Yf[0], Yf[-1], 0], facecolor = 'b', alpha = 0.2)
-            P.plot(x, Yt, 'r')
+                h   = P.plot(x, y, 'r', label ='_nolegend_')
+                hp['pdf'].extend(h)
+                if fill:
+                    #P.axvspan(-pval[c] + self.mu[c][0], pval[c] + self.mu[c][0], 
+                    #        facecolor = 'b', alpha = 0.2)
+                    id1 = -pval[c] + self.mu[c]
+                    id2 = pval[c] + self.mu[c]
+                    xc  = x[:, N.where(x>id1)[0]]
+                    xc  = xc[:, N.where(xc<id2)[0]]
+                    Yf  = self.w[c]/(N.sqrt(2*N.pi) * std[c]) * \
+                            N.exp(-(xc-self.mu[c][0])**2/(2*std[c]**2))
+                    xc  = N.concatenate(([xc[0]], xc, [xc[-1]]))
+                    Yf  = N.concatenate(([0], Yf, [0]))
+                    h   = P.fill(xc, Yf, 
+                            facecolor = 'b', alpha = 0.1, label='_nolegend_')
+                    hp['conf'].extend(h)
+                    #P.fill([xc[0], xc[0], xc[-1], xc[-1]], 
+                    #        [0, Yf[0], Yf[-1], 0], facecolor = 'b', alpha = 0.2)
+            if gpdf:
+                h           = P.plot(x, Yt, 'r:', label='_nolegend_')
+                hp['gpdf']  = h
+            return hp
         except ImportError:
             raise GmParamError("matplotlib not found, cannot plot...")
 

Modified: trunk/Lib/sandbox/pyem/pyem/gmm_em.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/gmm_em.py	2006-10-12 13:47:46 UTC (rev 2279)
+++ trunk/Lib/sandbox/pyem/pyem/gmm_em.py	2006-10-12 13:47:57 UTC (rev 2280)
@@ -1,9 +1,9 @@
 # /usr/bin/python
-# Last Change: Tue Aug 29 03:00 PM 2006 J
+# Last Change: Mon Sep 11 05:00 PM 2006 J
 
 # TODO:
-#   - which methods to avoid va shrinking to 0 ? There are several options, not sure which
-#   ones are appropriates
+#   - which methods to avoid va shrinking to 0 ? There are several options, 
+#   not sure which ones are appropriates
 #   - improve EM trainer
 #   - online EM
 
@@ -82,9 +82,10 @@
         else:
             raise GmmParamError("mode " + str(mode) + " not recognized")
 
-        self.gm.w   = w
-        self.gm.mu  = mu
-        self.gm.va  = va
+        self.gm.set_param(w, mu, va)
+        #self.gm.w   = w
+        #self.gm.mu  = mu
+        #self.gm.va  = va
 
     def init_random(self, data):
         """ Init the model at random."""
@@ -98,9 +99,7 @@
             raise GmmParamError("""init_random not implemented for
                     mode %s yet""", mode)
 
-        self.gm.w   = w
-        self.gm.mu  = mu
-        self.gm.va  = va
+        self.gm.set_param(w, mu, va)
 
     # TODO: 
     #   - format of parameters ? For variances, list of variances matrix,
@@ -265,9 +264,9 @@
     #   - mode: Mode of covariance matrix: full or diag
     #   - nframes: number of frames (frame = one data point = one
     #   row of d elements
-    k       = 3 
-    d       = 2         
-    mode    = 'full'        
+    k       = 2 
+    d       = 1
+    mode    = 'full'
     nframes = 1e3
 
     #+++++++++++++++++++++++++++++++++++++++++++
@@ -316,33 +315,47 @@
     print "drawing..."
     import pylab as P
     P.subplot(2, 1, 1)
-    # Draw what is happening
-    P.plot(data[:, 0], data[:, 1], '.', label = '_nolegend_')
 
-    # Real confidence ellipses
-    Xre, Yre  = gm.conf_ellipses()
-    P.plot(Xre[0], Yre[0], 'g', label = 'true confidence ellipsoides')
-    for i in range(1,k):
-        P.plot(Xre[i], Yre[i], 'g', label = '_nolegend_')
+    if not d == 1:
+        # Draw what is happening
+        P.plot(data[:, 0], data[:, 1], '.', label = '_nolegend_')
 
-    # Initial confidence ellipses as found by kmean
-    X0e, Y0e  = gm0.conf_ellipses()
-    P.plot(X0e[0], Y0e[0], 'k', label = 'initial confidence ellipsoides')
-    for i in range(1,k):
-        P.plot(X0e[i], Y0e[i], 'k', label = '_nolegend_')
+        # Real confidence ellipses
+        Xre, Yre  = gm.conf_ellipses()
+        P.plot(Xre[0], Yre[0], 'g', label = 'true confidence ellipsoides')
+        for i in range(1,k):
+            P.plot(Xre[i], Yre[i], 'g', label = '_nolegend_')
 
-    # Values found by EM
-    Xe, Ye  = lgm.conf_ellipses()
-    P.plot(Xe[0], Ye[0], 'r', label = 'confidence ellipsoides found by EM')
-    for i in range(1,k):
-        P.plot(Xe[i], Ye[i], 'r', label = '_nolegend_')
-    P.legend(loc = 0)
+        # Initial confidence ellipses as found by kmean
+        X0e, Y0e  = gm0.conf_ellipses()
+        P.plot(X0e[0], Y0e[0], 'k', label = 'initial confidence ellipsoides')
+        for i in range(1,k):
+            P.plot(X0e[i], Y0e[i], 'k', label = '_nolegend_')
 
-    #from scipy.cluster.vq import kmeans
-    #code    = kmeans(data, k)[0]
-    #print code
-    #P.plot(code[:,0], code[:, 1], 'oy')
-    #P.plot(gm0.mu[:,0], gm0.mu[:, 1], 'ok')
+        # Values found by EM
+        Xe, Ye  = lgm.conf_ellipses()
+        P.plot(Xe[0], Ye[0], 'r', label = 'confidence ellipsoides found by EM')
+        for i in range(1,k):
+            P.plot(Xe[i], Ye[i], 'r', label = '_nolegend_')
+        P.legend(loc = 0)
+    else:
+        # Real confidence ellipses
+        h   = gm.plot1d()
+        [i.set_color('g') for i in h['pdf']]
+        h['pdf'][0].set_label('true pdf')
+
+        # Initial confidence ellipses as found by kmean
+        h0  = gm0.plot1d()
+        [i.set_color('k') for i in h0['pdf']]
+        h0['pdf'][0].set_label('initial pdf')
+
+        # Values found by EM
+        hl  = lgm.plot1d(fill = 1, level = 0.66)
+        [i.set_color('r') for i in hl['pdf']]
+        hl['pdf'][0].set_label('pdf found by EM')
+
+        P.legend(loc = 0)
+
     P.subplot(2, 1, 2)
     P.plot(like)
     P.title('log likelihood')

Modified: trunk/Lib/sandbox/pyem/pyem/info.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/info.py	2006-10-12 13:47:46 UTC (rev 2279)
+++ trunk/Lib/sandbox/pyem/pyem/info.py	2006-10-12 13:47:57 UTC (rev 2280)
@@ -1,9 +1,35 @@
 """
-Routines for for Gaussian Mixture Models
-and Expectation Maximization learning 
-===================================
+Routines for Gaussian Mixture Models
+and learning with Expectation Maximization 
+==========================================
 
+This module contains classes and function to compute multivariate Gaussian densities
+(diagonal and full covariance matrices), Gaussian mixtures, Gaussian mixtures models
+and an Em trainer.
+
+More specifically, the module contains the following file:
+
+- densities.py: functions to compute multivariate Gaussian pdf and ellipsoides of
+confidence (gauss_den)
+- gauss_mix.py: defines the GM (Gaussian Mixture) class. A Gaussian Mixture can be
+created from its parameters weights, mean and variances, or from its meta parameters
+d (dimension of the Gaussian) and k (number of components in the mixture). A Gaussian
+Model can then be sampled or plot (if d>1, plot confidence ellipsoids projected on 
+2 chosen dimensions, if d == 1, plot the pdf of each component and fill the zone
+of confidence for a given level)
+- gmm_em.py: defines a class GMM (Gaussian Mixture Model). This class is constructed
+from a GM model gm, and can be used to train gm. The GMM can be initiated by
+kmean or at random, and can compute sufficient statistics, and update its parameters
+from the sufficient statistics.
+- kmean.py: implements a kmean algorithm. We cannot use scipy.cluster.vq kmeans, since
+its does not give membership of observations.
+
+Example of use: simply execute gmm_em.py for an example of training a GM and plotting
+the results compared to true values
+
 Copyright: David Cournapeau 2006
 License: BSD-style (see LICENSE.txt in main source directory)
 """
 
+depends = ['linalg', 'stats']
+ignore  = False

Added: trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py
===================================================================
--- trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py	2006-10-12 13:47:46 UTC (rev 2279)
+++ trunk/Lib/sandbox/pyem/pyem/tests/test_densities.py	2006-10-12 13:47:57 UTC (rev 2280)
@@ -0,0 +1,145 @@
+def generate_test_data(n, d, mode = 'diag', file='test.dat'):
+    """Generate a set of data of dimension d, with n frames,
+    that is input data, mean, var and output of gden, so that
+    other implementations can be tested against"""
+    mu  = randn(1, d)
+    if mode == 'diag':
+        va  = abs(randn(1, d))
+    elif mode == 'full':
+        va  = randn(d, d)
+        va  = dot(va, va.transpose())
+
+    input   = randn(n, d)
+    output  = gauss_den(input, mu, va)
+
+    import tables
+    h5file  = tables.openFile(file, "w")
+
+    h5file.createArray(h5file.root, 'input', input)
+    h5file.createArray(h5file.root, 'mu', mu)
+    h5file.createArray(h5file.root, 'va', va)
+    h5file.createArray(h5file.root, 'output', output)
+
+    h5file.close()
+
+def test_gauss_den():
+    """"""
+    # import tables
+    # import numpy as N
+    # 
+    # filename    = 'dendata.h5'
+
+    # # # Dimension 1
+    # # d   = 1
+    # # mu  = 1.0
+    # # va  = 2.0
+
+    # # X   = randn(1e3, 1)
+
+    # # Y   = gauss_den(X, mu, va)
+
+    # # h5file      = tables.openFile(filename, "w")
+
+    # # h5file.createArray(h5file.root, 'X', X)
+    # # h5file.createArray(h5file.root, 'mu', mu)
+    # # h5file.createArray(h5file.root, 'va', va)
+    # # h5file.createArray(h5file.root, 'Y', Y)
+
+    # # h5file.close()
+
+    # # # Dimension 2, diag
+    # # d   = 2
+    # # mu  = N.array([1.0, -2.0])
+    # # va  = N.array([1.0, 2.0])
+
+    # # X   = randn(1e3, 2)
+
+    # # Y   = gauss_den(X, mu, va)
+
+    # # h5file      = tables.openFile(filename, "w")
+
+    # # h5file.createArray(h5file.root, 'X', X)
+    # # h5file.createArray(h5file.root, 'mu', mu)
+    # # h5file.createArray(h5file.root, 'va', va)
+    # # h5file.createArray(h5file.root, 'Y', Y)
+
+    # # Dimension 2, full
+    # d   = 2
+    # mu  = N.array([[0.2, -1.0]])
+    # va  = N.array([[1.2, 0.1], [0.1, 0.5]])
+
+    # X   = randn(1e3, 2)
+
+    # Y   = gauss_den(X, mu, va)
+
+    # h5file      = tables.openFile(filename, "w")
+
+    # h5file.createArray(h5file.root, 'X', X)
+    # h5file.createArray(h5file.root, 'mu', mu)
+    # h5file.createArray(h5file.root, 'va', va)
+    # h5file.createArray(h5file.root, 'Y', Y)
+
+    # h5file.close()
+
+    import numpy.testing as testing
+    #=================
+    # Small test in 1d
+    #=================
+    va  = 2.0
+    mu  = 1.0
+    X   = N.linspace(-2, 2, 10)[:, N.NewAxis]
+
+    Yt  = N.array([0.02973257230591, 0.05512079811082, 0.09257745306945, 
+            0.14086453882683,
+            0.19418015562214, 0.24250166773127, 0.27436665745048, 0.28122547107069,
+            0.26114678964743, 0.21969564473386])
+
+    Y   = gauss_den(X, mu, va)
+    try:
+        testing.assert_array_almost_equal(Y, Yt)
+        print "1d test succeded"
+    except AssertionError:
+        print "test fails in 1d"
+
+    #============================
+    # Small test in 2d (diagonal)
+    #============================
+    mu  = N.atleast_2d([-1.0, 2.0])
+    va  = N.atleast_2d([2.0, 3.0])
+    X1  = N.linspace(-2, 2, 10)[:, N.NewAxis]
+    X2  = N.linspace(-1, 3, 10)[:, N.NewAxis]
+    X   = N.concatenate(([X1, X2]), 1)
+    
+    Yt  = N.array([0.01129091565384, 0.02025416837152, 0.03081845516786, 
+            0.03977576221540, 0.04354490552910, 0.04043592581117, 
+            0.03184994053539, 0.02127948225225, 0.01205937178755, 
+            0.00579694938623 ])
+
+    Y   = gauss_den(X, mu, va)
+    try:
+        testing.assert_array_almost_equal(Y, Yt)
+        print "2d diag test succeded"
+    except AssertionError:
+        print "test fails in 2d diag"
+
+    #============================
+    # Small test in 2d (full mat)
+    #============================
+    mu  = N.array([[0.2, -1.0]])
+    va  = N.array([[1.2, 0.1], [0.1, 0.5]])
+    X1  = N.linspace(-2, 2, 10)[:, N.NewAxis]
+    X2  = N.linspace(-3, 3, 10)[:, N.NewAxis]
+    X   = N.concatenate(([X1, X2]), 1)
+    
+    Yt  = N.array([0.00096157109751, 0.01368908714856,
+        0.07380823191162, 0.15072050533842, 
+        0.11656739937861, 0.03414436965525,
+        0.00378789836599, 0.00015915297541, 
+        0.00000253261067, 0.00000001526368])
+
+    Y   = gauss_den(X, mu, va)
+    try:
+        testing.assert_array_almost_equal(Y, Yt)
+        print "2d full test succeded"
+    except AssertionError:
+        print "test fails in 2d full"



More information about the Scipy-svn mailing list