← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to scipy/sandbox/pyem/densities2.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:
scipy/sandbox/pyem/densities2.py
 
1
#! /usr/bin/python
 
2
#
 
3
# Copyrighted David Cournapeau
 
4
# Last Change: Sat Jun 02 07:00 PM 2007 J
 
5
 
 
6
# New version, with default numpy ordering.
 
7
 
 
8
import numpy as N
 
9
import numpy.linalg as lin
 
10
from numpy.random import randn
 
11
from scipy.stats import chi2
 
12
 
 
13
# Error classes
 
14
class DenError(Exception):
 
15
    """Base class for exceptions in this module.
 
16
    
 
17
    Attributes:
 
18
        expression -- input expression in which the error occurred
 
19
        message -- explanation of the error"""
 
20
    def __init__(self, message):
 
21
        self.message    = message
 
22
    
 
23
    def __str__(self):
 
24
        return self.message
 
25
 
 
26
#============
 
27
# Public API
 
28
#============
 
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:
 
34
        
 
35
    requirements:
 
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)
 
41
    
 
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
 
44
    to use a rank 2 !
 
45
 
 
46
    If log is True, than the log density is returned 
 
47
    (useful for underflow ?)"""
 
48
 
 
49
    # If data is rank 1, then we have 1 dimension problem.
 
50
    if x.ndim == 1:
 
51
        d   = 1
 
52
        n   = x.size
 
53
        if not N.size(mu) == 1:
 
54
            raise DenError("for 1 dimension problem, mean must have only one element")
 
55
            
 
56
        if not N.size(va) == 1:
 
57
            raise DenError("for 1 dimension problem, mean must have only one element")
 
58
        
 
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
 
61
    elif x.ndim == 2:
 
62
        oaxis   = (axis + 1) % 2
 
63
        n       = x.shape[axis]
 
64
        d       = x.shape[oaxis]
 
65
 
 
66
        # Get away with 1d case now
 
67
        if d == 1:
 
68
            return _scalar_gauss_den(x, mu, va, log)
 
69
 
 
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)
 
74
            
 
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")
 
78
 
 
79
        if isfull:
 
80
            # Compute along rows
 
81
            if oaxis == 0:
 
82
                return  _full_gauss_den(x, mu[:, N.newaxis], va, log, axis)
 
83
            else:
 
84
                return  _full_gauss_den(x, mu, va, log, axis)
 
85
        else:
 
86
            return _diag_gauss_den(x, mu, va, log, axis)
 
87
    else:
 
88
        raise RuntimeError("Sorry, only rank up to 2 supported")
 
89
        
 
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)
 
96
    
 
97
    Returns the coordinate x and y of the ellipse"""
 
98
    
 
99
    c       = N.array(dim)
 
100
 
 
101
    if mu.size < 2:
 
102
        raise RuntimeError("this function only make sense for dimension 2 and more")
 
103
 
 
104
    if mu.size == va.size:
 
105
        mode    = 'diag'
 
106
    else:
 
107
        if va.ndim == 2:
 
108
            if va.shape[0] == va.shape[1]:
 
109
                mode    = 'full'
 
110
            else:
 
111
                raise DenError("variance not square")
 
112
        else:
 
113
            raise DenError("mean and variance are not dim conformant")
 
114
 
 
115
    # If X ~ N(mu, va), then [X` * va^(-1/2) * X] ~ Chi2
 
116
    chi22d  = chi2(2)
 
117
    mahal   = N.sqrt(chi22d.ppf(level))
 
118
    
 
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)])
 
122
 
 
123
    # Get the dimension which we are interested in:
 
124
    mu  = mu[dim]
 
125
    if mode == 'diag':
 
126
        va      = va[dim]
 
127
        elps    = N.outer(mu, N.ones(npoints))
 
128
        elps    += N.dot(N.diag(N.sqrt(va)), circle)
 
129
    elif mode == 'full':
 
130
        va  = va[c,:][:,c]
 
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)
 
141
    else:
 
142
        raise DenParam("var mode not recognized")
 
143
 
 
144
    return elps[0, :], elps[1, :]
 
145
 
 
146
#=============
 
147
# Private Api
 
148
#=============
 
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
 
