3
# Copyrighted David Cournapeau
4
# Last Change: Sat Jun 02 07:00 PM 2007 J
6
# New version, with default numpy ordering.
9
import numpy.linalg as lin
10
from numpy.random import randn
11
from scipy.stats import chi2
14
class DenError(Exception):
15
"""Base class for exceptions in this module.
18
expression -- input expression in which the error occurred
19
message -- explanation of the error"""
20
def __init__(self, message):
21
self.message = message
29
# The following function do all the fancy stuff to check that parameters
30
# are Ok, and call the right implementation if args are OK.
31
def gauss_den(x, mu, va, log = False, axis = -1):
32
""" Compute multivariate Gaussian density at points x for
33
mean mu and variance va along specified axis:
36
* mean must be rank 0 (1d) or rank 1 (multi variate gaussian)
37
* va must be rank 0 (1d), rank 1(multi variate, diag covariance) or rank 2
38
(multivariate, full covariance).
39
* in 1 dimension case, any rank for mean and va is ok, as long as their size
40
is 1 (eg they contain only 1 element)
42
Caution: if x is rank 1, it is assumed you have a 1d problem. You cannot compute
43
the gaussian densities of only one sample of dimension d; for this, you have
46
If log is True, than the log density is returned
47
(useful for underflow ?)"""
49
# If data is rank 1, then we have 1 dimension problem.
53
if not N.size(mu) == 1:
54
raise DenError("for 1 dimension problem, mean must have only one element")
56
if not N.size(va) == 1:
57
raise DenError("for 1 dimension problem, mean must have only one element")
59
return _scalar_gauss_den(x, mu, va, log)
60
# If data is rank 2, then we may have 1 dimension or multi-variate problem
62
oaxis = (axis + 1) % 2
66
# Get away with 1d case now
68
return _scalar_gauss_den(x, mu, va, log)
70
# Now, d > 1 (numpy attributes should be valid on mean and va now)
71
if not N.size(mu) == d or not mu.ndim == 1:
72
raise DenError("data is %d dimension, but mean's shape is %s" \
73
% (d, N.shape(mu)) + " (should be (%d,))" % d)
75
isfull = (va.ndim == 2)
76
if not (N.size(va) == d or (isfull and va.shape[0] == va.shape[1] == d)):
77
raise DenError("va has an invalid shape or number of elements")
82
return _full_gauss_den(x, mu[:, N.newaxis], va, log, axis)
84
return _full_gauss_den(x, mu, va, log, axis)
86
return _diag_gauss_den(x, mu, va, log, axis)
88
raise RuntimeError("Sorry, only rank up to 2 supported")
90
# To plot a confidence ellipse from multi-variate gaussian pdf
91
def gauss_ell(mu, va, dim = [0, 1], npoints = 100, level = 0.39):
92
""" Given a mean and covariance for multi-variate
93
gaussian, returns npoints points for the ellipse
94
of confidence given by level (all points will be inside
95
the ellipsoides with a probability equal to level)
97
Returns the coordinate x and y of the ellipse"""
102
raise RuntimeError("this function only make sense for dimension 2 and more")
104
if mu.size == va.size:
108
if va.shape[0] == va.shape[1]:
111
raise DenError("variance not square")
113
raise DenError("mean and variance are not dim conformant")
115
# If X ~ N(mu, va), then [X` * va^(-1/2) * X] ~ Chi2
117
mahal = N.sqrt(chi22d.ppf(level))
119
# Generates a circle of npoints
120
theta = N.linspace(0, 2 * N.pi, npoints)
121
circle = mahal * N.array([N.cos(theta), N.sin(theta)])
123
# Get the dimension which we are interested in:
127
elps = N.outer(mu, N.ones(npoints))
128
elps += N.dot(N.diag(N.sqrt(va)), circle)
131
# Method: compute the cholesky decomp of each cov matrix, that is
132
# compute cova such as va = cova * cova'
133
# WARN: scipy is different than matlab here, as scipy computes a lower
134
# triangular cholesky decomp:
135
# - va = cova * cova' (scipy)
136
# - va = cova' * cova (matlab)
137
# So take care when comparing results with matlab !
138
cova = lin.cholesky(va)
139
elps = N.outer(mu, N.ones(npoints))
140
elps += N.dot(cova, circle)
142
raise DenParam("var mode not recognized")
144
return elps[0, :], elps[1, :]
149
# Those 3 functions do almost all the actual computation
150
def _scalar_gauss_den(x, mu, va, log):
151
""" This function is the actual implementation
152
of gaussian pdf in scalar case. It assumes all args
153
are conformant, so it should not be used directly
155
Call gauss_den instead"""
157
fac = (2*N.pi) ** (-1/2.0) * N.sqrt(inva)
158
y = ((x-mu) ** 2) * -0.5 * inva
160
y = fac * N.exp(y.ravel())
166
def _diag_gauss_den(x, mu, va, log, axis):
167
""" This function is the actual implementation
168
of gaussian pdf in scalar case. It assumes all args
169
are conformant, so it should not be used directly
171
Call gauss_den instead"""
172
# Diagonal matrix case
175
x = N.swapaxes(x, 0, 1)
179
fac = (2*N.pi) ** (-d/2.0) * N.sqrt(inva)
180
y = (x[0] - mu[0]) ** 2 * inva * -0.5
181
for i in range(1, d):
184
y += (x[i] - mu[i]) ** 2 * inva * -0.5
187
y = _scalar_gauss_den(x[0], mu[0], va[0], log)
188
for i in range(1, d):
189
y += _scalar_gauss_den(x[i], mu[i], va[i], log)
193
def _full_gauss_den(x, mu, va, log, axis):
194
""" This function is the actual implementation
195
of gaussian pdf in full matrix case.
197
It assumes all args are conformant, so it should
198
not be used directly Call gauss_den instead
200
Does not check if va is definite positive (on inversible
201
for that matter), so the inverse computation and/or determinant
202
would throw an exception."""
205
fac = 1 / N.sqrt( (2*N.pi) ** d * N.fabs(lin.det(va)))
207
# # Slow version (does not work since version 0.6)
211
# y[i] = N.dot(x[i,:],
212
# N.dot(inva, N.transpose(x[i,:])))
215
# we are using a trick with sum to "emulate"
216
# the matrix multiplication inva * x without any explicit loop
218
y = N.dot(inva, (x-mu))
219
y = -0.5 * N.sum(y * (x-mu), 0)
221
y = N.dot((x-mu), inva)
222
y = -0.5 * N.sum(y * (x-mu), 1)
231
if __name__ == "__main__":
234
#=========================================
235
# Test plotting a simple diag 2d variance:
236
#=========================================
240
# Generate a multivariate gaussian of mean mu and covariance va
242
Yc = N.dot(N.diag(N.sqrt(va)), X)
243
Yc = Yc.transpose() + mu
246
Xe, Ye = gauss_ell(mu, va, npoints = 100)
248
pylab.plot(Yc[:, 0], Yc[:, 1], '.')
249
pylab.plot(Xe, Ye, 'r')
251
#=========================================
252
# Test plotting a simple full 2d variance:
253
#=========================================
254
va = N.array([[0.2, 0.1],[0.1, 0.5]])
257
# Generate a multivariate gaussian of mean mu and covariance va
259
Yc = N.dot(lin.cholesky(va), X.transpose())
260
Yc = Yc.transpose() + mu
263
Xe, Ye = gauss_ell(mu, va, npoints = 100, level=0.95)
265
pylab.plot(Yc[:, 0], Yc[:, 1], '.')
266
pylab.plot(Xe, Ye, 'r')