diff --git a/.flake8 b/.flake8
new file mode 100644
index 0000000..2b368e1
--- /dev/null
+++ b/.flake8
@@ -0,0 +1,11 @@
+[flake8]
+max-line-length = 99
+ignore = W504, W503, E266
+exclude =
+    .git,
+    __pycache__,
+    __init__.py,
+    scratch_*.py,
+    venv,
+    build,
+    dist
diff --git a/Examples.ipynb b/Examples.ipynb
index 7941650..248c94a 100644
--- a/Examples.ipynb
+++ b/Examples.ipynb
@@ -12,7 +12,7 @@
     "Example 3: Same, with two different lapse rates<br>\n",
     "Example 4: Fit data from .5 to 1, using weibull\n",
     "\n",
-    "### Advanced: \n",
+    "### Advanced:\n",
     "Add error bars to data, parameters, and fits"
    ]
   },
@@ -35,7 +35,7 @@
    "metadata": {},
    "source": [
     "## Example 1\n",
-    "Fit data from 0 to 1 and stimulus in log units, using erf "
+    "Fit data from 0 to 1 and stimulus in log units, using erf"
    ]
   },
   {
@@ -71,8 +71,8 @@
     "ntrials = 40\n",
     "dd = np.random.binomial(1, pp, size=(ntrials,nxx))\n",
     "\n",
-    "# data: \n",
-    "# 3 x n matrix where first row corrsponds to stim levels (log units), \n",
+    "# data:\n",
+    "# 3 x n matrix where first row corresponds to stim levels (log units),\n",
     "# the second to number of trials for each stim level (int),\n",
     "# the third to proportion correct (float between 0 and 1)\n",
     "data = np.vstack((np.log10(xx), 10 * np.ones((nxx,)), np.mean(dd, axis=0)))\n",
@@ -141,8 +141,8 @@
     "    'nfits': 10\n",
     "}\n",
     "\n",
-    "# data: \n",
-    "# 3 x n matrix where first row corrsponds to stim levels (% contrast), \n",
+    "# data:\n",
+    "# 3 x n matrix where first row corresponds to stim levels (% contrast),\n",
     "# the second to number of trials for each stim level (int),\n",
     "# the third to proportion rightward (float between 0 and 1)\n",
     "data = np.vstack((xx, ntrials * np.ones((nxx,)), np.mean(dd, axis=0)))\n",
@@ -215,8 +215,8 @@
     "    'nfits': 10\n",
     "}\n",
     "\n",
-    "# data: \n",
-    "# 3 x n matrix where first row corrsponds to stim levels (% contrast), \n",
+    "# data:\n",
+    "# 3 x n matrix where first row corresponds to stim levels (% contrast),\n",
     "# the second to number of trials for each stim level (int),\n",
     "# the third to proportion rightward (float between 0 and 1)\n",
     "data = np.vstack((xx, ntrials * np.ones((nxx,)), np.mean(dd, axis=0)))\n",
@@ -275,12 +275,12 @@
     "ntrials = 80\n",
     "dd = np.random.binomial(1., pp, size=(ntrials, nxx))\n",
     "\n",
-    "# data: \n",
-    "# 3 x n matrix where first row corrsponds to stim levels (notice we do NOT take log of x values), \n",
+    "# data:\n",
+    "# 3 x n matrix where first row corresponds to stim levels (notice we do NOT take log of x values),\n",
     "# the second to number of trials for each stim level (int),\n",
     "# the third to proportion correct (float between 0 and 1)\n",
     "data = np.vstack((xx, ntrials * np.ones((nxx,)), np.mean(dd, axis=0)))\n",
-    "# fit to reconstruct the parameters \n",
+    "# fit to reconstruct the parameters\n",
     "pars, L = psy.mle_fit_psycho(data, 'weibull50');\n",
     "\n",
     "# graphics\n",
@@ -316,4 +316,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
\ No newline at end of file
diff --git a/psychofit.py b/psychofit.py
index 4e15352..5199b27 100644
--- a/psychofit.py
+++ b/psychofit.py
@@ -7,17 +7,15 @@
 The data can be expressed in fraction correct (from .5 to 1) or in
 fraction of one specific choice (from 0 to 1). To fit them you can use
 these functions:
-  weibull50          - Weibull function from 0.5 to 1, with lapse rate
-  weibull            - Weibull function from 0 to 1, with lapse rate
-  erf_psycho         - erf function from 0 to 1, with lapse rate
-  erf_psycho_2gammas - erf function from 0 to 1, with two lapse rates
-
+  - weibull50:          Weibull function from 0.5 to 1, with lapse rate
+  - weibull:            Weibull function from 0 to 1, with lapse rate
+  - erf_psycho:         erf function from 0 to 1, with lapse rate
+  - erf_psycho_2gammas: erf function from 0 to 1, with two lapse rates
 Functions in the toolbox are:
-  mle_fit_psycho     - Maximumum likelihood fit of psychometric function
-  neg_likelihood     - Negative likelihood of a psychometric function
-
+  - mle_fit_psycho:     Maximum likelihood fit of psychometric function
+  - neg_likelihood:     Negative likelihood of a psychometric function
 For more info, see:
-  Examples           - Examples of use of psychofit toolbox
+  - Examples:           Examples of use of psychofit toolbox
 
 Matteo Carandini, 2000-2015
 """
@@ -85,7 +83,7 @@ def mle_fit_psycho(data, P_model='weibull', parstart=None, parmin=None, parmax=N
     if data.shape[0] != 3:
         raise ValueError('data must be m by 3 matrix')
 
-    rep = lambda x: (x, x) if P_model.endswith('2gammas') else (x,)
+    rep = lambda x: (x, x) if P_model.endswith('2gammas') else (x,)  # noqa
     if parstart is None:
         parstart = np.array([np.mean(data[0, :]), 3., *rep(.05)])
     if parmin is None:
@@ -103,7 +101,7 @@ def mle_fit_psycho(data, P_model='weibull', parstart=None, parmin=None, parmax=N
                           P_model=P_model, parmin=parmin, parmax=parmax)
     for ifit in range(nfits):
         pars[ifit, :] = scipy.optimize.fmin(f, parstart, disp=False)
-        parstart = parmin + np.random.rand(parmin.size) * (parmax-parmin)
+        parstart = parmin + np.random.rand(parmin.size) * (parmax - parmin)
         likelihoods[ifit] = -neg_likelihood(pars[ifit, :], data[:, ii], P_model, parmin, parmax)
 
     # the values to be output
@@ -128,7 +126,7 @@ def neg_likelihood(pars, data, P_model='weibull', parmin=None, parmax=None):
         parmax: Maximum bound for parameters.  If None, some reasonable defaults are used
 
     Returns:
-        l: The likelihood of the parameters.  The equation is:
+        ll: The likelihood of the parameters.  The equation is:
             - sum(nn.*(pp.*log10(P_model)+(1-pp).*log10(1-P_model)))
             See the the appendix of Watson, A.B. (1979). Probability
             summation over time. Vision Res 19, 515-522.
@@ -165,8 +163,8 @@ def neg_likelihood(pars, data, P_model='weibull', parmin=None, parmax=None):
 
     # here is where you effectively put the constraints.
     if (any(pars < parmin)) or (any(pars > parmax)):
-        l = 10000000
-        return l
+        ll = 10000000
+        return ll
 
     dispatcher = {
         'weibull': weibull,
@@ -186,8 +184,8 @@ def neg_likelihood(pars, data, P_model='weibull', parmin=None, parmax=None):
     probs[probs == 0] = np.finfo(float).eps
     probs[probs == 1] = 1 - np.finfo(float).eps
 
-    l = - sum(nn * (pp * np.log(probs) + (1 - pp) * np.log(1 - probs)))
-    return l
+    ll = - sum(nn * (pp * np.log(probs) + (1 - pp) * np.log(1 - probs)))
+    return ll
 
 
 def weibull(pars, xx):
@@ -218,7 +216,7 @@ def weibull(pars, xx):
         raise ValueError('pars must be of length 3')
 
     alpha, beta, gamma = pars
-    return (1 - gamma) - (1 - 2*gamma) * np.exp(-((xx / alpha)**beta))
+    return (1 - gamma) - (1 - 2 * gamma) * np.exp(-((xx / alpha) ** beta))
 
 
 def weibull50(pars, xx):
diff --git a/psychofit_tests.py b/psychofit_tests.py
index b854d98..20a03d2 100644
--- a/psychofit_tests.py
+++ b/psychofit_tests.py
@@ -110,9 +110,9 @@ def test_neg_likelihood(self):
         self.assertRaises(TypeError, psy.neg_likelihood, '(10, 20, .05)', None)
         self.assertRaises(ValueError, psy.neg_likelihood, (.5, 10, .05), data, P_model='foo')
 
-        l = psy.neg_likelihood((-20, 30, 2), data.tolist(), P_model='erf_psycho',
-                               parmin=np.array((-10, 20, 0)), parmax=np.array((10, 10, .05)))
-        self.assertTrue(l > 10000)
+        ll = psy.neg_likelihood((-20, 30, 2), data.tolist(), P_model='erf_psycho',
+                                parmin=np.array((-10, 20, 0)), parmax=np.array((10, 10, .05)))
+        self.assertTrue(ll > 10000)
 
     def test_mle_fit_psycho(self):
         expected = {