_base.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
  1. """Base class for mixture models."""
  2. # Author: Wei Xue <xuewei4d@gmail.com>
  3. # Modified by Thierry Guillemot <thierry.guillemot.work@gmail.com>
  4. # License: BSD 3 clause
  5. import warnings
  6. from abc import ABCMeta, abstractmethod
  7. from numbers import Integral, Real
  8. from time import time
  9. import numpy as np
  10. from scipy.special import logsumexp
  11. from .. import cluster
  12. from ..base import BaseEstimator, DensityMixin, _fit_context
  13. from ..cluster import kmeans_plusplus
  14. from ..exceptions import ConvergenceWarning
  15. from ..utils import check_random_state
  16. from ..utils._param_validation import Interval, StrOptions
  17. from ..utils.validation import check_is_fitted
  18. def _check_shape(param, param_shape, name):
  19. """Validate the shape of the input parameter 'param'.
  20. Parameters
  21. ----------
  22. param : array
  23. param_shape : tuple
  24. name : str
  25. """
  26. param = np.array(param)
  27. if param.shape != param_shape:
  28. raise ValueError(
  29. "The parameter '%s' should have the shape of %s, but got %s"
  30. % (name, param_shape, param.shape)
  31. )
  32. class BaseMixture(DensityMixin, BaseEstimator, metaclass=ABCMeta):
  33. """Base class for mixture models.
  34. This abstract class specifies an interface for all mixture classes and
  35. provides basic common methods for mixture models.
  36. """
  37. _parameter_constraints: dict = {
  38. "n_components": [Interval(Integral, 1, None, closed="left")],
  39. "tol": [Interval(Real, 0.0, None, closed="left")],
  40. "reg_covar": [Interval(Real, 0.0, None, closed="left")],
  41. "max_iter": [Interval(Integral, 0, None, closed="left")],
  42. "n_init": [Interval(Integral, 1, None, closed="left")],
  43. "init_params": [
  44. StrOptions({"kmeans", "random", "random_from_data", "k-means++"})
  45. ],
  46. "random_state": ["random_state"],
  47. "warm_start": ["boolean"],
  48. "verbose": ["verbose"],
  49. "verbose_interval": [Interval(Integral, 1, None, closed="left")],
  50. }
  51. def __init__(
  52. self,
  53. n_components,
  54. tol,
  55. reg_covar,
  56. max_iter,
  57. n_init,
  58. init_params,
  59. random_state,
  60. warm_start,
  61. verbose,
  62. verbose_interval,
  63. ):
  64. self.n_components = n_components
  65. self.tol = tol
  66. self.reg_covar = reg_covar
  67. self.max_iter = max_iter
  68. self.n_init = n_init
  69. self.init_params = init_params
  70. self.random_state = random_state
  71. self.warm_start = warm_start
  72. self.verbose = verbose
  73. self.verbose_interval = verbose_interval
  74. @abstractmethod
  75. def _check_parameters(self, X):
  76. """Check initial parameters of the derived class.
  77. Parameters
  78. ----------
  79. X : array-like of shape (n_samples, n_features)
  80. """
  81. pass
  82. def _initialize_parameters(self, X, random_state):
  83. """Initialize the model parameters.
  84. Parameters
  85. ----------
  86. X : array-like of shape (n_samples, n_features)
  87. random_state : RandomState
  88. A random number generator instance that controls the random seed
  89. used for the method chosen to initialize the parameters.
  90. """
  91. n_samples, _ = X.shape
  92. if self.init_params == "kmeans":
  93. resp = np.zeros((n_samples, self.n_components))
  94. label = (
  95. cluster.KMeans(
  96. n_clusters=self.n_components, n_init=1, random_state=random_state
  97. )
  98. .fit(X)
  99. .labels_
  100. )
  101. resp[np.arange(n_samples), label] = 1
  102. elif self.init_params == "random":
  103. resp = random_state.uniform(size=(n_samples, self.n_components))
  104. resp /= resp.sum(axis=1)[:, np.newaxis]
  105. elif self.init_params == "random_from_data":
  106. resp = np.zeros((n_samples, self.n_components))
  107. indices = random_state.choice(
  108. n_samples, size=self.n_components, replace=False
  109. )
  110. resp[indices, np.arange(self.n_components)] = 1
  111. elif self.init_params == "k-means++":
  112. resp = np.zeros((n_samples, self.n_components))
  113. _, indices = kmeans_plusplus(
  114. X,
  115. self.n_components,
  116. random_state=random_state,
  117. )
  118. resp[indices, np.arange(self.n_components)] = 1
  119. self._initialize(X, resp)
  120. @abstractmethod
  121. def _initialize(self, X, resp):
  122. """Initialize the model parameters of the derived class.
  123. Parameters
  124. ----------
  125. X : array-like of shape (n_samples, n_features)
  126. resp : array-like of shape (n_samples, n_components)
  127. """
  128. pass
  129. def fit(self, X, y=None):
  130. """Estimate model parameters with the EM algorithm.
  131. The method fits the model ``n_init`` times and sets the parameters with
  132. which the model has the largest likelihood or lower bound. Within each
  133. trial, the method iterates between E-step and M-step for ``max_iter``
  134. times until the change of likelihood or lower bound is less than
  135. ``tol``, otherwise, a ``ConvergenceWarning`` is raised.
  136. If ``warm_start`` is ``True``, then ``n_init`` is ignored and a single
  137. initialization is performed upon the first call. Upon consecutive
  138. calls, training starts where it left off.
  139. Parameters
  140. ----------
  141. X : array-like of shape (n_samples, n_features)
  142. List of n_features-dimensional data points. Each row
  143. corresponds to a single data point.
  144. y : Ignored
  145. Not used, present for API consistency by convention.
  146. Returns
  147. -------
  148. self : object
  149. The fitted mixture.
  150. """
  151. # parameters are validated in fit_predict
  152. self.fit_predict(X, y)
  153. return self
  154. @_fit_context(prefer_skip_nested_validation=True)
  155. def fit_predict(self, X, y=None):
  156. """Estimate model parameters using X and predict the labels for X.
  157. The method fits the model n_init times and sets the parameters with
  158. which the model has the largest likelihood or lower bound. Within each
  159. trial, the method iterates between E-step and M-step for `max_iter`
  160. times until the change of likelihood or lower bound is less than
  161. `tol`, otherwise, a :class:`~sklearn.exceptions.ConvergenceWarning` is
  162. raised. After fitting, it predicts the most probable label for the
  163. input data points.
  164. .. versionadded:: 0.20
  165. Parameters
  166. ----------
  167. X : array-like of shape (n_samples, n_features)
  168. List of n_features-dimensional data points. Each row
  169. corresponds to a single data point.
  170. y : Ignored
  171. Not used, present for API consistency by convention.
  172. Returns
  173. -------
  174. labels : array, shape (n_samples,)
  175. Component labels.
  176. """
  177. X = self._validate_data(X, dtype=[np.float64, np.float32], ensure_min_samples=2)
  178. if X.shape[0] < self.n_components:
  179. raise ValueError(
  180. "Expected n_samples >= n_components "
  181. f"but got n_components = {self.n_components}, "
  182. f"n_samples = {X.shape[0]}"
  183. )
  184. self._check_parameters(X)
  185. # if we enable warm_start, we will have a unique initialisation
  186. do_init = not (self.warm_start and hasattr(self, "converged_"))
  187. n_init = self.n_init if do_init else 1
  188. max_lower_bound = -np.inf
  189. self.converged_ = False
  190. random_state = check_random_state(self.random_state)
  191. n_samples, _ = X.shape
  192. for init in range(n_init):
  193. self._print_verbose_msg_init_beg(init)
  194. if do_init:
  195. self._initialize_parameters(X, random_state)
  196. lower_bound = -np.inf if do_init else self.lower_bound_
  197. if self.max_iter == 0:
  198. best_params = self._get_parameters()
  199. best_n_iter = 0
  200. else:
  201. for n_iter in range(1, self.max_iter + 1):
  202. prev_lower_bound = lower_bound
  203. log_prob_norm, log_resp = self._e_step(X)
  204. self._m_step(X, log_resp)
  205. lower_bound = self._compute_lower_bound(log_resp, log_prob_norm)
  206. change = lower_bound - prev_lower_bound
  207. self._print_verbose_msg_iter_end(n_iter, change)
  208. if abs(change) < self.tol:
  209. self.converged_ = True
  210. break
  211. self._print_verbose_msg_init_end(lower_bound)
  212. if lower_bound > max_lower_bound or max_lower_bound == -np.inf:
  213. max_lower_bound = lower_bound
  214. best_params = self._get_parameters()
  215. best_n_iter = n_iter
  216. # Should only warn about convergence if max_iter > 0, otherwise
  217. # the user is assumed to have used 0-iters initialization
  218. # to get the initial means.
  219. if not self.converged_ and self.max_iter > 0:
  220. warnings.warn(
  221. "Initialization %d did not converge. "
  222. "Try different init parameters, "
  223. "or increase max_iter, tol "
  224. "or check for degenerate data." % (init + 1),
  225. ConvergenceWarning,
  226. )
  227. self._set_parameters(best_params)
  228. self.n_iter_ = best_n_iter
  229. self.lower_bound_ = max_lower_bound
  230. # Always do a final e-step to guarantee that the labels returned by
  231. # fit_predict(X) are always consistent with fit(X).predict(X)
  232. # for any value of max_iter and tol (and any random_state).
  233. _, log_resp = self._e_step(X)
  234. return log_resp.argmax(axis=1)
  235. def _e_step(self, X):
  236. """E step.
  237. Parameters
  238. ----------
  239. X : array-like of shape (n_samples, n_features)
  240. Returns
  241. -------
  242. log_prob_norm : float
  243. Mean of the logarithms of the probabilities of each sample in X
  244. log_responsibility : array, shape (n_samples, n_components)
  245. Logarithm of the posterior probabilities (or responsibilities) of
  246. the point of each sample in X.
  247. """
  248. log_prob_norm, log_resp = self._estimate_log_prob_resp(X)
  249. return np.mean(log_prob_norm), log_resp
  250. @abstractmethod
  251. def _m_step(self, X, log_resp):
  252. """M step.
  253. Parameters
  254. ----------
  255. X : array-like of shape (n_samples, n_features)
  256. log_resp : array-like of shape (n_samples, n_components)
  257. Logarithm of the posterior probabilities (or responsibilities) of
  258. the point of each sample in X.
  259. """
  260. pass
  261. @abstractmethod
  262. def _get_parameters(self):
  263. pass
  264. @abstractmethod
  265. def _set_parameters(self, params):
  266. pass
  267. def score_samples(self, X):
  268. """Compute the log-likelihood of each sample.
  269. Parameters
  270. ----------
  271. X : array-like of shape (n_samples, n_features)
  272. List of n_features-dimensional data points. Each row
  273. corresponds to a single data point.
  274. Returns
  275. -------
  276. log_prob : array, shape (n_samples,)
  277. Log-likelihood of each sample in `X` under the current model.
  278. """
  279. check_is_fitted(self)
  280. X = self._validate_data(X, reset=False)
  281. return logsumexp(self._estimate_weighted_log_prob(X), axis=1)
  282. def score(self, X, y=None):
  283. """Compute the per-sample average log-likelihood of the given data X.
  284. Parameters
  285. ----------
  286. X : array-like of shape (n_samples, n_dimensions)
  287. List of n_features-dimensional data points. Each row
  288. corresponds to a single data point.
  289. y : Ignored
  290. Not used, present for API consistency by convention.
  291. Returns
  292. -------
  293. log_likelihood : float
  294. Log-likelihood of `X` under the Gaussian mixture model.
  295. """
  296. return self.score_samples(X).mean()
  297. def predict(self, X):
  298. """Predict the labels for the data samples in X using trained model.
  299. Parameters
  300. ----------
  301. X : array-like of shape (n_samples, n_features)
  302. List of n_features-dimensional data points. Each row
  303. corresponds to a single data point.
  304. Returns
  305. -------
  306. labels : array, shape (n_samples,)
  307. Component labels.
  308. """
  309. check_is_fitted(self)
  310. X = self._validate_data(X, reset=False)
  311. return self._estimate_weighted_log_prob(X).argmax(axis=1)
  312. def predict_proba(self, X):
  313. """Evaluate the components' density for each sample.
  314. Parameters
  315. ----------
  316. X : array-like of shape (n_samples, n_features)
  317. List of n_features-dimensional data points. Each row
  318. corresponds to a single data point.
  319. Returns
  320. -------
  321. resp : array, shape (n_samples, n_components)
  322. Density of each Gaussian component for each sample in X.
  323. """
  324. check_is_fitted(self)
  325. X = self._validate_data(X, reset=False)
  326. _, log_resp = self._estimate_log_prob_resp(X)
  327. return np.exp(log_resp)
  328. def sample(self, n_samples=1):
  329. """Generate random samples from the fitted Gaussian distribution.
  330. Parameters
  331. ----------
  332. n_samples : int, default=1
  333. Number of samples to generate.
  334. Returns
  335. -------
  336. X : array, shape (n_samples, n_features)
  337. Randomly generated sample.
  338. y : array, shape (nsamples,)
  339. Component labels.
  340. """
  341. check_is_fitted(self)
  342. if n_samples < 1:
  343. raise ValueError(
  344. "Invalid value for 'n_samples': %d . The sampling requires at "
  345. "least one sample." % (self.n_components)
  346. )
  347. _, n_features = self.means_.shape
  348. rng = check_random_state(self.random_state)
  349. n_samples_comp = rng.multinomial(n_samples, self.weights_)
  350. if self.covariance_type == "full":
  351. X = np.vstack(
  352. [
  353. rng.multivariate_normal(mean, covariance, int(sample))
  354. for (mean, covariance, sample) in zip(
  355. self.means_, self.covariances_, n_samples_comp
  356. )
  357. ]
  358. )
  359. elif self.covariance_type == "tied":
  360. X = np.vstack(
  361. [
  362. rng.multivariate_normal(mean, self.covariances_, int(sample))
  363. for (mean, sample) in zip(self.means_, n_samples_comp)
  364. ]
  365. )
  366. else:
  367. X = np.vstack(
  368. [
  369. mean
  370. + rng.standard_normal(size=(sample, n_features))
  371. * np.sqrt(covariance)
  372. for (mean, covariance, sample) in zip(
  373. self.means_, self.covariances_, n_samples_comp
  374. )
  375. ]
  376. )
  377. y = np.concatenate(
  378. [np.full(sample, j, dtype=int) for j, sample in enumerate(n_samples_comp)]
  379. )
  380. return (X, y)
  381. def _estimate_weighted_log_prob(self, X):
  382. """Estimate the weighted log-probabilities, log P(X | Z) + log weights.
  383. Parameters
  384. ----------
  385. X : array-like of shape (n_samples, n_features)
  386. Returns
  387. -------
  388. weighted_log_prob : array, shape (n_samples, n_component)
  389. """
  390. return self._estimate_log_prob(X) + self._estimate_log_weights()
  391. @abstractmethod
  392. def _estimate_log_weights(self):
  393. """Estimate log-weights in EM algorithm, E[ log pi ] in VB algorithm.
  394. Returns
  395. -------
  396. log_weight : array, shape (n_components, )
  397. """
  398. pass
  399. @abstractmethod
  400. def _estimate_log_prob(self, X):
  401. """Estimate the log-probabilities log P(X | Z).
  402. Compute the log-probabilities per each component for each sample.
  403. Parameters
  404. ----------
  405. X : array-like of shape (n_samples, n_features)
  406. Returns
  407. -------
  408. log_prob : array, shape (n_samples, n_component)
  409. """
  410. pass
  411. def _estimate_log_prob_resp(self, X):
  412. """Estimate log probabilities and responsibilities for each sample.
  413. Compute the log probabilities, weighted log probabilities per
  414. component and responsibilities for each sample in X with respect to
  415. the current state of the model.
  416. Parameters
  417. ----------
  418. X : array-like of shape (n_samples, n_features)
  419. Returns
  420. -------
  421. log_prob_norm : array, shape (n_samples,)
  422. log p(X)
  423. log_responsibilities : array, shape (n_samples, n_components)
  424. logarithm of the responsibilities
  425. """
  426. weighted_log_prob = self._estimate_weighted_log_prob(X)
  427. log_prob_norm = logsumexp(weighted_log_prob, axis=1)
  428. with np.errstate(under="ignore"):
  429. # ignore underflow
  430. log_resp = weighted_log_prob - log_prob_norm[:, np.newaxis]
  431. return log_prob_norm, log_resp
  432. def _print_verbose_msg_init_beg(self, n_init):
  433. """Print verbose message on initialization."""
  434. if self.verbose == 1:
  435. print("Initialization %d" % n_init)
  436. elif self.verbose >= 2:
  437. print("Initialization %d" % n_init)
  438. self._init_prev_time = time()
  439. self._iter_prev_time = self._init_prev_time
  440. def _print_verbose_msg_iter_end(self, n_iter, diff_ll):
  441. """Print verbose message on initialization."""
  442. if n_iter % self.verbose_interval == 0:
  443. if self.verbose == 1:
  444. print(" Iteration %d" % n_iter)
  445. elif self.verbose >= 2:
  446. cur_time = time()
  447. print(
  448. " Iteration %d\t time lapse %.5fs\t ll change %.5f"
  449. % (n_iter, cur_time - self._iter_prev_time, diff_ll)
  450. )
  451. self._iter_prev_time = cur_time
  452. def _print_verbose_msg_init_end(self, ll):
  453. """Print verbose message on the end of iteration."""
  454. if self.verbose == 1:
  455. print("Initialization converged: %s" % self.converged_)
  456. elif self.verbose >= 2:
  457. print(
  458. "Initialization converged: %s\t time lapse %.5fs\t ll %.5f"
  459. % (self.converged_, time() - self._init_prev_time, ll)
  460. )