# 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.