docs: def generate_states(self, initial_condition=None, with_noise=True, stateseq=None):
          """
          Jointly sample the discrete and continuous states
          """
          from pybasicbayes.util.stats import sample_discrete
          # Generate from the prior and raise exception if unstable
          T, K, n = self.T, self.num_states, self.D_latent
          A = self.trans_matrix
  
          # Initialize discrete state sequence
          dss = -1 * np.ones(T, dtype=np.int32) if stateseq is None else stateseq.astype(np.int32)
          assert dss.shape == (T,)
          gss = np.empty((T,n), dtype='double')
  
          if initial_condition is None:
              if dss[0] == -1:
                  dss[0] = sample_discrete(self.pi_0)
              gss[0] = self.init_dynamics_distns[dss[0]].rvs()
          else:
              dss[0] = initial_condition[0]
              gss[0] = initial_condition[1]
  
          for t in range(1,T):
              # Sample discrete state given previous continuous state
              if with_noise:
                  # Sample discre=te state from recurrent transition matrix
                  if dss[t] == -1:
                      dss[t] = sample_discrete(A[dss[t-1], :])
  
                  # Sample continuous state given current discrete state
                  gss[t] = self.dynamics_distns[dss[t-1]].\
                      rvs(x=np.hstack((gss[t-1][None,:], self.inputs[t-1][None,:])),
                          return_xy=False)
              else:
                  # Pick the most likely next discrete state and continuous state
                  if dss[t] == -1:
                      dss[t] = np.argmax(A[dss[t-1], :])
  
                  gss[t] = self.dynamics_distns[dss[t-1]]. \
                      predict(np.hstack((gss[t-1][None,:], self.inputs[t-1][None,:])))
              assert np.all(np.isfinite(gss[t])), "SLDS appears to be unstable!"
  
          self.stateseq = dss
          self.gaussian_states = gss