[Scipy-svn] r2280 - in trunk/Lib/sandbox/pyem: . pyem pyem/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Oct 12 09:48:05 EDT 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