diff --git a/glimix_core/__init__.py b/glimix_core/__init__.py
index 55c6e870..e7acf5df 100644
--- a/glimix_core/__init__.py
+++ b/glimix_core/__init__.py
@@ -12,7 +12,7 @@
 from . import cov, ggp, glmm, gp, lik, link, lmm, mean, random, util
 from .testit import test
 
-__version__ = "1.3.8"
+__version__ = "1.4.0"
 
 __all__ = [
     "__version__", "test", "ggp", "gp", "lmm", "glmm", "cov", "lik", "mean",
diff --git a/glimix_core/glmm/expfam.py b/glimix_core/glmm/expfam.py
index 0c040412..5e47cb42 100644
--- a/glimix_core/glmm/expfam.py
+++ b/glimix_core/glmm/expfam.py
@@ -2,11 +2,10 @@
 
 from copy import copy
 
+from glimix_core.ep import EPLinearKernel
 from liknorm import LikNormMachine
-from numpy import exp
-from numpy import ascontiguousarray
+from numpy import ascontiguousarray, exp
 
-from glimix_core.ep import EPLinearKernel
 from .glmm import GLMM
 
 
@@ -67,12 +66,12 @@ class GLMMExpFam(GLMM):
         After: -13.43
     """
 
-    def __init__(self, y, lik_name, X, QS):
+    def __init__(self, y, lik_name, X, QS, n_int=1000):
         GLMM.__init__(self, y, lik_name, X, QS)
 
         self._ep = EPLinearKernel(self._X.shape[0])
         self._ep.set_compute_moments(self.compute_moments)
-        self._machine = LikNormMachine(self._lik_name, 1000)
+        self._machine = LikNormMachine(self._lik_name, n_int)
         self.update_approx = True
 
         self.variables().get('beta').listen(self.set_update_approx)
diff --git a/glimix_core/glmm/glmm.py b/glimix_core/glmm/glmm.py
index bd593bc6..bc3d665b 100644
--- a/glimix_core/glmm/glmm.py
+++ b/glimix_core/glmm/glmm.py
@@ -7,8 +7,8 @@
 from numpy_sugar.linalg import ddot, sum2diag
 from optimix import Function, Scalar, Vector
 
-from ..util import check_covariates, check_economic_qs, check_outcome
-from ..util import normalise_outcome
+from ..util import (check_covariates, check_economic_qs, check_outcome,
+                    normalise_outcome)
 
 
 class GLMM(Function):
@@ -151,7 +151,7 @@ def fix(self, var_name):
         """
         Function.fix(self, _to_internal_name(var_name))
 
-    def fit(self, verbose=True):
+    def fit(self, verbose=True, factr=1e5, pgtol=1e-7):
         r"""Maximise the marginal likelihood.
 
         Parameters
@@ -160,7 +160,10 @@ def fit(self, verbose=True):
             ``True`` for progress output; ``False`` otherwise.
             Defaults to ``True``.
         """
-        self.feed().maximize(verbose=verbose)
+        f = self.feed()
+        f.factr = factr
+        f.pgtol = pgtol
+        f.maximize(verbose=verbose)
 
     def lml(self):
         r"""Log of the marginal likelihood.