Generate samples from the posterior latent function mixtures of the GLM
        for query inputs, X*.

        Parameters
        ----------
        X : ndarray
            (N*, d) array query input dataset (N* samples, D dimensions).
        nsamples : int
            Number of samples for sampling the latent function.
        genaxis : int
            Axis to return samples from, i.e.
            - ``genaxis=1`` will give you one sample at a time of f for ALL
                observations (so it will iterate over nsamples).
            - ``genaxis=0`` will give you all samples of f for ONE
                observation at a time (so it will iterate through X*, row by
                row)

        Yields
        ------
        fsamples : ndarray
            of shape (N*,) if ``genaxis=1`` with each call being a sample
            from the mixture of latent functions over all N*. Or of shape
            (nsamples,) if ``genaxis=0``, with each call being a all
            samples for an observation, n in N*.

def _sample_func(self, X, nsamples, genaxis=1):
        """
        Generate samples from the posterior latent function mixtures of the GLM
        for query inputs, X*.

        Parameters
        ----------
        X : ndarray
            (N*, d) array query input dataset (N* samples, D dimensions).
        nsamples : int
            Number of samples for sampling the latent function.
        genaxis : int
            Axis to return samples from, i.e.
            - ``genaxis=1`` will give you one sample at a time of f for ALL
                observations (so it will iterate over nsamples).
            - ``genaxis=0`` will give you all samples of f for ONE
                observation at a time (so it will iterate through X*, row by
                row)

        Yields
        ------
        fsamples : ndarray
            of shape (N*,) if ``genaxis=1`` with each call being a sample
            from the mixture of latent functions over all N*. Or of shape
            (nsamples,) if ``genaxis=0``, with each call being a all
            samples for an observation, n in N*.
        """
        check_is_fitted(self, ['weights_', 'covariance_', 'basis_hypers_',
                               'like_hypers_', 'regularizer_'])
        X = check_array(X)
        D, K = self.weights_.shape

        # Generate weight samples from all mixture components
        k = self.random_.randint(0, K, size=(nsamples,))
        w = self.weights_[:, k] + self.random_.randn(D, nsamples) \
            * np.sqrt(self.covariance_[:, k])

        # Do this here for *massive* speed improvements
        Phi = self.basis.transform(X, *atleast_list(self.basis_hypers_))

        # Now generate latent functions samples either colwise or rowwise
        if genaxis == 1:
            fs = (Phi.dot(ws) for ws in w.T)
        elif genaxis == 0:
            fs = (phi_n.dot(w) for phi_n in Phi)
        else:
            raise ValueError("Invalid axis to generate samples from")

        return fs