← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/pyem/densities.py

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pyem/densities.py
1
 
#! /usr/bin/python
2
 
#
3
 
# Copyrighted David Cournapeau
4
 
# Last Change: Fri Nov 10 10:00 AM 2006 J
5
 
 
6
 
import numpy as N
7
 
import numpy.linalg as lin
8
 
from numpy.random import randn
9
 
from scipy.stats import chi2
10
 
 
11
 
# Error classes
12
 
class DenError(Exception):
13
 
    """Base class for exceptions in this module.
14
 
    
15
 
    Attributes:
16
 
        expression -- input expression in which the error occurred
17
 
        message -- explanation of the error"""
18
 
    def __init__(self, message):
19
 
        self.message    = message
20
 
    
21
 
    def __str__(self):
22
 
        return self.message
23
 
 
24
 
# The following function do all the fancy stuff to check that parameters
25
 
# are Ok, and call the right implementation if args are OK.
26
 
def gauss_den(x, mu, va, log = False):
27
 
    """ Compute multivariate Gaussian density at points x for 
28
 
    mean mu and variance va.
29
 
    
30
 
    Vector are row vectors, except va which can be a matrix
31
 
    (row vector variance for diagonal variance)
32
 
    
33
 
    If log is True, than the log density is returned 
34
 
    (useful for underflow ?)"""
35
 
    mu  = N.atleast_2d(mu)
36
 
    va  = N.atleast_2d(va)
37
 
    x   = N.atleast_2d(x)
38
 
    
39
 
    #=======================#
40
 
    # Checking parameters   #
41
 
    #=======================#
42
 
    if len(N.shape(mu)) != 2:
43
 
        raise DenError("mu is not rank 2")
44
 
        
45
 
    if len(N.shape(va)) != 2:
46
 
        raise DenError("va is not rank 2")
47
 
        
48
 
    if len(N.shape(x)) != 2:
49
 
        raise DenError("x is not rank 2")
50
 
        
51
 
    (n, d)      = x.shape
52
 
    (dm0, dm1)  = mu.shape
53
 
    (dv0, dv1)  = va.shape
54
 
    
55
 
    # Check x and mu same dimension
56
 
    if dm0 != 1:
57
 
        msg = "mean must be a row vector!"
58
 
        raise DenError(msg)
59
 
    if dm1 != d:
60
 
        msg = "x and mu not same dim"
61
 
        raise DenError(msg)
62
 
    # Check va and mu same size
63
 
    if dv1 != d:
64
 
        msg = "mu and va not same dim"
65
 
        raise DenError(msg)
66
 
    if dv0 != 1 and dv0 != d:
67
 
        msg = "va not square"
68
 
        raise DenError(msg)
69
 
 
70
 
    #===============#
71
 
    # Computation   #
72
 
    #===============#
73
 
    if d == 1:
74
 
        # scalar case
75
 
        return _scalar_gauss_den(x[:, 0], mu[0, 0], va[0, 0], log)
76
 
    elif dv0 == 1:
77
 
        # Diagonal matrix case
78
 
        return _diag_gauss_den(x, mu, va, log)
79
 
    elif dv1 == dv0:
80
 
        # full case
81
 
        return  _full_gauss_den(x, mu, va, log)
82
 
    else:
83
 
        raise DenError("variance mode not recognized, this is a bug")
84
 
 
85
 
# Those 3 functions do almost all the actual computation
86
 
def _scalar_gauss_den(x, mu, va, log):
87
 
    """ This function is the actual implementation
88
 
    of gaussian pdf in scalar case. It assumes all args
89
 
    are conformant, so it should not be used directly
90
 
    
91
 
    Call gauss_den instead"""
92
 
    d       = mu.size
93
 
    inva    = 1/va
94
 
    fac     = (2*N.pi) ** (-d/2.0) * N.sqrt(inva)
95
 
    y       = ((x-mu) ** 2) * -0.5 * inva
96
 
    if not log:
97
 
        y   = fac * N.exp(y)
98
 
    else:
99
 
        y   = y + log(fac)
100
 
 
101
 
    return y
102
 
    
103
 
#from ctypes import cdll, c_uint, c_int, c_double, POINTER
104
 
#_gden   = cdll.LoadLibrary('src/libgden.so')
105
 
#_gden.gden_diag.restype     = c_int
106
 
#_gden.gden_diag.argtypes    = [POINTER(c_double), c_uint, c_uint,
107
 
#        POINTER(c_double), POINTER(c_double), POINTER(c_double)]
108
 
 
109
 
def _diag_gauss_den(x, mu, va, log):
110
 
    """ This function is the actual implementation
111
 
    of gaussian pdf in scalar case. It assumes all args
112
 
    are conformant, so it should not be used directly
113
 
    
114
 
    Call gauss_den instead"""
115
 
    # Diagonal matrix case
116
 
    d   = mu.size
117
 
    n   = x.shape[0]
118
 
    if not log:
119
 
        inva    = 1/va[0,0]
120
 
        fac     = (2*N.pi) ** (-d/2.0) * N.sqrt(inva)
121
 
        y       =  (x[:,0] - mu[0,0]) ** 2 * inva * -0.5
122
 
        for i in range(1, d):
123
 
            inva    = 1/va[0,i]
124
 
            fac     *= N.sqrt(inva)
125
 
            y       += (x[:,i] - mu[0,i]) ** 2 * inva * -0.5
126
 
        y   = fac * N.exp(y)
127
 
    else:
128
 
        y   = _scalar_gauss_den(x[:,0], mu[0,0], va[0,0], log)
129
 
        for i in range(1, d):
130
 
            y    +=  _scalar_gauss_den(x[:,i], mu[0,i], va[0,i], log)
131
 
    return y
132
 
 
133
 
