r"""
        Learn the parameters of a Bayesian generalized linear model (GLM).

        Parameters
        ----------
        X : ndarray
            (N, d) array input dataset (N samples, d dimensions).
        y : ndarray
            (N,) array targets (N samples)
        likelihood : Object
            A likelihood object, see the likelihoods module.
        likelihood_args : sequence, optional
            sequence of arguments to pass to the likelihood function. These are
            non-learnable parameters. They can be scalars or arrays of length
            N.

def fit(self, X, y, likelihood_args=()):
        r"""
        Learn the parameters of a Bayesian generalized linear model (GLM).

        Parameters
        ----------
        X : ndarray
            (N, d) array input dataset (N samples, d dimensions).
        y : ndarray
            (N,) array targets (N samples)
        likelihood : Object
            A likelihood object, see the likelihoods module.
        likelihood_args : sequence, optional
            sequence of arguments to pass to the likelihood function. These are
            non-learnable parameters. They can be scalars or arrays of length
            N.
        """
        X, y = check_X_y(X, y)

        # Batch magnification factor
        N, _ = X.shape
        self.B_ = X.shape[0] / self.batch_size
        self.D_ = self.basis.get_dim(X)

        # Pack data
        likelihood_args = _reshape_likelihood_args(likelihood_args, N)
        data = (X, y) + likelihood_args

        # Pack params
        params = [Parameter(WGTRND, Bound(), shape=(self.D_, self.K)),
                  Parameter(COVRND, Positive(), shape=(self.D_, self.K)),
                  self.basis.regularizer,
                  self.likelihood.params,
                  self.basis.params]

        log.info("Optimising parameters...")
        self.__it = -self.nstarts  # Keeping track of iterations for logging
        nsgd = structured_sgd(logtrick_sgd(sgd))
        res = nsgd(self._elbo,
                   params,
                   data,
                   eval_obj=True,
                   maxiter=self.maxiter,
                   updater=self.updater,
                   batch_size=self.batch_size,
                   random_state=self.random_,
                   nstarts=self.nstarts
                   )

        # Unpack params
        (self.weights_,
         self.covariance_,
         self.regularizer_,
         self.like_hypers_,
         self.basis_hypers_
         ) = res.x

        log.info("Finished! reg = {}, likelihood_hypers = {}, "
                 "basis_hypers = {}, message: {}."
                 .format(self.regularizer_,
                         self.like_hypers_,
                         self.basis_hypers_,
                         res.message))

        return self