hmm.py

# Hidden Markov Models
#
# Author: Ron Weiss <ronweiss@gmail.com>
#         Shiqiao Du <lucidfrontier.45@gmail.com>
# API changes: Jaques Grobler <jaquesgrobler@gmail.com>
# Modifications to create of the HMMLearn module: Gael Varoquaux
# More API changes: Sergei Lebedev <superbobry@gmail.com>
# FSHMM code: Elizabeth Fons <elifons@gmail.com>
#			  Alejandro Sztrajman <asztrajman@gmail.com>
# FSHMM algorithm from:
# Adams, Stephen & Beling, Peter & Cogill, Randy. (2016). 
# Feature Selection for Hidden Markov Models and Hidden Semi-Markov Models. 
"""
The :mod:`hmmlearn.hmm` module implements hidden Markov models.
"""

import numpy as np
import pandas as pd
from sklearn import cluster
from sklearn.mixture import (
    GMM, sample_gaussian,
    log_multivariate_normal_density,
    distribute_covar_matrix_to_match_covariance_type, _validate_covars)

from scipy.misc import logsumexp
from sklearn.base import BaseEstimator, _pprint
from sklearn.utils import check_array, check_random_state
from sklearn.utils.validation import check_is_fitted

from base import _BaseHMM
from hmmlearn.utils import normalize, iter_from_X_lengths, normalize

from base import ConvergenceMonitor

__all__ = ["GMMHMM", "GaussianHMM", "MultinomialHMM", "GaussianFSHHM"]

COVARIANCE_TYPES = frozenset(("spherical", "diag", "full", "tied"))