def _full_gauss_den(x, mu, va, log):
134
 
    """ This function is the actual implementation
135
 
    of gaussian pdf in full matrix case. 
136
 
    
137
 
    It assumes all args are conformant, so it should 
138
 
    not be used directly Call gauss_den instead
139
 
    
140
 
    Does not check if va is definite positive (on inversible 
141
 
    for that matter), so the inverse computation and/or determinant
142
 
    would throw an exception."""
143
 
    d       = mu.size
144
 
    inva    = lin.inv(va)
145
 
    fac     = 1 / N.sqrt( (2*N.pi) ** d * N.fabs(lin.det(va)))
146
 
 
147
 
    # # Slow version
148
 
    # n       = N.size(x, 0)
149
 
    # y       = N.zeros(n)
150
 
    # for i in range(n):
151
 
    #     y[i] = N.dot(x[i,:],
152
 
    #              N.dot(inva, N.transpose(x[i,:])))
153
 
    # y *= -0.5
154
 
 
155
 
    # we are using a trick with sum to "emulate" 
156
 
    # the matrix multiplication inva * x without any explicit loop
157
 
    y   = N.dot((x-mu), inva)
158
 
    y   = -0.5 * N.sum(y * (x-mu), 1)
159
 
 
160
 
    if not log:
161
 
        y   = fac * N.exp(y)
162
 
    else:
163
 
        y   = y + N.log(fac)
164
 
 
165
 
    return y
166
 
 
167
 
# To plot a confidence ellipse from multi-variate gaussian pdf
168
 
def gauss_ell(mu, va, dim = [0, 1], npoints = 100, level = 0.39):
169
 
    """ Given a mean and covariance for multi-variate
170
 
    gaussian, returns npoints points for the ellipse
171
 
    of confidence given by level (all points will be inside
172
 
    the ellipsoides with a probability equal to level)
173
 
    
174
 
    Returns the coordinate x and y of the ellipse"""
175
 
    
176
 
    mu      = N.atleast_1d(mu)
177
 
    va      = N.atleast_1d(va)
178
 
    c       = N.array(dim)
179
 
 
180
 
    if mu.size == va.size:
181
 
        mode    = 'diag'
182
 
    else:
183
 
        if va.ndim == 2:
184
 
            if va.shape[0] == va.shape[1]:
185
 
                mode    = 'full'
186
 
            else:
187
 
                raise DenError("variance not square")
188
 
        else:
189
 
            raise DenError("mean and variance are not dim conformant")
190
 
 
191
 
    chi22d  = chi2(2)
192
 
    mahal   = N.sqrt(chi22d.ppf(level))
193
 
    
194
 
    # Generates a circle of npoints
195
 
    theta   = N.linspace(0, 2 * N.pi, npoints)
196
 
    circle  = mahal * N.array([N.cos(theta), N.sin(theta)])
197
 
 
198
 
    # Get the dimension which we are interested in:
199
 
    mu  = mu[dim]
200
 
    if mode == 'diag':
201
 
        va      = va[dim]
202
 
        elps    = N.outer(mu, N.ones(npoints))
203
 
        elps    += N.dot(N.diag(N.sqrt(va)), circle)
204
 
    elif mode == 'full':
205
 
        va  = va[c,:][:,c]
206
 
        # Method: compute the cholesky decomp of each cov matrix, that is
207
 
        # compute cova such as va = cova * cova' 
208
 
        # WARN: scipy is different than matlab here, as scipy computes a lower
209
 
        # triangular cholesky decomp: 
210
 
        #   - va = cova * cova' (scipy)
211
 
        #   - va = cova' * cova (matlab)
212
 
        # So take care when comparing results with matlab !
213
 
        cova    = lin.cholesky(va)
214
 
        elps    = N.outer(mu, N.ones(npoints))
215
 
        elps    += N.dot(cova, circle)
216
 
    else:
217
 
        raise DenParam("var mode not recognized")
218
 
 
219
 
    return elps[0, :], elps[1, :]
220
 
 
221
 
if __name__ == "__main__":
222
 
    import pylab
223
 
 
224
 
    #=========================================
225
 
    # Test plotting a simple diag 2d variance:
226
 
    #=========================================
227
 
    va  = N.array([5, 3])
228
 
    mu  = N.array([2, 3])
229
 
 
230
 
    # Generate a multivariate gaussian of mean mu and covariance va
231
 
    X       = randn(1e3, 2)
232
 
    Yc      = N.dot(N.diag(N.sqrt(va)), X.transpose())
233
 
    Yc      = Yc.transpose() + mu
234
 
 
235
 
    # Plotting
236
 
    Xe, Ye  = gauss_ell(mu, va, npoints = 100)
237
 
    pylab.figure()
238
 
    pylab.plot(Yc[:, 0], Yc[:, 1], '.')
239
 
    pylab.plot(Xe, Ye, 'r')
240
 
 
241
 
    #=========================================
242
 
    # Test plotting a simple full 2d variance:
243
 
    #=========================================
244
 
    va  = N.array([[0.2, 0.1],[0.1, 0.5]])
245
 
    mu  = N.array([0, 3])
246
 
 
247
 
    # Generate a multivariate gaussian of mean mu and covariance va
248
 
    X       = randn(1e3, 2)
249
 
    Yc      = N.dot(lin.cholesky(va), X.transpose())
250
 
    Yc      = Yc.transpose() + mu
251
 
 
252
 
    # Plotting
253
 
    Xe, Ye  = gauss_ell(mu, va, npoints = 100, level=0.95)
254
 
    pylab.figure()
255
 
    pylab.plot(Yc[:, 0], Yc[:, 1], '.')
256
 
    pylab.plot(Xe, Ye, 'r')
257
 
    pylab.show()