source: orange/Orange/clustering/mixture.py @ 9976:e80c153a825a

Revision 9976:e80c153a825a, 14.1 KB checked in by ales_erjavec, 2 years ago (diff)

weightId -> weight_id

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.preprocess import Preprocessor_impute, DomainContinuizer
293#        data = Preprocessor_impute(data)
294        dc = DomainContinuizer()
295        dc.multinomial_treatment = DomainContinuizer.AsOrdinal
296        dc.continuous_treatment = DomainContinuizer.NormalizeByVariance
297        dc.class_treatment = DomainContinuizer.Ignore
298        domain = dc(data)
299        data = data.translate(domain)
300       
301        means, correlations = self.init_function(data, self.n)
302        means = numpy.asmatrix(means)
303        array, _, _ = data.to_numpy_MA()
304#        avg = numpy.mean(array, axis=0)
305#        array -= avg.reshape((1, -1))
306#        means -= avg.reshape((1, -1))
307#        std = numpy.std(array, axis=0)
308#        array /= std.reshape((1, -1))
309#        means /= std.reshape((1, -1))
310        solver = EMSolver(array, numpy.ones((self.n)) / self.n,
311                          means, correlations)
312        solver.run()
313        norm_model = GMModel(solver.weights, solver.means, solver.covariances)
314        return GMClusterModel(domain, norm_model)
315   
316       
317class GMClusterModel(object):
318    """
319    """
320    def __init__(self, domain, norm_model):
321        self.domain = domain
322        self.norm_model = norm_model
323        self.cluster_vars = [Orange.feature.Continuous("cluster %i" % i)\
324                             for i in range(len(norm_model))]
325        self.weights = self.norm_model.weights
326        self.means = self.norm_model.means
327        self.covariances = self.norm_model.covariances
328        self.inv_covariances = self.norm_model.inv_covariances
329        self.cov_determinants = self.norm_model.cov_determinants
330       
331    def __call__(self, instance, *args):
332        data = Orange.data.Table(self.domain, [instance])
333        data,_,_ = data.to_numpy_MA()
334#        print data
335       
336        p = prob_est(data, self.norm_model.weights,
337                     self.norm_model.means,
338                     self.norm_model.covariances,
339                     self.norm_model.inv_covariances,
340                     self.norm_model.cov_determinants)
341#        print p
342        p /= numpy.sum(p)
343        vals = []
344        for p, var in zip(p[0], self.cluster_vars):
345            vals.append(var(p))
346        return vals
347       
348       
349def plot_model(data_array, mixture, axis=(0, 1), samples=20, contour_lines=20):
350    """ Plot the scaterplot of data_array and the contour lines of the
351    probability for the mixture.
352     
353    """
354    import matplotlib
355    import matplotlib.pylab as plt
356    import matplotlib.cm as cm
357   
358    axis = list(axis)
359   
360    if isinstance(mixture, GMClusterModel):
361        mixture = mixture.norm_model
362   
363    if isinstance(data_array, Orange.data.Table):
364        data_array, _, _ = data_array.to_numpy_MA()
365    array = data_array[:, axis]
366   
367    weights = mixture.weights
368    means = mixture.means[:, axis]
369   
370    covariances = [cov[axis,:][:, axis] for cov in mixture.covariances] # TODO: Need the whole marginal distribution.
371   
372    gmm = GMModel(weights, means, covariances)
373   
374    min = numpy.min(array, 0)
375    max = numpy.max(array, 0)
376    extent = (min[0], max[0], min[1], max[1])
377   
378    X = numpy.linspace(min[0], max[0], samples)
379    Y = numpy.linspace(min[1], max[1], samples)
380   
381    Z = numpy.zeros((X.shape[0], Y.shape[0]))
382   
383    for i, x in enumerate(X):
384        for j, y in enumerate(Y):
385            Z[i, j] = gmm([x, y])
386           
387    plt.plot(array[:,0], array[:,1], "ro")
388    plt.contour(X, Y, Z.T, contour_lines,
389                extent=extent)
390   
391    im = plt.imshow(Z.T, interpolation='bilinear', origin='lower',
392                cmap=cm.gray, extent=extent)
393   
394    plt.plot(means[:, 0], means[:, 1], "b+")
395    plt.show()
396   
397def test(seed=0):
398#    data = Orange.data.Table(os.path.expanduser("brown-selected.tab"))
399#    data = Orange.data.Table(os.path.expanduser("~/Documents/brown-selected-fss.tab"))
400#    data = Orange.data.Table(os.path.expanduser("~/Documents/brown-selected-fss-1.tab"))
401#    data = Orange.data.Table(os.path.expanduser("~/ozone1"))
402    data = Orange.data.Table("iris.tab")
403#    data = Orange.data.Table(Orange.data.Domain(data.domain[:2], None), data)
404    numpy.random.seed(seed)
405    random.seed(seed)
406    gmm = GaussianMixture(data, n=3, init_function=init_kmeans)
407    data = data.translate(gmm.domain)
408    plot_model(data, gmm, axis=(0, 2), samples=40, contour_lines=100)
409
410   
411   
412if __name__ == "__main__":
413    test()
414   
415   
Note: See TracBrowser for help on using the repository browser.