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.