Repository: glm-tools/pyglmnet

Issue:
Make bleeding edge version installation easier for non developers
According to @jasmainak 's suggestion; provide something like this and modify `doc/install.rst`:
```
$ pip install https://api.github.com/repos/glm-tools/pyglmnet/zipball/master
```

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/glm-tools/pyglmnet/pull/322?src=pr&el=h1) Report
> Merging [#322](https://codecov.io/gh/glm-tools/pyglmnet/pull/322?src=pr&el=desc) into [master](https://codecov.io/gh/glm-tools/pyglmnet/commit/83b3883167caf96116850b5499ceba50d2224818?src=pr&el=desc) will **not change** coverage.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/glm-tools/pyglmnet/pull/322/graphs/tree.svg?width=650&token=kkw3ydHRLW&height=150&src=pr)](https://codecov.io/gh/glm-tools/pyglmnet/pull/322?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master     #322   +/-   ##
=======================================
  Coverage   74.56%   74.56%           
=======================================
  Files           4        4           
  Lines         629      629           
  Branches      128      128           
=======================================
  Hits          469      469           
  Misses        125      125           
  Partials       35       35
```



------

[Continue to review full report at Codecov](https://codecov.io/gh/glm-tools/pyglmnet/pull/322?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/glm-tools/pyglmnet/pull/322?src=pr&el=footer). Last update [83b3883...7d89aa4](https://codecov.io/gh/glm-tools/pyglmnet/pull/322?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

pavanramkumar: @jasmainak LGTM!

Patch:
diff --git a/README.rst b/README.rst
index 1e4da9f..67e2beb 100644
--- a/README.rst
+++ b/README.rst
@@ -108,13 +108,7 @@ Clone the repository.
 
 .. code:: bash
 
-    $ git clone http://github.com/glm-tools/pyglmnet
-
-Install ``pyglmnet`` using ``setup.py`` as follows
-
-.. code:: bash
-
-    $ python setup.py develop
+    $ pip install https://api.github.com/repos/glm-tools/pyglmnet/zipball/master
 
 Getting Started
 ~~~~~~~~~~~~~~~
@@ -132,15 +126,15 @@ Here is an example on how to use the ``GLM`` estimator.
    distr = 'poisson'
 
    # sample a sparse model
-   beta0 = np.random.normal(0.0, 1.0, 1)
-   beta = sps.rand(n_features, 1, 0.1)
-   beta = np.array(beta.todense())
+   beta0 = np.random.rand()
+   beta = np.random.random(n_features)
+   beta[beta < 0.9] = 0
 
    # simulate data
    Xtrain = np.random.normal(0.0, 1.0, [n_samples, n_features])
-   ytrain = simulate_glm('poisson', beta0, beta, Xtrain)[:, 0]
+   ytrain = simulate_glm('poisson', beta0, beta, Xtrain)
    Xtest = np.random.normal(0.0, 1.0, [n_samples, n_features])
-   ytest = simulate_glm('poisson', beta0, beta, Xtest)[:, 0]
+   ytest = simulate_glm('poisson', beta0, beta, Xtest)
 
    # create an instance of the GLM class
    glm = GLM(distr='poisson', score_metric='deviance')
diff --git a/doc/install.rst b/doc/install.rst
index 54ecbbd..ed6a095 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -23,8 +23,7 @@ the dependencies using ``pip``.
 
 .. code-block:: bash
 
-   pip install numpy
-   pip install scipy
+   pip install numpy scipy
    pip install -U scikit-learn
 
 Get pyglmnet
@@ -39,11 +38,4 @@ To install the development version of ``pyglmnet``, first clone the repository
 
 .. code-block:: bash
 
-    git clone http://github.com/glm-tools/pyglmnet
-
-Then run ``setup.py`` to install the library
-
-.. code-block:: bash
-
-    cd pyglmnet
-    python setup.py develop
+	pip install https://api.github.com/repos/glm-tools/pyglmnet/zipball/master
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
index 8624804..d489e8f 100644
--- a/doc/whats_new.rst
+++ b/doc/whats_new.rst
@@ -20,6 +20,7 @@ Changelog
     - Improved testing of gradients for optimization by `Mainak Jas`_ and `Pavan Ramkumar`_.
     - Improved testing of estimated coefficients at convergence by `Mainak Jas`_ and `Pavan Ramkumar`_
     - Add `n_iter_` attribute to :class:`pyglmnet.GLM` by `Mainak Jas`_, `Olivier Pieters`_, `Peter Foley`_ and `Chris Rogers`_
+    - Add option to fit intercept or not with a `fit_intercept` boolean by `Drew Harris`_, `Mainak Jas`_ and `Pavan Ramkumar`_
 
 BUG
 ~~~
@@ -28,11 +29,12 @@ BUG
     - Fixed incorrect usage of the random_state parameter by `Giovanni De Toni`_.
     - Fixed incorrect proximal operator for group lasso by `Yu Umegaki`_.
     - Changed stopping criteria for convergence (a threshold on the change in objective value) which
-      stopped too learly. The new criteria is a threshold on the norm of the gradient, by `Pavan Ramkumar`_.
+      stopped too early. The new criteria is a threshold on the norm of the gradient, by `Pavan Ramkumar`_.
     - Fixed `group` parameter not being passed for Group Lasso by `Beibin Li`_.
     - Fix URL paths when running on Windows by `Scott Otterson`_.
     - Made temporary file handling OS independent and usable on a cluster by `Scott Otterson`_.
-    - Replace StratifiedKFold() with KFold(), as StratifiedKFold() doesn't work for continuous values by `Scott Otterson`_.
+    - Replace `StratifiedKFold()` with `KFold()`, as `StratifiedKFold()` doesn't work for continuous values by `Scott Otterson`_.
+    - :func:`pyglmnet.simulate_glm` accepts `beta0` as float and `beta` of shape (n_features,) only, by `Mainak Jas`_
 
 API
 ~~~
@@ -76,3 +78,4 @@ Changelog
 .. _Peter Foley: https://github.com/peterfoley
 .. _Olivier Pieters: https://olivierpieters.be/
 .. _Chris Rogers: https://github.com/cxrodgers
+.. _Drew Harris: https://github.com/DreHar
diff --git a/pyglmnet/pyglmnet.py b/pyglmnet/pyglmnet.py
index 497de60..ba66e46 100644
--- a/pyglmnet/pyglmnet.py
+++ b/pyglmnet/pyglmnet.py
@@ -6,7 +6,7 @@ from copy import deepcopy
 import numpy as np
 from scipy.special import expit
 from scipy.stats import norm
-from .utils import logger, set_log_level
+from .utils import logger, set_log_level, _check_params
 from .base import BaseEstimator, is_classifier, check_version
 
 
@@ -71,20 +71,29 @@ def _probit_g6(z, pdfz, cdfz, thresh=5):
     return res
 
 
-def _lmb(distr, beta0, beta, X, eta):
+def _z(beta0, beta, X, fit_intercept):
+    """Compute z to be passed through non-linearity"""
+    if fit_intercept:
+        z = beta0 + np.dot(X, beta)
+    else:
+        z = np.dot(X, np.r_[beta0, beta])
+    return z
+
+
+def _lmb(distr, beta0, beta, X, eta, fit_intercept=True):
     """Conditional intensity function."""
-    z = beta0 + np.dot(X, beta)
-    return _mu(distr, z, eta)
+    z = _z(beta0, beta, X, fit_intercept)
+    return _mu(distr, z, eta, fit_intercept)
 
 
-def _mu(distr, z, eta):
+def _mu(distr, z, eta, fit_intercept):
     """The non-linearity (inverse link)."""
     if distr in ['softplus', 'gamma']:
         mu = np.log1p(np.exp(z))
     elif distr == 'poisson':
         mu = z.copy()
-        intercept = (1 - eta) * np.exp(eta)
-        mu[z > eta] = z[z > eta] * np.exp(eta) + intercept
+        beta0 = (1 - eta) * np.exp(eta) if fit_intercept else 0.
+        mu[z > eta] = z[z > eta] * np.exp(eta) + beta0
         mu[z <= eta] = np.exp(z[z <= eta])
     elif distr == 'gaussian':
         mu = z
@@ -184,76 +193,99 @@ def _L1penalty(beta, group=None):
     return L1penalty
 
 
-def _loss(distr, alpha, Tau, reg_lambda, X, y, eta, group, beta):
+def _loss(distr, alpha, Tau, reg_lambda, X, y, eta, group, beta,
+          fit_intercept=True):
     """Define the objective function for elastic net."""
-    n_samples = X.shape[0]
-    z = beta[0] + np.dot(X, beta[1:])
-    y_hat = _mu(distr, z, eta)
+    n_samples, n_features = X.shape
+    z = _z(beta[0], beta[1:], X, fit_intercept)
+    y_hat = _mu(distr, z, eta, fit_intercept)
     L = 1. / n_samples * _logL(distr, y, y_hat, z)
-    P = _penalty(alpha, beta[1:], Tau, group)
+    if fit_intercept:
+        P = _penalty(alpha, beta[1:], Tau, group)
+    else:
+        P = _penalty(alpha, beta, Tau, group)
     J = -L + reg_lambda * P
     return J
 
 
-def _L2loss(distr, alpha, Tau, reg_lambda, X, y, eta, group, beta):
+def _L2loss(distr, alpha, Tau, reg_lambda, X, y, eta, group, beta,
+            fit_intercept=True):
     """Define the objective function for elastic net."""
-    n_samples = X.shape[0]
-    z = beta[0] + np.dot(X, beta[1:])
-    y_hat = _mu(distr, z, eta)
+    n_samples, n_features = X.shape
+    z = _z(beta[0], beta[1:], X, fit_intercept)
+    y_hat = _mu(distr, z, eta, fit_intercept)
     L = 1. / n_samples * _logL(distr, y, y_hat, z)
-    P = 0.5 * (1 - alpha) * _L2penalty(beta[1:], Tau)
+    if fit_intercept:
+        P = 0.5 * (1 - alpha) * _L2penalty(beta[1:], Tau)
+    else:
+        P = 0.5 * (1 - alpha) * _L2penalty(beta, Tau)
     J = -L + reg_lambda * P
     return J
 
 
-def _grad_L2loss(distr, alpha, Tau, reg_lambda, X, y, eta, beta):
+def _grad_L2loss(distr, alpha, Tau, reg_lambda, X, y, eta, beta,
+                 fit_intercept=True):
     """The gradient."""
-    n_samples = np.float(X.shape[0])
+    n_samples, n_features = X.shape
+    n_samples = np.float(n_samples)
 
     if Tau is None:
-        Tau = np.eye(beta[1:].shape[0])
+        if fit_intercept:
+            Tau = np.eye(beta[1:].shape[0])
+        else:
+            Tau = np.eye(beta.shape[0])
     InvCov = np.dot(Tau.T, Tau)
 
-    z = beta[0] + np.dot(X, beta[1:])
-    mu = _mu(distr, z, eta)
+    z = _z(beta[0], beta[1:], X, fit_intercept)
+    mu = _mu(distr, z, eta, fit_intercept)
     grad_mu = _grad_mu(distr, z, eta)
 
+    grad_beta0 = 0.
     if distr in ['poisson', 'softplus']:
-        grad_beta0 = np.sum(grad_mu) - np.sum(y * grad_mu / mu)
+        if fit_intercept:
+            grad_beta0 = np.sum(grad_mu) - np.sum(y * grad_mu / mu)
         grad_beta = ((np.dot(grad_mu.T, X) -
                       np.dot((y * grad_mu / mu).T, X)).T)
 
     elif distr == 'gaussian':
-        grad_beta0 = np.sum((mu - y) * grad_mu)
+        if fit_intercept:
+            grad_beta0 = np.sum((mu - y) * grad_mu)
         grad_beta = np.dot((mu - y).T, X * grad_mu[:, None]).T
 
     elif distr == 'binomial':
-        grad_beta0 = np.sum(mu - y)
+        if fit_intercept:
+            grad_beta0 = np.sum(mu - y)
         grad_beta = np.dot((mu - y).T, X).T
 
     elif distr == 'probit':
         grad_logl = (y * _probit_g3(z, grad_mu, mu) -
                      (1 - y) * _probit_g4(z, grad_mu, mu))
-        grad_beta0 = -np.sum(grad_logl)
+        if fit_intercept:
+            grad_beta0 = -np.sum(grad_logl)
         grad_beta = -np.dot(grad_logl.T, X).T
 
     elif distr == 'gamma':
         nu = 1.
         grad_logl = (y / mu ** 2 - 1 / mu) * grad_mu
-        grad_beta0 = -nu * np.sum(grad_logl)
+        if fit_intercept:
+            grad_beta0 = -nu * np.sum(grad_logl)
         grad_beta = -nu * np.dot(grad_logl.T, X).T
 
     grad_beta0 *= 1. / n_samples
     grad_beta *= 1. / n_samples
-    grad_beta += reg_lambda * (1 - alpha) * np.dot(InvCov, beta[1:])
-    n_features = X.shape[1]
-    g = np.zeros((n_features + 1, ))
-    g[0] = grad_beta0
-    g[1:] = grad_beta
+    if fit_intercept:
+        grad_beta += reg_lambda * (1 - alpha) * np.dot(InvCov, beta[1:])
+        g = np.zeros((n_features + 1, ))
+        g[0] = grad_beta0
+        g[1:] = grad_beta
+    else:
+        grad_beta += reg_lambda * (1 - alpha) * np.dot(InvCov, beta)
+        g = grad_beta
+
     return g
 
 
-def _gradhess_logloss_1d(distr, xk, y, z, eta):
+def _gradhess_logloss_1d(distr, xk, y, z, eta, fit_intercept=True):
     """
     Compute gradient (1st derivative)
     and Hessian (2nd derivative)
@@ -278,7 +310,7 @@ def _gradhess_logloss_1d(distr, xk, y, z, eta):
     n_samples = xk.shape[0]
 
     if distr == 'softplus':
-        mu = _mu(distr, z, eta)
+        mu = _mu(distr, z, eta, fit_intercept)
         s = expit(z)
         gk = np.sum(s * xk) - np.sum(y * s / mu * xk)
 
@@ -287,7 +319,7 @@ def _gradhess_logloss_1d(distr, xk, y, z, eta):
         hk = np.sum(grad_s * xk ** 2) - np.sum(y * grad_s_by_mu * xk ** 2)
 
     elif distr == 'poisson':
-        mu = _mu(distr, z, eta)
+        mu = _mu(distr, z, eta, fit_intercept)
         s = expit(z)
         gk = np.sum((mu[z <= eta] - y[z <= eta]) *
                     xk[z <= eta]) + \
@@ -304,7 +336,7 @@ def _gradhess_logloss_1d(distr, xk, y, z, eta):
         hk = np.sum(xk * xk)
 
     elif distr == 'binomial':
-        mu = _mu(distr, z, eta)
+        mu = _mu(distr, z, eta, fit_intercept)
         gk = np.sum((mu - y) * xk)
         hk = np.sum(mu * (1.0 - mu) * xk * xk)
 
@@ -354,6 +386,12 @@ def simulate_glm(distr, beta0, beta, X, eta=2.0, random_state=None,
         raise ValueError("'distr' must be in %s, got %s"
                          % (repr(ALLOWED_DISTRS), distr))
 
+    if not isinstance(beta0, float):
+        raise ValueError("'beta0' must be float, got %s" % type(beta0))
+
+    if beta.ndim != 1:
+        raise ValueError("'beta' must be 1D, got %dD" % beta.ndim)
+
     if not sample:
         return _lmb(distr, beta0, beta, X, eta)
 
@@ -447,7 +485,11 @@ class GLM(BaseEstimator):
         specifies the scoring metric.
         one of either 'deviance' or 'pseudo_R2'.
         default: 'deviance'
-    random_state: int
+    fit_intercept: boolean
+        specifies if a constant (a.k.a. bias or intercept) should be
+        added to the decision function.
+        default: True
+    random_state : int
         seed of the random number generator used to initialize the solution.
         default: 0
     verbose: boolean or int
@@ -492,14 +534,11 @@ class GLM(BaseEstimator):
                  solver='batch-gradient',
                  learning_rate=2e-1, max_iter=1000,
                  tol=1e-3, eta=2.0, score_metric='deviance',
+                 fit_intercept=True,
                  random_state=0, callback=None, verbose=False):
 
-        if not isinstance(max_iter, int):
-            raise ValueError('max_iter must be of type int')
-
-        if distr not in ALLOWED_DISTRS:
-            raise ValueError('distr must be one of %s. Got:'
-                             '%s' % (', '.join(ALLOWED_DISTRS), distr))
+        _check_params(distr=distr, max_iter=max_iter,
+                      fit_intercept=fit_intercept)
 
         self.distr = distr
         self.alpha = alpha
@@ -516,6 +555,7 @@ class GLM(BaseEstimator):
         self.tol = tol
         self.eta = eta
         self.score_metric = score_metric
+        self.fit_intercept = fit_intercept
         self.random_state = random_state
         self.rng = np.random.RandomState(self.random_state)
         self.callback = callback
@@ -627,7 +667,7 @@ class GLM(BaseEstimator):
 
             return result
 
-    def _cdfast(self, X, y, ActiveSet, beta, rl):
+    def _cdfast(self, X, y, ActiveSet, beta, rl, fit_intercept=True):
         """
         Perform one cycle of Newton updates for all coordinates.
 
@@ -640,10 +680,10 @@ class GLM(BaseEstimator):
             Labels to the data
             n_samples x 1
         ActiveSet: array
-            n_features + 1 x 1
+            n_features + 1 x 1, or n_features
             Active set storing which betas are non-zero
         beta: array
-            n_features + 1 x 1
+            n_features + 1 x 1, or n_features
             Parameters to be updated
         rl: float
             Regularization lambda
@@ -651,35 +691,44 @@ class GLM(BaseEstimator):
         Returns
         -------
         beta: array
-            (n_features + 1) x 1
+            (n_features + 1) x 1, or (n_features)
             Updated parameters
         """
-        n_samples = X.shape[0]
-        n_features = X.shape[1]
+        n_samples, n_features = X.shape
         reg_scale = rl * (1 - self.alpha)
-        z = beta[0] + np.dot(X, beta[1:])
+        z = _z(beta[0], beta[1:], X, fit_intercept)
 
-        for k in range(0, n_features + 1):
+        for k in range(0, n_features + int(fit_intercept)):
             # Only update parameters in active set
             if ActiveSet[k] != 0:
-                if k > 0:
-                    xk = X[:, k - 1]
+                if fit_intercept:
+                    if k == 0:
+                        xk = np.ones((n_samples, ))
+                    else:
+                        xk = X[:, k - 1]
                 else:
-                    xk = np.ones((n_samples, ))
+                    xk = X[:, k]
 
                 # Calculate grad and hess of log likelihood term
-                gk, hk = _gradhess_logloss_1d(self.distr, xk, y, z, self.eta)
+                gk, hk = _gradhess_logloss_1d(self.distr, xk, y, z, self.eta,
+                                              fit_intercept)
 
                 # Add grad and hess of regularization term
                 if self.Tau is None:
-                    gk_reg = beta[k]
-                    hk_reg = 1.0
+                    if k == 0 and fit_intercept:
+                        gk_reg, hk_reg = 0.0, 0.0
+                    else:
+                        gk_reg, hk_reg = beta[k], 1.0
                 else:
                     InvCov = np.dot(self.Tau.T, self.Tau)
-                    gk_reg = np.sum(InvCov[k - 1, :] * beta[1:])
-                    hk_reg = InvCov[k - 1, k - 1]
-                gk += np.ravel([reg_scale * gk_reg if k > 0 else 0.0])
-                hk += np.ravel([reg_scale * hk_reg if k > 0 else 0.0])
+                    if fit_intercept:
+                        gk_reg = np.sum(InvCov[k - 1, :] * beta[1:])
+                        hk_reg = InvCov[k - 1, k - 1]
+                    else:
+                        gk_reg = np.sum(InvCov[k, :] * beta)
+                        hk_reg = InvCov[k, k]
+                gk += reg_scale * gk_reg
+                hk += reg_scale * hk_reg
 
                 # Update parameters, z
                 update = 1. / hk * gk
@@ -702,7 +751,6 @@ class GLM(BaseEstimator):
         self: instance of GLM
             The fitted model.
         """
-
         # checks for group
         if self.group is not None:
             self.group = np.array(self.group)
@@ -720,9 +768,11 @@ class GLM(BaseEstimator):
                    .format(type(X), type(y)))
             raise ValueError(msg)
 
-        if X.ndim != 2 or y.ndim != 1:
-            msg = "X must be a 2D ndarray, and y must be 1D"
-            raise ValueError(msg)
+        if X.ndim != 2:
+            raise ValueError("X must be a 2D array, got %sD" % X.ndim)
+
+        if y.ndim != 1:
+            raise ValueError("y must be 1D, got %sD" % y.ndim)
 
         n_observations, n_features = X.shape
 
@@ -732,15 +782,22 @@ class GLM(BaseEstimator):
                              .format(n_observations, len(y)))
 
         # Initialize parameters
-        beta = np.zeros((n_features + 1,))
-        if self.beta0_ is None and self.beta_ is None:
-            beta[0] = 1 / (n_features + 1) * \
-                self.rng.normal(0.0, 1.0, 1)
-            beta[1:] = 1 / (n_features + 1) * \
-                self.rng.normal(0.0, 1.0, (n_features, ))
+        beta = np.zeros((n_features + int(self.fit_intercept),))
+        if self.fit_intercept:
+            if self.beta0_ is None and self.beta_ is None:
+                beta[0] = 1 / (n_features + 1) * \
+                    self.rng.normal(0.0, 1.0, 1)
+                beta[1:] = 1 / (n_features + 1) * \
+                    self.rng.normal(0.0, 1.0, (n_features, ))
+            else:
+                beta[0] = self.beta0_
+                beta[1:] = self.beta_
         else:
-            beta[0] = self.beta0_
-            beta[1:] = self.beta_
+            if self.beta0_ is None and self.beta_ is None:
+                beta = 1 / (n_features + 1) * \
+                    self.rng.normal(0.0, 1.0, (n_features, ))
+            else:
+                beta = self.beta_
 
         logger.info('Lambda: %6.4f' % self.reg_lambda)
 
@@ -749,7 +806,8 @@ class GLM(BaseEstimator):
         reg_lambda = self.reg_lambda
 
         if self.solver == 'cdfast':
-            ActiveSet = np.ones(n_features + 1)     # init active set
+            # init active set
+            ActiveSet = np.ones_like(beta)
 
         # Iterative updates
         for t in range(0, self.max_iter):
@@ -759,7 +817,7 @@ class GLM(BaseEstimator):
                 grad = _grad_L2loss(self.distr,
                                     alpha, self.Tau,
                                     reg_lambda, X, y, self.eta,
-                                    beta)
+                                    beta, self.fit_intercept)
                 # Converged if the norm(gradient) < tol
                 if (t > 1) and (np.linalg.norm(grad) < tol):
                     msg = ('\tConverged in {0:d} iterations'.format(t))
@@ -770,7 +828,8 @@ class GLM(BaseEstimator):
             elif self.solver == 'cdfast':
                 beta_old = deepcopy(beta)
                 beta = \
-                    self._cdfast(X, y, ActiveSet, beta, reg_lambda)
+                    self._cdfast(X, y, ActiveSet, beta, reg_lambda,
+                                 self.fit_intercept)
                 # Converged if the norm(update) < tol
                 if (t > 1) and (np.linalg.norm(beta - beta_old) < tol):
                     msg = ('\tConverged in {0:d} iterations'.format(t))
@@ -781,12 +840,16 @@ class GLM(BaseEstimator):
                                  "'('batch-gradient', 'cdfast'), got %s."
                                  % (self.solver))
             # Apply proximal operator
-            beta[1:] = self._prox(beta[1:], reg_lambda * alpha)
+            if self.fit_intercept:
+                beta[1:] = self._prox(beta[1:], reg_lambda * alpha)
+            else:
+                beta = self._prox(beta, reg_lambda * alpha)
 
             # Update active set
             if self.solver == 'cdfast':
                 ActiveSet[beta == 0] = 0
-                ActiveSet[0] = 1.
+                if self.fit_intercept:
+                    ActiveSet[0] = 1.
 
             # Compute and save loss if callbacks are requested
             if callable(self.callback):
@@ -797,8 +860,12 @@ class GLM(BaseEstimator):
                 "Reached max number of iterations without convergence.")
 
         # Update the estimated variables
-        self.beta0_ = beta[0]
-        self.beta_ = beta[1:]
+        if self.fit_intercept:
+            self.beta0_ = beta[0]
+            self.beta_ = beta[1:]
+        else:
+            self.beta0_ = 0
+            self.beta_ = beta
         self.ynull_ = np.mean(y)
         return self
 
@@ -819,8 +886,8 @@ class GLM(BaseEstimator):
             raise ValueError('Input data should be of type ndarray (got %s).'
                              % type(X))
 
-        yhat = _lmb(self.distr, self.beta0_,
-                    self.beta_, X, self.eta)
+        yhat = _lmb(self.distr, self.beta0_, self.beta_, X, self.eta,
+                    fit_intercept=True)
 
         if self.distr == 'binomial':
             yhat = (yhat > 0.5).astype(int)
@@ -854,8 +921,8 @@ class GLM(BaseEstimator):
             raise ValueError('Input data should be of type ndarray (got %s).'
                              % type(X))
 
-        yhat = _lmb(self.distr,
-                    self.beta0_, self.beta_, X, self.eta)
+        yhat = _lmb(self.distr, self.beta0_, self.beta_, X, self.eta,
+                    fit_intercept=True)
         yhat = np.asarray(yhat)
         return yhat
 
@@ -1008,7 +1075,11 @@ class GLMCV(object):
         specifies the scoring metric.
         one of either 'deviance' or 'pseudo_R2'.
         default: 'deviance'
-    random_state: int
+    fit_intercept: boolean
+        specifies if a constant (a.k.a. bias or intercept) should be
+        added to the decision function.
+        default: True
+    random_state : int
         seed of the random number generator used to initialize the solution.
         default: 0
     verbose: boolean or int
@@ -1047,6 +1118,7 @@ class GLMCV(object):
                  solver='batch-gradient',
                  learning_rate=2e-1, max_iter=1000,
                  tol=1e-3, eta=2.0, score_metric='deviance',
+                 fit_intercept=True,
                  random_state=0, verbose=False):
 
         if reg_lambda is None:
@@ -1055,12 +1127,8 @@ class GLMCV(object):
         if not isinstance(reg_lambda, (list, np.ndarray)):
             reg_lambda = [reg_lambda]
 
-        if distr not in ALLOWED_DISTRS:
-            raise ValueError('distr must be one of %s, Got '
-                             '%s' % (', '.join(ALLOWED_DISTRS), distr))
-
-        if not isinstance(max_iter, int):
-            raise ValueError('max_iter must be of type int')
+        _check_params(distr=distr, max_iter=max_iter,
+                      fit_intercept=fit_intercept)
 
         self.distr = distr
         self.alpha = alpha
@@ -1080,6 +1148,7 @@ class GLMCV(object):
         self.tol = tol
         self.eta = eta
         self.score_metric = score_metric
+        self.fit_intercept = fit_intercept
         self.random_state = random_state
         self.verbose = verbose
         set_log_level(verbose)
@@ -1150,6 +1219,7 @@ class GLMCV(object):
                       tol=self.tol,
                       eta=self.eta,
                       score_metric=self.score_metric,
+                      fit_intercept=self.fit_intercept,
                       random_state=self.random_state,
                       verbose=self.verbose)
             logger.info('Lambda: %6.4f' % rl)
@@ -1172,6 +1242,7 @@ class GLMCV(object):
                 glm.beta0_, glm.beta_ = self.beta0_, self.beta_
             else:
                 glm.beta0_, glm.beta_ = glms[-1].beta0_, glms[-1].beta_
+
             glm.fit(X, y)
             glms.append(glm)
         # Update the estimated variables
diff --git a/pyglmnet/utils.py b/pyglmnet/utils.py
index 14e618f..a46d2f0 100644
--- a/pyglmnet/utils.py
+++ b/pyglmnet/utils.py
@@ -90,6 +90,21 @@ def tikhonov_from_prior(prior_cov, n_samples, threshold=0.0001):
     return Tau
 
 
+def _check_params(distr, max_iter, fit_intercept):
+    from .pyglmnet import ALLOWED_DISTRS
+
+    if distr not in ALLOWED_DISTRS:
+        raise ValueError('distr must be one of %s, Got '
+                         '%s' % (', '.join(ALLOWED_DISTRS), distr))
+
+    if not isinstance(max_iter, int):
+        raise ValueError('max_iter must be of type int')
+
+    if not isinstance(fit_intercept, bool):
+        raise ValueError('fit_intercept must be bool, got %s'
+                         % type(fit_intercept))
+
+
 def set_log_level(verbose):
     """Convenience function for setting the log level.
 


Test patch:
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
index 6138938..daf4199 100644
--- a/tests/test_metrics.py
+++ b/tests/test_metrics.py
@@ -8,7 +8,7 @@ def test_deviance():
     """Test deviance."""
     n_samples, n_features = 1000, 100
 
-    beta0 = np.random.normal(0.0, 1.0, 1)
+    beta0 = np.random.rand()
     beta = np.random.normal(0.0, 1.0, n_features)
 
     # sample train and test data
@@ -26,7 +26,7 @@ def test_pseudoR2():
     """Test pseudo r2."""
     n_samples, n_features = 1000, 100
 
-    beta0 = np.random.normal(0.0, 1.0, 1)
+    beta0 = np.random.rand()
     beta = np.random.normal(0.0, 1.0, n_features)
 
     # sample train and test data
@@ -44,13 +44,15 @@ def test_accuracy():
     """Testing accuracy."""
     n_samples, n_features, n_classes = 1000, 100, 2
 
-    beta0 = np.random.normal(0.0, 1.0, 1)
-    beta = np.random.normal(0.0, 1.0, (n_features, n_classes))
+    beta0 = np.random.rand()
+    betas = np.random.normal(0.0, 1.0, (n_features, n_classes))
 
     # sample train and test data
     glm_sim = GLM(distr='binomial', score_metric='accuracy')
     X = np.random.randn(n_samples, n_features)
-    y = simulate_glm(glm_sim.distr, beta0, beta, X)
+    y = np.zeros((n_samples, 2))
+    for idx, beta in enumerate(betas.T):
+        y[:, idx] = simulate_glm(glm_sim.distr, beta0, beta, X)
     y = np.argmax(y, axis=1)
     glm_sim.fit(X, y)
     score = glm_sim.score(X, y)
diff --git a/tests/test_pyglmnet.py b/tests/test_pyglmnet.py
index 84b9726..43a098f 100644
--- a/tests/test_pyglmnet.py
+++ b/tests/test_pyglmnet.py
@@ -120,7 +120,7 @@ def test_group_lasso():
     groups[60:] = 3
 
     # sample random coefficients
-    beta0 = np.random.normal(0.0, 1.0, 1)
+    beta0 = np.random.rand()
     beta = np.random.normal(0.0, 1.0, n_features)
     beta[groups == 2] = 0.
 
@@ -155,7 +155,8 @@ def test_group_lasso():
 
 @pytest.mark.parametrize("distr", ALLOWED_DISTRS)
 @pytest.mark.parametrize("reg_lambda", [0.0, 0.1])
-def test_glmnet(distr, reg_lambda):
+@pytest.mark.parametrize("fit_intercept", [True, False])
+def test_glmnet(distr, reg_lambda, fit_intercept):
     """Test glmnet."""
     raises(ValueError, GLM, distr='blah')
     raises(ValueError, GLM, distr='gaussian', max_iter=1.8)
@@ -163,9 +164,11 @@ def test_glmnet(distr, reg_lambda):
     n_samples, n_features = 100, 10
 
     # coefficients
-    beta0 = 1. / (np.float(n_features) + 1.) * \
-        np.random.normal(0.0, 1.0)
-    beta = 1. / (np.float(n_features) + 1.) * \
+    beta0 = 0.
+    if fit_intercept:
+        beta0 = 1. / (np.float(n_features) + 1.) * \
+            np.random.normal(0.0, 1.0)
+    beta = 1. / (np.float(n_features) + int(fit_intercept)) * \
         np.random.normal(0.0, 1.0, (n_features,))
 
     solvers = ['batch-gradient', 'cdfast']
@@ -196,12 +199,14 @@ def test_glmnet(distr, reg_lambda):
             Tau = None
             loss_trace.append(
                 _loss(distr, alpha, Tau, reg_lambda,
-                      X_train, y_train, eta, group, beta))
+                      X_train, y_train, eta, group, beta,
+                      fit_intercept=fit_intercept))
 
         glm = GLM(distr, learning_rate=learning_rate,
                   reg_lambda=reg_lambda, tol=1e-3, max_iter=5000,
                   alpha=alpha, solver=solver, score_metric=score_metric,
-                  random_state=random_state, callback=callback)
+                  random_state=random_state, callback=callback,
+                  fit_intercept=fit_intercept)
         assert(repr(glm))
 
         glm.fit(X_train, y_train)
@@ -324,6 +329,13 @@ def test_compare_sklearn(solver):
     y_pg = glm.predict(X)
     assert abs(rmse(Y, y_sk) - rmse(Y, y_pg)) < 1.0
 
+    glm = GLM(distr='gaussian', alpha=l1_ratio, reg_lambda=alpha,
+              solver=solver, tol=1e-5, max_iter=5, fit_intercept=False)
+    glm.fit(X, Y)
+    assert glm.beta0_ == 0.
+
+    glm.predict(X)
+
 
 @pytest.mark.parametrize("distr", ALLOWED_DISTRS)
 def test_cdfast(distr):
@@ -430,15 +442,21 @@ def test_simulate_glm(distr):
     n_samples, n_features = 10, 3
 
     # sample random coefficients
-    beta0 = state.normal(0.0, 1.0, 1)
+    beta0 = state.rand()
     beta = state.normal(0.0, 1.0, n_features)
 
     X = state.normal(0.0, 1.0, [n_samples, n_features])
     simulate_glm(distr, beta0, beta, X, random_state=random_state)
 
+    with pytest.raises(ValueError, match="'beta0' must be float"):
+        simulate_glm(distr, np.array([1.0]), beta, X, random_state)
+
+    with pytest.raises(ValueError, match="'beta' must be 1D"):
+        simulate_glm(distr, 1.0, np.atleast_2d(beta), X, random_state)
+
     # If the distribution name is garbage it will fail
     distr = 'multivariate_gaussian_poisson'
-    with pytest.raises(ValueError):
+    with pytest.raises(ValueError, match="'distr' must be in"):
         simulate_glm(distr, 1.0, 1.0, np.array([[1.0]]))
 
 
@@ -471,6 +489,9 @@ def test_api_input():
     with pytest.raises(ValueError, match="solver must be one of"):
         glm.fit(X, y)
 
+    with pytest.raises(ValueError, match="fit_intercept must be"):
+        glm = GLM(distr='gaussian', fit_intercept='blah')
+
     glm = GLM(distr='gaussian', max_iter=2)
     with pytest.warns(UserWarning, match='Reached max number of iterat'):
         glm.fit(X, y)