# source:orange/Orange/clustering/mixture.py@10542:7dde0640e266

Revision 10542:7dde0640e266, 14.2 KB checked in by anzeh <anze.staric@…>, 2 years ago (diff)

Moved preprocess from Orange to Orange.data.

Line
1"""
2*******************************
3Gaussian Mixtures (``mixture``)
4*******************************
5
6This module implements a Gaussian mixture model.
7
8Example ::
9
10    >>> mixture = GaussianMixture(data, n=3)
11    >>> print mixture.means
12    >>> print mixture.weights
13    >>> print mixture.covariances
14    >>> plot_model(data, mixture, samples=40)
15
16"""
17
18import sys, os
19import numpy
20import random
21import Orange
22
23class GMModel(object):
24    """ Gaussian mixture model
25    """
26    def __init__(self, weights, means, covariances, inv_covariances=None,
27                 cov_determinants=None):
28        self.weights = weights
29        self.means = means
30        self.covariances = covariances
31        if inv_covariances is None:
32            self.inv_covariances = [numpy.linalg.pinv(cov) for cov in covariances]
33        else:
34            self.inv_covariances = inv_covariances
35
36        if cov_determinants is None:
37            self.cov_determinants = [numpy.linalg.det(cov) for cov in covariances]
38        else:
39            self.cov_determinants = cov_determinants
40
41
42    def __call__(self, instance):
43        """ Return the probability of instance.
44        """
45        return numpy.sum(prob_est([instance], self.weights, self.means,
46                                  self.covariances,
47                                  self.inv_covariances,
48                                  self.cov_determinants))
49
50    def __getitem__(self, index):
51        """ Return the index-th gaussian.
52        """
53        return GMModel([1.0], self.means[index: index + 1],
54                       self.covariances[index: index + 1],
55                       self.inv_covariances[index: index + 1],
56                       self.cov_determinants[index: index + 1])
57
58    def __len__(self):
59        return len(self.weights)
60
61
62def init_random(data, n, *args, **kwargs):
63    """ Init random means and correlations from a data table.
64
65    :param data: data table
66    :type data: :class:`Orange.data.Table`
67    :param n: Number of centers and correlations to return.
68    :type n: int
69
70    """
71    if isinstance(data, Orange.data.Table):
72        array, w, c = data.toNumpyMA()
73    else:
74        array = numpy.asarray(data)
75
76    min, max = array.max(0), array.min(0)
77    dim = array.shape[1]
78    means = numpy.zeros((n, dim))
79    for i in range(n):
80        means[i] = [numpy.random.uniform(low, high) for low, high in zip(min, max)]
81
82    correlations = [numpy.asmatrix(numpy.eye(dim)) for i in range(n)]
83    return means, correlations
84
85def init_kmeans(data, n, *args, **kwargs):
86    """ Init with k-means algorithm.
87
88    :param data: data table
89    :type data: :class:`Orange.data.Table`
90    :param n: Number of centers and correlations to return.
91    :type n: int
92
93    """
94    if not isinstance(data, Orange.data.Table):
95        raise TypeError("Orange.data.Table instance expected!")
96    from Orange.clustering.kmeans import Clustering
97    km = Clustering(data, centroids=n, maxiters=20, nstart=3)
98    centers = Orange.data.Table(km.centroids)
99    centers, w, c = centers.toNumpyMA()
100    dim = len(data.domain.attributes)
101    correlations = [numpy.asmatrix(numpy.eye(dim)) for i in range(n)]
102    return centers, correlations
103
104def prob_est1(data, mean, covariance, inv_covariance=None, det=None):
105    """ Return the probability of data given mean and covariance matrix
106    """
107    data = numpy.asmatrix(data)
108    mean = numpy.asmatrix(mean)
109    if inv_covariance is None:
110        inv_covariance = numpy.linalg.pinv(covariance)
111
112    inv_covariance = numpy.asmatrix(inv_covariance)
113
114    diff = data - mean
115    p = numpy.zeros(data.shape[0])
116    for i in range(data.shape[0]):
117        d = diff[i]
118        p[i] = d * inv_covariance * d.T
119
120    p *= -0.5
121    p = numpy.exp(p)
122    p /= numpy.power(2 * numpy.pi, numpy.rank(covariance) / 2.0)
123    if det is None:
124        det = numpy.linalg.det(covariance)
125    assert(det != 0.0)
126    p /= det
127    return p
128
129
130def prob_est(data, weights, means, covariances, inv_covariances=None, cov_determinants=None):
131    """ Return the probability estimation of data for each
132    gausian given weights, means and covariances.
133
134    """
135    if inv_covariances is None:
136        inv_covariances = [numpy.linalg.pinv(cov) for cov in covariances]
137
138    if cov_determinants is None:
139        cov_determinants = [numpy.linalg.det(cov) for cov in covariances]
140
141    data = numpy.asmatrix(data)
142    probs = numpy.zeros((data.shape[0], len(weights)))
143
144    for i, (w, mean, cov, inv_cov, det) in enumerate(zip(weights, means,
145                                        covariances, inv_covariances,
146                                        cov_determinants)):
147        probs[:, i] = w * prob_est1(data, mean, cov, inv_cov, det)
148
149    return probs
150
151
152class EMSolver(object):
153    """ An EM solver for gaussian mixture model
154    """
155    _TRACE_MEAN = False
156    def __init__(self, data, weights, means, covariances):
157        self.data = data
158        self.weights = weights
159        self.means = means
160        self.covariances = covariances
161        self.inv_covariances = [numpy.matrix(numpy.linalg.pinv(cov)) for cov in covariances]
162        self.cov_determinants = [numpy.linalg.det(cov) for cov in self.covariances]
163        self.n_clusters = len(self.weights)
164        self.data_dim = self.data.shape[1]
165
166        self.probs = prob_est(data, weights, means, covariances,
167                              self.inv_covariances, self.cov_determinants)
168
169        self.log_likelihood = self._log_likelihood()
170#        print "W", self.weights
171#        print "P", self.probs
172#        print "L", self.log_likelihood
173#        print "C", self.covariances
174#        print "Det", self.cov_determinants
175
176    def _log_likelihood(self):
177        """ Compute the log likelihood of the current solution.
178        """
179        log_p = numpy.log(numpy.sum(self.weights * self.probs, axis=0))
180        return 1.0 / len(self.data) * numpy.sum(log_p)
181
182    def E_step(self):
183        """ E Step
184        """
185        self.probs = prob_est(self.data, self.weights, self.means,
186                         self.covariances, self.inv_covariances)
187
188#        print "PPP", self.probs
189#        print "P sum", numpy.sum(self.probs, axis=1).reshape((-1, 1))
190        self.probs /= numpy.sum(self.probs, axis=1).reshape((-1, 1))
191
192        # Update the Q
193#        self.Q = 0.0
194#        prob_sum = numpy.sum(self.probs, axis=0)
195#        self.Q = sum([p*(numpy.log(w) - 0.5 * numpy.linalg.det(cov)) \
196#                      for p, w, cov in zip(prob_sum, self.weights,
197#                                           self.covariances)]) * \
198#                      len(self.data)
199#
200#        for i in range(len(data)):
201#            for j in range(self.n_clusters):
202#                diff = numpy.asmatrix(self.data[i] - self.means[j])
203#                self.Q += - 0.5 * self.probs[i, j] * diff.T * self.inv_covariances[j] * diff
204#        print self.Q
205
206        self.log_likelihood = self._log_likelihood()
207#        print "Prob:", self.probs
208#        print "Log like.:", self.log_likelihood
209
210    def M_step(self):
211        """ M step
212        """
213        # Compute the new weights
214        prob_sum = numpy.sum(self.probs, axis=0)
215        weights = prob_sum / numpy.sum(prob_sum)
216
217        # Compute the new means
218        means = []
219        for j in range(self.n_clusters):
220            means.append(numpy.sum(self.probs[:, j].reshape((-1, 1)) * self.data, axis=0) /  prob_sum[j])
221
222        # Compute the new covariances
223        covariances = []
224        cov_determinants = []
225        for j in range(self.n_clusters):
226            cov = numpy.zeros(self.covariances[j].shape)
227            diff = self.data - means[j]
228            diff = numpy.asmatrix(diff)
229            for i in range(len(self.data)): # TODO: speed up
230                cov += self.probs[i, j] * diff[i].T * diff[i]
231
232            cov *= 1.0 / prob_sum[j]
233            det = numpy.linalg.det(cov)
234
235            covariances.append(numpy.asmatrix(cov))
236            cov_determinants.append(det)
237#            self.inv_covariances[j] = numpy.linalg.pinv(cov)
238#            self.cov_determinants[j] = det
239        self.weights = weights
240        self.means = numpy.asmatrix(means)
241        self.covariances = covariances
242        self.cov_determinants = cov_determinants
243
244    def one_step(self):
245        """ One iteration of both M and E step.
246        """
247        self.E_step()
248        self.M_step()
249
250    def run(self, max_iter=sys.maxint, eps=1e-5):
251        """ Run the EM algorithm.
252        """
253        if self._TRACE_MEAN:
254            from pylab import plot, show, draw, ion
255            ion()
256            plot(self.data[:, 0], self.data[:, 1], "ro")
257            vec_plot = plot(self.means[:, 0], self.means[:, 1], "bo")[0]
258
259        curr_iter = 0
260
261        while True:
262            old_objective = self.log_likelihood
263            self.one_step()
264
265            if self._TRACE_MEAN:
266                vec_plot.set_xdata(self.means[:, 0])
267                vec_plot.set_ydata(self.means[:, 1])
268                draw()
269
270            curr_iter += 1
271#            print curr_iter
272#            print abs(old_objective - self.log_likelihood)
273            if abs(old_objective - self.log_likelihood) < eps or curr_iter > max_iter:
274                break
275
276class GaussianMixture(object):
277    """ Computes the gaussian mixture model from an Orange data-set.
278    """
279    def __new__(cls, data=None, weight_id=None, **kwargs):
280        self = object.__new__(cls)
281        if data is not None:
282            self.__init__(**kwargs)
283            return self.__call__(data, weight_id)
284        else:
285            return self
286
287    def __init__(self, n=3, init_function=init_kmeans):
288        self.n = n
289        self.init_function = init_function
290
291    def __call__(self, data, weight_id=None):
292        from Orange.data import preprocess
293        #import Preprocessor_impute, DomainContinuizer
294#        data = Preprocessor_impute(data)
295        dc = preprocess.DomainContinuizer()
296        dc.multinomial_treatment = preprocess.DomainContinuizer.AsOrdinal
297        dc.continuous_treatment = preprocess.DomainContinuizer.NormalizeByVariance
298        dc.class_treatment = preprocess.DomainContinuizer.Ignore
299        domain = dc(data)
300        data = data.translate(domain)
301
302        means, correlations = self.init_function(data, self.n)
303        means = numpy.asmatrix(means)
304        array, _, _ = data.to_numpy_MA()
305#        avg = numpy.mean(array, axis=0)
306#        array -= avg.reshape((1, -1))
307#        means -= avg.reshape((1, -1))
308#        std = numpy.std(array, axis=0)
309#        array /= std.reshape((1, -1))
310#        means /= std.reshape((1, -1))
311        solver = EMSolver(array, numpy.ones(self.n) / self.n,
312                          means, correlations)
313        solver.run()
314        norm_model = GMModel(solver.weights, solver.means, solver.covariances)
315        return GMClusterModel(domain, norm_model)
316
317
318class GMClusterModel(object):
319    """
320    """
321    def __init__(self, domain, norm_model):
322        self.domain = domain
323        self.norm_model = norm_model
324        self.cluster_vars = [Orange.feature.Continuous("cluster %i" % i)\
325                             for i in range(len(norm_model))]
326        self.weights = self.norm_model.weights
327        self.means = self.norm_model.means
328        self.covariances = self.norm_model.covariances
329        self.inv_covariances = self.norm_model.inv_covariances
330        self.cov_determinants = self.norm_model.cov_determinants
331
332    def __call__(self, instance, *args):
333        data = Orange.data.Table(self.domain, [instance])
334        data,_,_ = data.to_numpy_MA()
335#        print data
336
337        p = prob_est(data, self.norm_model.weights,
338                     self.norm_model.means,
339                     self.norm_model.covariances,
340                     self.norm_model.inv_covariances,
341                     self.norm_model.cov_determinants)
342#        print p
343        p /= numpy.sum(p)
344        vals = []
345        for p, var in zip(p[0], self.cluster_vars):
346            vals.append(var(p))
347        return vals
348
349
350def plot_model(data_array, mixture, axis=(0, 1), samples=20, contour_lines=20):
351    """ Plot the scaterplot of data_array and the contour lines of the
352    probability for the mixture.
353
354    """
355    import matplotlib
356    import matplotlib.pylab as plt
357    import matplotlib.cm as cm
358
359    axis = list(axis)
360
361    if isinstance(mixture, GMClusterModel):
362        mixture = mixture.norm_model
363
364    if isinstance(data_array, Orange.data.Table):
365        data_array, _, _ = data_array.to_numpy_MA()
366    array = data_array[:, axis]
367
368    weights = mixture.weights
369    means = mixture.means[:, axis]
370
371    covariances = [cov[axis,:][:, axis] for cov in mixture.covariances] # TODO: Need the whole marginal distribution.
372
373    gmm = GMModel(weights, means, covariances)
374
375    min = numpy.min(array, 0)
376    max = numpy.max(array, 0)
377    extent = (min[0], max[0], min[1], max[1])
378
379    X = numpy.linspace(min[0], max[0], samples)
380    Y = numpy.linspace(min[1], max[1], samples)
381
382    Z = numpy.zeros((X.shape[0], Y.shape[0]))
383
384    for i, x in enumerate(X):
385        for j, y in enumerate(Y):
386            Z[i, j] = gmm([x, y])
387
388    plt.plot(array[:,0], array[:,1], "ro")
389    plt.contour(X, Y, Z.T, contour_lines,
390                extent=extent)
391
392    im = plt.imshow(Z.T, interpolation='bilinear', origin='lower',
393                cmap=cm.gray, extent=extent)
394
395    plt.plot(means[:, 0], means[:, 1], "b+")
396    plt.show()
397
398def test(seed=0):
399#    data = Orange.data.Table(os.path.expanduser("brown-selected.tab"))
400#    data = Orange.data.Table(os.path.expanduser("~/Documents/brown-selected-fss.tab"))
401#    data = Orange.data.Table(os.path.expanduser("~/Documents/brown-selected-fss-1.tab"))
402#    data = Orange.data.Table(os.path.expanduser("~/ozone1"))
403    data = Orange.data.Table("iris.tab")
404#    data = Orange.data.Table(Orange.data.Domain(data.domain[:2], None), data)
405    numpy.random.seed(seed)
406    random.seed(seed)
407    gmm = GaussianMixture(data, n=3, init_function=init_kmeans)
408    data = data.translate(gmm.domain)
409    plot_model(data, gmm, axis=(0, 2), samples=40, contour_lines=100)
410
411
412
413if __name__ == "__main__":
414    test()
415
416
Note: See TracBrowser for help on using the repository browser.