class GaussianHMM(_BaseHMM):
    """Hidden Markov Model with Gaussian emissions.

    Parameters
    ----------
    n_components : int
        Number of states.

    covariance_type : string
        String describing the type of covariance parameters to
        use.  Must be one of

        * "spherical" --- each state uses a single variance value that
          applies to all features;
        * "diag" --- each state uses a diagonal covariance matrix;
        * "full" --- each state uses a full (i.e. unrestricted)
          covariance matrix;
        * "tied" --- all states use **the same** full covariance matrix.

        Defaults to "diag".

    min_covar : float
        Floor on the diagonal of the covariance matrix to prevent
        overfitting. Defaults to 1e-3.

    startprob_prior : array, shape (n_components, )
        Initial state occupation prior distribution.

    transmat_prior : array, shape (n_components, n_components)
        Matrix of prior transition probabilities between states.

    algorithm : string
        Decoder algorithm. Must be one of "viterbi" or "map".
        Defaults to "viterbi".

    random_state: RandomState or an int seed
        A random number generator instance.

    n_iter : int, optional
        Maximum number of iterations to perform.

    tol : float, optional
        Convergence threshold. EM will stop if the gain in log-likelihood
        is below this value.

    verbose : bool, optional
        When ``True`` per-iteration convergence reports are printed
        to :data:`sys.stderr`. You can diagnose convergence via the
        :attr:`monitor_` attribute.

    params : string, optional
        Controls which parameters are updated in the training
        process.  Can contain any combination of 's' for startprob,
        't' for transmat, 'm' for means and 'c' for covars. Defaults
        to all parameters.

    init_params : string, optional
        Controls which parameters are initialized prior to
        training.  Can contain any combination of 's' for
        startprob, 't' for transmat, 'm' for means and 'c' for covars.
        Defaults to all parameters.

    Attributes
    ----------
    n_features : int
        Dimensionality of the Gaussian emissions.

    monitor\_ : ConvergenceMonitor
        Monitor object used to check the convergence of EM.

    transmat\_ : array, shape (n_components, n_components)
        Matrix of transition probabilities between states.

    startprob\_ : array, shape (n_components, )
        Initial state occupation distribution.

    means\_ : array, shape (n_components, n_features)
        Mean parameters for each state.

    covars\_ : array
        Covariance parameters for each state.

        The shape depends on ``covariance_type``::

            (n_components, )                        if 'spherical',
            (n_features, n_features)                if 'tied',
            (n_components, n_features)              if 'diag',
            (n_components, n_features, n_features)  if 'full'

    Examples
    --------
    >>> from hmmlearn.hmm import GaussianHMM
    >>> GaussianHMM(n_components=2)
    ...                             #doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
    GaussianHMM(algorithm='viterbi',...
    """
    def __init__(self, n_components=1, covariance_type='diag',
                 min_covar=1e-3,
                 startprob_prior=1.0, transmat_prior=1.0,
                 means_prior=0, means_weight=0,
                 covars_prior=1e-2, covars_weight=1,
                 algorithm="viterbi", random_state=None,
                 n_iter=400, tol=1e-9, verbose=False,
        
                 params="stmc", init_params="stmc"):
        _BaseHMM.__init__(self, n_components,
                          startprob_prior=startprob_prior,
                          transmat_prior=transmat_prior, algorithm=algorithm,
                          random_state=random_state, n_iter=n_iter,
                          tol=tol, params=params, verbose=verbose,
                          init_params=init_params)

        self.covariance_type = covariance_type
        self.min_covar = min_covar
        self.means_prior = means_prior
        self.means_weight = means_weight
        self.covars_prior = covars_prior
        self.covars_weight = covars_weight

    def _get_covars(self):
        """Return covars as a full matrix."""
        if self.covariance_type == 'full':
            return self._covars_
        elif self.covariance_type == 'diag':
            return np.array([np.diag(cov) for cov in self._covars_])
        elif self.covariance_type == 'tied':
            return np.array([self._covars_] * self.n_components)
        elif self.covariance_type == 'spherical':
            return np.array(
                [np.eye(self.n_features) * cov for cov in self._covars_])

    def _set_covars(self, covars):
        self._covars_ = np.asarray(covars).copy()
        covars_ = property(_get_covars, _set_covars)

    def _check(self):
        super(GaussianHMM, self)._check()

        self.means_ = np.asarray(self.means_)
        self.n_features = self.means_.shape[1]
        if self.covariance_type not in COVARIANCE_TYPES:
            raise ValueError('covariance_type must be one of {0}'
                             .format(COVARIANCE_TYPES))
        _validate_covars(self._covars_, self.covariance_type,
                         self.n_components)

    def _init(self, X, lengths=None):
        super(GaussianHMM, self)._init(X, lengths=lengths)

        _, n_features = X.shape
        if hasattr(self, 'n_features') and self.n_features != n_features:
            raise ValueError('Unexpected number of dimensions, got %s but '
                             'expected %s' % (n_features, self.n_features))

        self.n_features = n_features
        if 'm' in self.init_params or not hasattr(self, "means_"): #paper: mu initialized randomly
            kmeans = cluster.KMeans(n_clusters=self.n_components,
                                    random_state=self.random_state)
            kmeans.fit(X)
            self.means_ = kmeans.cluster_centers_
        if 'c' in self.init_params or not hasattr(self, "covars_"): #paper: sigma initialized with 4
            cv = np.cov(X.T) + self.min_covar * np.eye(X.shape[1])
            if not cv.shape:
                cv.shape = (1, 1)
            self._covars_ = distribute_covar_matrix_to_match_covariance_type(
                cv, self.covariance_type, self.n_components).copy()
            self._covars_ = np.ones(self._covars_.shape)*4.0 #villavilla

    def _compute_log_likelihood(self, X):
        return log_multivariate_normal_density(
            X, self.means_, self._covars_, self.covariance_type)

    def _generate_sample_from_state(self, state, random_state=None):
        if self.covariance_type == 'tied':
            cv = self._covars_
        else:
            cv = self._covars_[state]
        return sample_gaussian(self.means_[state], cv, self.covariance_type,
                               random_state=random_state)

    def _initialize_sufficient_statistics(self):
        stats = super(GaussianHMM, self)._initialize_sufficient_statistics()
        stats['post'] = np.zeros(self.n_components)
        stats['obs'] = np.zeros((self.n_components, self.n_features))
        stats['obs**2'] = np.zeros((self.n_components, self.n_features))
        if self.covariance_type in ('tied', 'full'):
            stats['obs*obs.T'] = np.zeros((self.n_components, self.n_features,
                                           self.n_features))
        # print(stats)
        return stats

    def _accumulate_sufficient_statistics(self, stats, obs, framelogprob, posteriors, fwdlattice, bwdlattice):
        super(GaussianHMM, self)._accumulate_sufficient_statistics(
            stats, obs, framelogprob, posteriors, fwdlattice, bwdlattice)

        if 'm' in self.params or 'c' in self.params:
            stats['post'] += posteriors.sum(axis=0)
            stats['obs'] += np.dot(posteriors.T, obs)

        if 'c' in self.params:
            if self.covariance_type in ('spherical', 'diag'):
                stats['obs**2'] += np.dot(posteriors.T, obs ** 2)
            elif self.covariance_type in ('tied', 'full'):
                # posteriors: (nt, nc); obs: (nt, nf); obs: (nt, nf)
                # -> (nc, nf, nf)
                stats['obs*obs.T'] += np.einsum(
                    'ij,ik,il->jkl', posteriors, obs, obs)

    def _do_mstep(self, stats):
        super(GaussianHMM, self)._do_mstep(stats)

        means_prior = self.means_prior
        means_weight = self.means_weight

        # TODO: find a proper reference for estimates for different
        #       covariance models.
        # Based on Huang, Acero, Hon, "Spoken Language Processing",
        # p. 443 - 445
        denom = stats['post'][:, np.newaxis]
        if 'm' in self.params:
            self.means_ = ((means_weight * means_prior + stats['obs'])
                           / (means_weight + denom))

        if 'c' in self.params:
            covars_prior = self.covars_prior
            covars_weight = self.covars_weight
            meandiff = self.means_ - means_prior

            if self.covariance_type in ('spherical', 'diag'):
                cv_num = (means_weight * meandiff**2
                          + stats['obs**2']
                          - 2 * self.means_ * stats['obs']
                          + self.means_**2 * denom)
                cv_den = max(covars_weight - 1, 0) + denom
                self._covars_ = \
                    (covars_prior + cv_num) / np.maximum(cv_den, 1e-5)
                if self.covariance_type == 'spherical':
                    self._covars_ = np.tile(
                        self._covars_.mean(1)[:, np.newaxis],
                        (1, self._covars_.shape[1]))
            elif self.covariance_type in ('tied', 'full'):
                cv_num = np.empty((self.n_components, self.n_features,
                                  self.n_features))
                for c in range(self.n_components):
                    obsmean = np.outer(stats['obs'][c], self.means_[c])

                    cv_num[c] = (means_weight * np.outer(meandiff[c],
                                                         meandiff[c])
                                 + stats['obs*obs.T'][c]
                                 - obsmean - obsmean.T
                                 + np.outer(self.means_[c], self.means_[c])
                                 * stats['post'][c])
                cvweight = max(covars_weight - self.n_features, 0)
                if self.covariance_type == 'tied':
                    self._covars_ = ((covars_prior + cv_num.sum(axis=0)) /
                                     (cvweight + stats['post'].sum()))
                elif self.covariance_type == 'full':
                    self._covars_ = ((covars_prior + cv_num) /
                                     (cvweight + stats['post'][:, None, None]))