154
    
 
155
    Call gauss_den instead"""
 
156
    inva    = 1/va
 
157
    fac     = (2*N.pi) ** (-1/2.0) * N.sqrt(inva)
 
158
    y       = ((x-mu) ** 2) * -0.5 * inva
 
159
    if not log:
 
160
        y   = fac * N.exp(y.ravel())
 
161
    else:
 
162
        y   = y + log(fac)
 
163
 
 
164
    return y
 
165
    
 
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
 
170
    
 
171
    Call gauss_den instead"""
 
172
    # Diagonal matrix case
 
173
    d   = mu.size
 
174
    if axis % 2 == 0:
 
175
        x  = N.swapaxes(x, 0, 1)
 
176
 
 
177
    if not log:
 
178
        inva    = 1/va[0]
 
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):
 
182
            inva    = 1/va[i]
 
183
            fac     *= N.sqrt(inva)
 
184
            y       += (x[i] - mu[i]) ** 2 * inva * -0.5
 
185
        y   = fac * N.exp(y)
 
186
    else:
 
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)
 
190
 
 
191
    return y
 
192
 
 
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. 
 
196
    
 
197
    It assumes all args are conformant, so it should 
 
198
    not be used directly Call gauss_den instead
 
199
    
 
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."""
 
203
    d       = mu.size
 
204
    inva    = lin.inv(va)
 
205
    fac     = 1 / N.sqrt( (2*N.pi) ** d * N.fabs(lin.det(va)))
 
206
 
 
207
    # # Slow version (does not work since version 0.6)
 
208
    # n       = N.size(x, 0)
 
209
    # y       = N.zeros(n)
 
210
    # for i in range(n):
 
211
    #     y[i] = N.dot(x[i,:],
 
212
    #              N.dot(inva, N.transpose(x[i,:])))
 
213
    # y *= -0.5
 
214
 
 
215
    # we are using a trick with sum to "emulate" 
 
216
    # the matrix multiplication inva * x without any explicit loop
 
217
    if axis % 2 == 1:
 
218
        y   = N.dot(inva, (x-mu))
 
219
        y   = -0.5 * N.sum(y * (x-mu), 0)
 
220
    else:
 
221
        y   = N.dot((x-mu), inva)
 
222
        y   = -0.5 * N.sum(y * (x-mu), 1)
 
223
 
 
224
    if not log:
 
225
        y   = fac * N.exp(y)
 
226
    else:
 
227
        y   = y + N.log(fac)
 
228
 
 
229
    return y
 
230
 
 
231
if __name__ == "__main__":
 
232
    import pylab
 
233
 
 
234
    #=========================================
 
235
    # Test plotting a simple diag 2d variance:
 
236
    #=========================================
 
237
    va  = N.array([5, 3])
 
238
    mu  = N.array([2, 3])
 
239
 
 
240
    # Generate a multivariate gaussian of mean mu and covariance va
 
241
    X       = randn(2, 1e3)
 
242
    Yc      = N.dot(N.diag(N.sqrt(va)), X)
 
243
    Yc      = Yc.transpose() + mu
 
244
 
 
245
    # Plotting
 
246
    Xe, Ye  = gauss_ell(mu, va, npoints = 100)
 
247
    pylab.figure()
 
248
    pylab.plot(Yc[:, 0], Yc[:, 1], '.')
 
249
    pylab.plot(Xe, Ye, 'r')
 
250
 
 
251
    #=========================================
 
252
    # Test plotting a simple full 2d variance:
 
253
    #=========================================
 
254
    va  = N.array([[0.2, 0.1],[0.1, 0.5]])
 
255
    mu  = N.array([0, 3])
 
256
 
 
257
    # Generate a multivariate gaussian of mean mu and covariance va
 
258
    X       = randn(1e3, 2)
 
259
    Yc      = N.dot(lin.cholesky(va), X.transpose())
 
260
    Yc      = Yc.transpose() + mu
 
261
 
 
262
    # Plotting
 
263
    Xe, Ye  = gauss_ell(mu, va, npoints = 100, level=0.95)
 
264
    pylab.figure()
 
265
    pylab.plot(Yc[:, 0], Yc[:, 1], '.')
 
266
    pylab.plot(Xe, Ye, 'r')
 
267
    pylab.show()