def max_likelihood(self, data, weights=None, stats=None, lmbda=0.1):
        """
        As an alternative to MCMC with Polya-gamma augmentation,
        we also implement maximum likelihood learning via gradient
        descent with autograd. This follows the pybasicbayes
        convention.

        :param data: list of tuples, (x,y), for each dataset.
        :param weights: Not used in this implementation.
        :param stats: Not used in this implementation.
        """
        import autograd.numpy as anp
        from autograd import value_and_grad, hessian_vector_product
        from scipy.optimize import minimize

        assert weights is None
        assert stats is None
        if not isinstance(data, list):
            assert isinstance(data, tuple) and len(data) == 2
            data = [data]

        # Define a helper function for the log of the logistic fn
        def loglogistic(psi):
            return psi - anp.log(1+anp.exp(psi))

        # optimize each row of A and b
        for n in range(self.D_out):

            # Define an objective function for the n-th row of hstack((A, b))
            # This is the negative log likelihood of the n-th column of data.
            def nll(abn):
                an, bn = abn[:-1], abn[-1]
                T = 0
                ll = 0
                for (x, y) in data:
                    T += x.shape[0]
                    yn = y[:, n]
                    psi = anp.dot(x, an) + bn
                    ll += anp.sum(yn * loglogistic(psi))
                    ll += anp.sum((1 - yn) * loglogistic(-1. * psi))

                # Include a penalty on the weights
                ll -= lmbda * T * anp.sum(an**2)
                ll -= lmbda * T * bn**2

                return -1 * ll / T

            abn0 = np.concatenate((self.A[n], self.b[n]))
            res = minimize(value_and_grad(nll), abn0,
                           tol=1e-3,
                           method="Newton-CG",
                           jac=True,
                           hessp=hessian_vector_product(nll))

            assert res.success
            self.A[n] = res.x[:-1]
            self.b[n] = res.x[-1]
As an alternative to MCMC with Polya-gamma augmentation,
        we also implement maximum likelihood learning via gradient
        descent with autograd. This follows the pybasicbayes
        convention.

        :param data: list of tuples, (x,y), for each dataset.
        :param weights: Not used in this implementation.
        :param stats: Not used in this implementation.