#init values
#self.startprob_ -> pi_i (i:estado)
#self.transmat_ -> a_ij (i,j:estado)
#self.means_ -> mu_il (i:estado, l:feature/serie)
#self._covars_ -> sigma2_il (i:estado, l:feature/serie)
#self.rho_ -> rho_l (l:feature/serie)
#self.epsilon_ -> epsilon_l (l:feature/serie)
#self.tau_ -> tau_l (l:feature/serie)

#self.n_features = X.shape[1], con dimension L
#self.n_components -> number of states

#FIXME: revisar que cuando le asigné tau o std() a una variable esté bien eso, y que no tendría que haberle asignado _covars_ o algo distinto.
#FIXME: revisar que las cuentas finales con np.dot(np.sum()) y np.sum(np.sum()) hagan lo que tienen que hacer.


def gaussiana(x, mu, sigma2):
    return (1.0/np.sqrt(2*np.pi*sigma2))*np.exp(-0.5*((x - mu)**2)/sigma2)

class GaussianFSHMM(GaussianHMM):
    def __init__(self, k, **kwargs):
        super(GaussianFSHMM, self).__init__(**kwargs)
        self.k_factor_ = k

    def init_values_FS(self, X, epsilon=None, tau=None, rho=None): #default initialization is for the first example of the paper
        self.rho_ = np.ones(self.n_features)*0.5 if (rho is None) else rho

        self.epsilon_ = pd.DataFrame(X).mean().values if (epsilon is None) else epsilon

        self.tau_ = pd.DataFrame(X).std().values if (tau is None) else tau

    def select_hyperparams(self, X): #pass hyperparameters as arguments of this function
        self.p_ = np.ones(self.n_components)*2
        self.a_ = np.ones((self.n_components, self.n_components))*2

        self.b_ = pd.DataFrame(X).mean().values # self.b_ = X.mean().values

        self.m_ = self.means_.copy()
        for l in range(self.n_features):
            self.m_[0, l] = self.b_[l] - self.tau_[l]
            self.m_[1, l] = self.b_[l] + self.tau_[l]

        self.s_ = np.ones(self.m_.shape)*0.5
        self.zeta_ = np.ones(self.m_.shape)*0.5
        self.eta_ = np.ones(self.m_.shape)*0.5

        self.c_ = np.ones(self.n_features)*1.0
        self.nu_ = np.ones(self.n_features)*0.5
        self.psi_ = np.ones(self.n_features)*0.5
        self.k_ = np.ones(self.n_features)*self.k_factor_

    def compute_FS_ESTEP(self, X, gamma):
        I = self.n_components
        L = self.n_features
        T = X.shape[0]
        self.uilt = np.zeros((I, L, T))
        self.vilt = np.zeros((I, L, T))
        for i in range(I):
            for l in range(L):
                for t in range(T):
                    eilt = self.rho_[l]*gaussiana(X[t, l], self.means_[i, l], self._covars_[i, l])
                    hilt = (1.0-self.rho_[l])*gaussiana(X[t, l], self.epsilon_[l], self.tau_[l])
                    gilt = eilt + hilt
                    uilt = gamma[t, i]*eilt/gilt #tiene este orden gamma?
                    vilt = gamma[t, i] - uilt
                    self.uilt[i, l, t] = uilt
                    self.vilt[i, l, t] = vilt
        



    def fit(self, X, lengths=None):
        """Estimate model parameters.

        An initialization step is performed before entering the
        EM algorithm. If you want to avoid this step for a subset of
        the parameters, pass proper ``init_params`` keyword argument
        to estimator's constructor.

        Parameters
        ----------
        X : array-like, shape (n_samples, n_features)
            Feature matrix of individual samples.

        lengths : array-like of integers, shape (n_sequences, )
            Lengths of the individual sequences in ``X``. The sum of
            these should be ``n_samples``.

        Returns
        -------
        self : object
            Returns self.
        """
        X = check_array(X)
        self._init(X, lengths=lengths)
        self._check()
        self.init_values_FS(X)
        self.select_hyperparams(X)

        self.monitor_ = ConvergenceMonitor(self.tol, self.n_iter, self.verbose)
        for iter in range(self.n_iter):
            stats = self._initialize_sufficient_statistics()
            curr_logprob = 0
            for i, j in iter_from_X_lengths(X, lengths):
                framelogprob = self._compute_log_likelihood(X[i:j])
                logprob, fwdlattice = self._do_forward_pass(framelogprob)
                curr_logprob += logprob
                bwdlattice = self._do_backward_pass(framelogprob)
                posteriors = self._compute_posteriors(fwdlattice, bwdlattice) #posteriors <- gamma
                self._accumulate_sufficient_statistics(
                    stats, X[i:j], framelogprob, posteriors, fwdlattice,
                    bwdlattice)
                self.compute_FS_ESTEP(X, posteriors)

            # XXX must be before convergence check, because otherwise
            #     there won't be any updates for the case ``n_iter=1``.
            self._do_mstep(X, stats)

            self.monitor_.report(curr_logprob)
            if self.monitor_.converged:
                break

        return self

    def _do_mstep(self, X, stats):
        """Performs the M-step of EM algorithm.

        Parameters
        ----------
        stats : dict
            Sufficient statistics updated from all available samples.
        """
        # The ``np.where`` calls guard against updating forbidden states
        # or transitions in e.g. a left-right HMM.
        I = self.n_components
        L = self.n_features
        T = X.shape[0]

        if 's' in self.params:
            startprob_ = self.startprob_prior - 1.0 + stats['start'] #dimensions? #estos asumimos que estan bien
            self.startprob_ = np.where(self.startprob_ == 0.0,
                                       self.startprob_, startprob_)
            normalize(self.startprob_)
        if 't' in self.params:
            transmat_ = self.transmat_prior - 1.0 + stats['trans'] #dimensions? #estos asumimos que estan bien
            self.transmat_ = np.where(self.transmat_ == 0.0,
                                      self.transmat_, transmat_)
            normalize(self.transmat_, axis=1)

        means_prior = self.means_prior
        means_weight = self.means_weight

        # TODO: find a proper reference for estimates for different
        #       covariance models.
        # Based on Huang, Acero, Hon, "Spoken Language Processing",
        # p. 443 - 445
        denom = stats['post'][:, np.newaxis]
        if 'm' in self.params:
            for i in range(I):
                for l in range(L):
                    sil2 = self.s_[i,l]**2
                    sigmail2 = self._covars_[i,l]
                    term1 = sil2*np.dot(self.uilt[i, l, :], X[:,l])
                    num = term1 + self.m_[i, l]*sigmail2
                    den = sil2*np.sum(self.uilt[i, l, :]) + sigmail2
                    self.means_[i, l] = num/den
            #self.means_ = ((means_weight * means_prior + stats['obs'])
            #              / (means_weight + denom))

        if 'c' in self.params:
            covars_prior = self.covars_prior
            covars_weight = self.covars_weight
            meandiff = self.means_ - means_prior

            if self.covariance_type in ('spherical', 'diag'):
                for i in range(I):
                    for l in range(L):
                        term1 = np.dot(self.uilt[i, l, :], (X[:,l] - self.means_[i, l])**2)
                        num = term1 + 2 * self.eta_[i, l]
                        den = np.sum(self.uilt[i, l, :]) + 2*(self.zeta_[i, l] + 1.0)
                        self._covars_[i, l] = num/den

        for l in range(L):
            cl2 = self.c_[l]**2
            bl = self.b_[l]
            taul2 = self.tau_[l]**2
            psil = self.psi_[l]
            nul = self.nu_[l]
            epsilonl = self.epsilon_[l]
            kl = self.k_[l]
            hatT = T + 1 + kl

            epsilonl_num = cl2*np.dot(np.sum(self.vilt[:, l, :], axis=0), X[:, l]) + taul2*bl #????
            epsilonl_den = cl2*np.sum(np.sum(self.vilt[:, l, :], axis=0), axis=-1) + taul2
            self.epsilon_[l] = epsilonl_num / epsilonl_den

            taul_num = np.dot(np.sum(self.vilt[:, l, :], axis=0), (X[:, l] - epsilonl)**2) + 2*psil #????
            taul_den = np.sum(np.sum(self.vilt[:, l, :], axis=0), axis=-1) + 2 * (nul + 1.0)
            self.tau_[l] = np.sqrt(taul_num / taul_den)

            self.rho_[l] = (hatT - np.sqrt(hatT**2 - 4*kl*np.sum(np.sum(self.uilt[:, l, :], axis=0), axis=-1))) / (2*kl)