def simulate_indices(self, ts_length, init=None, num_reps=None,
                         random_state=None):
        """
        Simulate time series of state transitions, where state indices
        are returned.

        Parameters
        ----------
        ts_length : scalar(int)
            Length of each simulation.

        init : int or array_like(int, ndim=1), optional
            Initial state(s). If None, the initial state is randomly
            drawn.

        num_reps : scalar(int), optional(default=None)
            Number of repetitions of simulation.

        random_state : int or np.random.RandomState, optional
            Random seed (integer) or np.random.RandomState instance to
            set the initial state of the random number generator for
            reproducibility. If None, a randomly initialized RandomState
            is used.

        Returns
        -------
        X : ndarray(ndim=1 or 2)
            Array containing the state values of the sample path(s). See
            the `simulate` method for more information.

        """
        random_state = check_random_state(random_state)
        dim = 1  # Dimension of the returned array: 1 or 2

        msg_out_of_range = 'index {init} is out of the state space'

        try:
            k = len(init)  # init is an array
            dim = 2
            init_states = np.asarray(init, dtype=int)
            # Check init_states are in the state space
            if (init_states >= self.n).any() or (init_states < -self.n).any():
                idx = np.where(
                    (init_states >= self.n) + (init_states < -self.n)
                )[0][0]
                raise ValueError(msg_out_of_range.format(init=idx))
            if num_reps is not None:
                k *= num_reps
                init_states = np.tile(init_states, num_reps)
        except TypeError:  # init is a scalar(int) or None
            k = 1
            if num_reps is not None:
                dim = 2
                k = num_reps
            if init is None:
                init_states = random_state.randint(self.n, size=k)
            elif isinstance(init, numbers.Integral):
                # Check init is in the state space
                if init >= self.n or init < -self.n:
                    raise ValueError(msg_out_of_range.format(init=init))
                init_states = np.ones(k, dtype=int) * init
            else:
                raise ValueError(
                    'init must be int, array_like of ints, or None'
                )

        # === set up array to store output === #
        X = np.empty((k, ts_length), dtype=int)

        # Random values, uniformly sampled from [0, 1)
        random_values = random_state.random_sample(size=(k, ts_length-1))

        # Generate sample paths and store in X
        if not self.is_sparse:  # Dense
            _generate_sample_paths(
                self.cdfs, init_states, random_values, out=X
            )
        else:  # Sparse
            _generate_sample_paths_sparse(
                self.cdfs1d, self.P.indices, self.P.indptr, init_states,
                random_values, out=X
            )

        if dim == 1:
            return X[0]
        else:
            return X