def sample_bitstrings(self, n_samples):
        """
        Sample bitstrings from the distribution defined by the wavefunction.

        :param n_samples: The number of bitstrings to sample
        :return: An array of shape (n_samples, n_qubits)
        """
        possible_bitstrings = np.array(list(itertools.product((0, 1), repeat=len(self))))
        inds = np.random.choice(2 ** len(self), n_samples, p=self.probabilities())
        bitstrings = possible_bitstrings[inds, :]
        return bitstrings