Included random_state parameter for reproducibility

pyhhmm/base.py

@@ -71,6 +71,7 @@ class BaseHMM(object):
     :param verbose: flag to be set to True if per-iteration convergence reports should 
         be printed. Defaults to True.
     :type verbose: bool, optional 
+    :param random_state: seed for the random number generator
     """
 
     def __init__(
@@ -83,6 +84,7 @@ def __init__(
         A_prior=1.0,
         learning_rate=0.,
         verbose=True,
+        random_state=None,
     ):
         """Constructor method."""
 
@@ -97,6 +99,8 @@ def __init__(
         self.A_prior = A_prior
         self.learning_rate = learning_rate
         self.verbose = verbose
+        self.random_state = random_state
+        self.rng = np.random.default_rng(random_state)
 
     def __str__(self):
         """Function to allow directly printing the object."""
@@ -329,13 +333,13 @@ def sample(self, n_sequences=1, n_samples=1, return_states=False):
         transmat_cdf = np.cumsum(self.A, axis=1)
 
         for _ in range(n_sequences):
-            currstate = (startprob_cdf > np.random.rand()).argmax()
+            currstate = (startprob_cdf > self.rng.rand()).argmax()
             state_sequence = [currstate]
             X = [self._generate_sample_from_state(currstate)]
 
             for _ in range(n_samples - 1):
                 currstate = (transmat_cdf[currstate]
-                             > np.random.rand()).argmax()
+                             > self.rng.rand()).argmax()
                 state_sequence.append(currstate)
                 X.append(self._generate_sample_from_state(currstate))
             samples.append(np.vstack(X))
@@ -373,12 +377,12 @@ def _init_model_params(self):
                 self.A = np.full((self.n_states, self.n_states), init)
         else:
             if 's' in self.init_params:
-                self.pi = np.random.dirichlet(
+                self.pi = self.rng.dirichlet(
                     alpha=self.pi_prior * np.ones(self.n_states), size=1
                 )[0]
 
             if 't' in self.init_params:
-                self.A = np.random.dirichlet(
+                self.A = self.rng.dirichlet(
                     alpha=self.A_prior * np.ones(self.n_states), size=self.n_states
                 )
 

pyhhmm/gaussian.py

@@ -66,6 +66,7 @@ class GaussianHMM(BaseHMM):
     :type learning_rate: float, optional
     :param verbose: flag to be set to True if per-iteration convergence reports should be printed. Defaults to True.
     :type verbose: bool, optional 
+    :param random_state: seed for the random number generator
     """
 
     def __init__(
@@ -84,6 +85,7 @@ def __init__(
         min_covar=1e-3,
         learning_rate=0.,
         verbose=False,
+        random_state=None,
     ):
         if covariance_type not in COVARIANCE_TYPES:
             raise ValueError(
@@ -99,6 +101,7 @@ def __init__(
             A_prior=A_prior,
             learning_rate=learning_rate,
             verbose=verbose,
+            random_state=random_state,
         )
 
         self.n_emissions = n_emissions
@@ -168,7 +171,7 @@ def _init_model_params(self, X):
         X_concat = concatenate_observation_sequences(X)
 
         if 'm' in self.init_params:
-            kmeans = cluster.KMeans(n_clusters=self.n_states)
+            kmeans = cluster.KMeans(n_clusters=self.n_states, random_state=self.random_state)
             kmeans.fit(X_concat)
             self.means = kmeans.cluster_centers_
         if 'c' in self.init_params:
@@ -481,4 +484,4 @@ def _generate_sample_from_state(self, state):
             from the emission distribution corresponding to a given state
         :rtype: array_like
         """
-        return np.random.multivariate_normal(self.means[state], self.covars[state])
+        return self.rng.multivariate_normal(self.means[state], self.covars[state])

pyhhmm/heterogeneous.py

@@ -76,6 +76,7 @@ class HeterogeneousHMM(BaseHMM):
     :type learning_rate: float, optional
     :param verbose: flag to be set to True if per-iteration convergence reports should be printed, defaults to True
     :type verbose: bool, optional
+    :param random_state: seed for the random number generator
     """
 
     def __init__(
@@ -98,6 +99,7 @@ def __init__(
         min_covar=1e-3,
         learning_rate=0,
         verbose=False,
+        random_state=None,
     ):
         """Constructor method.
 
@@ -125,6 +127,7 @@ def __init__(
             A_prior=A_prior,
             learning_rate=learning_rate,
             verbose=verbose,
+            random_state=random_state,
         )
 
         self.n_g_emissions = n_g_emissions
@@ -203,7 +206,7 @@ def _init_model_params(self, X):
             X, gidx=self.n_g_emissions)
 
         if 'm' in self.init_params:
-            kmeans = cluster.KMeans(n_clusters=self.n_states, random_state=0)
+            kmeans = cluster.KMeans(n_clusters=self.n_states, random_state=self.random_state)
             kmeans.fit(X_concat)
             self.means = kmeans.cluster_centers_
         if 'c' in self.init_params:
@@ -629,20 +632,20 @@ def _pdf(self, x, mean, covar):
         return multivariate_normal.pdf(x, mean=mean, cov=covar, allow_singular=True)
 
     def _generate_sample_from_state(self, state):
-        """ Generates a random sample from a given component.
+        """ Generates a random sample from fa given component.
         :param state: index of the component to condition on
         :type state: int
         :return: array of shape (n_g_features+n_d_features, ) containing a random sample
             from the emission distribution corresponding to a given state
         :rtype: array_like
         """
-        gauss_sample = np.random.multivariate_normal(
+        gauss_sample = self.rng.multivariate_normal(
             self.means[state], self.covars[state]
         )
 
         cat_sample = []
         for e in range(self.n_d_emissions):
             cdf = np.cumsum(self.B[e][state, :])
-            cat_sample.append((cdf > np.random.rand()).argmax())
+            cat_sample.append((cdf > self.rng.rand()).argmax())
 
         return np.concatenate([gauss_sample, cat_sample])

pyhhmm/multinomial.py

@@ -54,6 +54,7 @@ class MultinomialHMM(BaseHMM):
     :type learning_rate: float, optional
     :param verbose: flag to be set to True if per-iteration convergence reports should be printed, defaults to True
     :type verbose: bool, optional
+    :param random_state: seed for the random number generator
     """
 
     def __init__(
@@ -71,6 +72,7 @@ def __init__(
         state_no_train_de=None,
         learning_rate=0.1,
         verbose=True,
+        random_state=None,
     ):
         """Constructor method
 
@@ -91,6 +93,7 @@ def __init__(
             A_prior=A_prior,
             verbose=verbose,
             learning_rate=learning_rate,
+            random_state=random_state,
         )
         self.n_emissions = n_emissions
         self.n_features = n_features
@@ -157,7 +160,7 @@ def _init_model_params(self):
             else:
                 if self.nr_no_train_de == 0:
                     self.B = [
-                        np.random.rand(self.n_states, self.n_features[i])
+                        self.rng.rand(self.n_states, self.n_features[i])
                         for i in range(self.n_emissions)
                     ]
                     for i in range(self.n_emissions):
@@ -301,5 +304,5 @@ def _generate_sample_from_state(self, state):
         res = []
         for e in range(self.n_emissions):
             cdf = np.cumsum(self.B[e][state, :])
-            res.append((cdf > np.random.rand()).argmax())
+            res.append((cdf > self.rng.rand()).argmax())
         return np.asarray(res)