diff --git a/docs/source/pybambi.neuralnetworks.rst b/docs/source/pybambi.neuralnetworks.rst
index a39a0a8..cd4bcd9 100644
--- a/docs/source/pybambi.neuralnetworks.rst
+++ b/docs/source/pybambi.neuralnetworks.rst
@@ -20,10 +20,10 @@ pybambi.neuralnetworks.kerasnet module
     :undoc-members:
     :show-inheritance:
 
-pybambi.neuralnetworks.nearestneighbor module
----------------------------------------------
+pybambi.neuralnetworks.nearestneighbour module
+----------------------------------------------
 
-.. automodule:: pybambi.neuralnetworks.nearestneighbor
+.. automodule:: pybambi.neuralnetworks.nearestneighbour
     :members:
     :undoc-members:
     :show-inheritance:
diff --git a/docs/source/pybambi.rst b/docs/source/pybambi.rst
index b3eeb5c..babafe8 100644
--- a/docs/source/pybambi.rst
+++ b/docs/source/pybambi.rst
@@ -19,10 +19,10 @@ pybambi.bambi module
     :undoc-members:
     :show-inheritance:
 
-pybambi.dumper module
----------------------
+pybambi.manager module
+----------------------
 
-.. automodule:: pybambi.dumper
+.. automodule:: pybambi.manager
     :members:
     :undoc-members:
     :show-inheritance:
diff --git a/pybambi/bambi.py b/pybambi/bambi.py
index d1fdf0e..85bb395 100644
--- a/pybambi/bambi.py
+++ b/pybambi/bambi.py
@@ -4,7 +4,7 @@
 Date: November 2018
 """
 import os
-from pybambi.dumper import dumper
+from pybambi.manager import BambiManager
 
 
 def run_pyBAMBI(loglikelihood, prior, nDims, **kwargs):
@@ -32,6 +32,22 @@ def run_pyBAMBI(loglikelihood, prior, nDims, **kwargs):
         efficiency for multinest.
         Default `0.5**nDims`
 
+    learner: object
+        information indicating what learning algorithm to use for approximating
+        the likelihood. Can be the string `'keras'`, or a `keras.models.Model`
+        Default `'keras'`
+
+    ntrain: int
+        Number of training points to use
+        Default `nlive`
+
+    proxy_tolerance: float
+        Required accuracy of proxy.
+        Default `0.01`
+
+    ns_output: int
+        Nested sampling output level.
+
     """
     # Process kwargs
     nested_sampler = kwargs.pop('nested_sampler', 'polychord')
@@ -39,20 +55,28 @@ def run_pyBAMBI(loglikelihood, prior, nDims, **kwargs):
     root = kwargs.pop('root', os.path.join('chains', nested_sampler))
     num_repeats = kwargs.pop('num_repeats', nDims*5)
     eff = kwargs.pop('eff', 0.5**nDims)
+    learner = kwargs.pop('learner', 'keras')
+    proxy_tolerance = kwargs.pop('proxy_tolerance', 0.1)
+    failure_tolerance = kwargs.pop('failure_tolerance', 0.5)
+    ntrain = kwargs.pop('ntrain', nlive)
 
     if kwargs:
         raise TypeError('Unexpected **kwargs: %r' % kwargs)
 
+    # Set up the global manager of the BAMBI session.
+    thumper = BambiManager(loglikelihood, learner, proxy_tolerance,
+                           failure_tolerance, ntrain)
+
     # Choose and run sampler
     if nested_sampler == 'polychord':
         from pybambi.polychord import run_polychord
-        run_polychord(loglikelihood, prior, dumper, nDims,
-                      nlive, root, num_repeats)
+        run_polychord(thumper.loglikelihood, prior, thumper.dumper, nDims,
+                      nlive, root, ntrain//2, num_repeats)
 
     elif nested_sampler == 'multinest':
         from pybambi.multinest import run_multinest
-        run_multinest(loglikelihood, prior, dumper, nDims,
-                      nlive, root, eff)
+        run_multinest(thumper.loglikelihood, prior, thumper.dumper, nDims,
+                      nlive, root, ntrain//2, eff)
 
     else:
         raise NotImplementedError('nested sampler %s is not implemented'
diff --git a/pybambi/dumper.py b/pybambi/dumper.py
deleted file mode 100644
index 3ae59b1..0000000
--- a/pybambi/dumper.py
+++ /dev/null
@@ -1,25 +0,0 @@
-"""Function giving access to live points.
-
-Author: Will Handley (wh260@cam.ac.uk)
-Date: November 2018
-"""
-
-
-def dumper(live):
-    """Dumper function giving access to the live points.
-
-    Parameters
-    ----------
-    live:
-        `numpy.array` of live parameters and loglikelihoods,
-        `shape=(nlive,nDims+1)`
-
-    """
-    params = live[:, :-1]
-    loglikes = live[:, -1]
-
-    print("-----------------------------")
-    print("Call neural network code here")
-    print(params.shape)
-    print(loglikes)
-    print("-----------------------------")
diff --git a/pybambi/manager.py b/pybambi/manager.py
new file mode 100644
index 0000000..3c23659
--- /dev/null
+++ b/pybambi/manager.py
@@ -0,0 +1,98 @@
+"""BAMBI management object.
+
+Author: Pat Scott (p.scott@imperial.ac.uk)
+Date: Feb 2019
+"""
+
+import numpy as np
+
+from pybambi.neuralnetworks.kerasnet import KerasNetInterpolation
+from pybambi.neuralnetworks.nearestneighbour \
+        import NearestNeighbourInterpolation
+import keras.models
+
+
+class BambiManager(object):
+    """Does all the talking for BAMBI.
+
+    Takes a new set of training data from the dumper and trains (or retrains) a
+    neural net, and assesses whether or not it can be used for a given
+    parameter combination.
+
+    Parameters
+    ----------
+    ntrain: int
+        Number of training points to use
+
+    """
+
+    def __init__(self, loglikelihood, learner, proxy_tolerance,
+                 failure_tolerance, ntrain):
+        """Construct bambi object."""
+        self.proxy_tolerance = proxy_tolerance
+        self._loglikelihood = loglikelihood
+        self._learner = learner
+        self._proxy_tolerance = proxy_tolerance
+        self._failure_tolerance = failure_tolerance
+        self._ntrain = ntrain
+        self._proxy_trained = False
+        self.old_learners = []
+
+    def make_learner(self, params, loglikes):
+        """Construct a Predictor."""
+        if self._learner == 'keras':
+            return KerasNetInterpolation(params, loglikes)
+        elif self._learner == 'nearestneighbour':
+            return NearestNeighbourInterpolation(params, loglikes)
+        elif issubclass(type(self._learner), keras.models.Model):
+            return KerasNetInterpolation(params, loglikes, model=self._learner)
+        else:
+            raise NotImplementedError('learner %s is not implemented.'
+                                      % self._learner)
+
+    def dumper(self, live_params, live_loglks, dead_params, dead_loglks):
+        """Respond to signal from nested sampler."""
+        if not self._proxy_trained:
+            params = np.concatenate((live_params, dead_params))
+            loglikes = np.concatenate((live_loglks, dead_loglks))
+            self.train_new_learner(params[:self._ntrain, :],
+                                   loglikes[:self._ntrain])
+        if self._proxy_trained:
+            print("Using trained proxy")
+        else:
+            print("Unable to use proxy")
+
+    def loglikelihood(self, params):
+        """Bambi Proxy wrapper for original loglikelihood."""
+        # Short circuit to the full likelihood if proxy not yet fully trained
+        if not self._proxy_trained:
+            return self._loglikelihood(params)
+
+        # Call the learner
+        candidate_loglikelihood = self._current_learner(params)
+
+        # If the learner can be trusted, use its estimate,
+        # otherwise use the original like and update the failure status
+        if self._current_learner.valid(candidate_loglikelihood):
+            return candidate_loglikelihood
+        else:
+            self._rolling_failure_fraction = (1.0 + (self._ntrain - 1.0) *
+                                              self._rolling_failure_fraction
+                                              ) / self._ntrain
+            if self._rolling_failure_fraction > self._failure_tolerance:
+                self._proxy_trained = False
+            return self._loglikelihood(params)
+
+    def train_new_learner(self, params, loglikes):
+        """Train a new Predictor."""
+        try:
+            self.old_learners.append(self._current_learner)
+        except AttributeError:
+            pass
+        self._current_learner = self.make_learner(params, loglikes)
+        sigma = self._current_learner.uncertainty()
+        print("Current uncertainty in network log-likelihood predictions: %s"
+              % sigma)
+        if sigma < self._proxy_tolerance:
+            self._proxy_trained = True
+            self._rolling_failure_fraction = 0.0
diff --git a/pybambi/multinest.py b/pybambi/multinest.py
index 8c05a24..9d576a6 100644
--- a/pybambi/multinest.py
+++ b/pybambi/multinest.py
@@ -6,7 +6,8 @@
 from numpy.ctypeslib import as_array
 
 
-def run_multinest(loglikelihood, prior, dumper, nDims, nlive, root, eff):
+def run_multinest(loglikelihood, prior, dumper, nDims, nlive, root, ndump,
+                  eff):
     """Run MultiNest.
 
     See https://arxiv.org/abs/0809.3437 for more detail
@@ -46,6 +47,9 @@ def run_multinest(loglikelihood, prior, dumper, nDims, nlive, root, eff):
     root: str
         base name for output files
 
+    ndump: int
+        How many iterations between dumper function calls
+
     eff: float
         Efficiency of MultiNest
 
@@ -53,7 +57,9 @@ def run_multinest(loglikelihood, prior, dumper, nDims, nlive, root, eff):
     import pymultinest
 
     def multinest_prior(cube, ndim, nparams):
-        return prior(as_array(cube, shape=(nparams,)))
+        theta = prior(as_array(cube, shape=(nparams,)))
+        for i, elem in enumerate(theta):
+            cube[i] = elem
 
     def multinest_loglikelihood(cube, ndim, nparams):
         return loglikelihood(as_array(cube, shape=(nparams,)))
@@ -61,9 +67,11 @@ def multinest_loglikelihood(cube, ndim, nparams):
     def multinest_dumper(nSamples, nlive, nPar,
                          physLive, posterior, paramConstr,
                          maxLogLike, logZ, logZerr, nullcontext):
-        dumper(physLive)
+        dumper(physLive[:, :-1], physLive[:, -1],
+               posterior[:, :-2], posterior[:, -2])
 
     pymultinest.run(multinest_loglikelihood, multinest_prior, nDims,
                     resume=False, verbose=True, dump_callback=multinest_dumper,
-                    n_iter_before_update=nlive//10, n_live_points=nlive,
-                    outputfiles_basename=root, sampling_efficiency=eff)
+                    n_iter_before_update=ndump//10, n_live_points=nlive,
+                    outputfiles_basename=root, sampling_efficiency=eff,
+                    evidence_tolerance=0.01)
diff --git a/pybambi/neuralnetworks/base.py b/pybambi/neuralnetworks/base.py
index 922eaa6..ea8b3ec 100644
--- a/pybambi/neuralnetworks/base.py
+++ b/pybambi/neuralnetworks/base.py
@@ -15,7 +15,7 @@ class Predictor(object):
     Parameters
     ----------
     params:
-        `numpy.array of` physical parameters to train on
+        `numpy.array` of physical parameters to train on
         shape (ntrain, ndims)
 
     logL:
@@ -24,10 +24,11 @@ class Predictor(object):
 
     """
 
-    def __init__(self, params, logL):
+    def __init__(self, params, logL, split=0.8):
         """Construct predictor from training data."""
         params = numpy.array(params)
         logL = numpy.array(logL)
+
         if len(params) != len(logL):
             raise ValueError("input and target must be the same length")
         elif params.ndim != 2:
@@ -35,6 +36,18 @@ def __init__(self, params, logL):
         elif logL.ndim != 1:
             raise ValueError("target must be one-dimensional")
 
+        nparams = len(params)
+        randomize = numpy.random.permutation(nparams)
+        params = params[randomize]
+        logL = logL[randomize]
+
+        self._maxLogL = numpy.max(logL)
+        self._minLogL = numpy.min(logL)
+        ntrain = int(split*nparams)
+        indx = [ntrain]
+        self.params_training, self.params_testing = numpy.split(params, indx)
+        self.logL_training, self.logL_testing = numpy.split(logL, indx)
+
     def __call__(self, x):
         """Calculate proxy loglikelihood.
 
@@ -50,3 +63,26 @@ def __call__(self, x):
         """
         err = "Predictor: You need to implement a call function"
         raise NotImplementedError(err)
+
+    def uncertainty(self):
+        """Uncertainty value for the trained model."""
+        err = "Predictor: You need to implement an uncertainty function"
+        raise NotImplementedError(err)
+
+    def valid(self, loglikelihood):
+        """Check validity of proxy.
+
+        Checks to see if the supplied log likelihood value is within the
+        current range of likelihoods, including the uncertainty
+
+        Parameters
+        ----------
+        loglikelihood:
+            Value of the log likelihood that needs checking
+
+        """
+        inRange = True
+        if loglikelihood > self._maxLogL + self.uncertainty() \
+                or loglikelihood < self._minLogL - self.uncertainty():
+            inRange = False
+        return inRange
diff --git a/pybambi/neuralnetworks/kerasnet.py b/pybambi/neuralnetworks/kerasnet.py
index 27160c9..c6eafc9 100644
--- a/pybambi/neuralnetworks/kerasnet.py
+++ b/pybambi/neuralnetworks/kerasnet.py
@@ -10,6 +10,7 @@
 from pybambi.neuralnetworks.base import Predictor
 from keras.models import Sequential
 from keras.layers import Dense
+from keras.callbacks import EarlyStopping
 
 
 class KerasNetInterpolation(Predictor):
@@ -31,38 +32,38 @@ class KerasNetInterpolation(Predictor):
 
     """
 
-    def __init__(self, params, logL, split=0.8):
+    def __init__(self, params, logL, split=0.8, model=None):
         """Construct predictor from training data."""
-        super(KerasNetInterpolation, self).__init__(params, logL)
-        self._params = params[:]
-        self._logL = logL[:]
-
-        # This function takes a parameter set, trains a neural net and returns
-        # the interpolated likelihood
-        # It is intended to be cannibalised for later pyBAMBI functionality
+        super(KerasNetInterpolation, self).__init__(params, logL, split)
+
+        if model is None:
+            self.model = self._default_architecture()
+        else:
+            self.model = model
+
+        callbacks = [EarlyStopping(monitor='val_loss', mode='min',
+                                   min_delta=0.001,
+                                   patience=10,
+                                   restore_best_weights=True)]
+
+        self.history = self.model.fit(self.params_training,
+                                      self.logL_training,
+                                      validation_data=(self.params_testing,
+                                                       self.logL_testing),
+                                      epochs=300,
+                                      callbacks=callbacks)
+
+    def _default_architecture(self):
+        # Create model
+        model = Sequential()
 
         # Number of neurons in each hidden layer, could make this configurable?
         numNeurons = 200
 
-        # Shuffle the params and logL in unison
-        # (important for splitting data into training and test sets)
-        ntot = len(self._params)
-        randomize = numpy.random.permutation(ntot)
-        params = self._params[randomize]
-        logL = self._logL[randomize]
-
-        # Now split into training and test sets
-        ntrain = int(split*ntot)
-        params_training, params_test = numpy.split(params, [ntrain])
-        logL_training, logL_test = numpy.split(logL, [ntrain])
-
-        # Create model
-        model = Sequential()
-
         # Get number of input parameters
         # Note: if params contains extra quantities (ndim+others),
         # we need to change this
-        n_cols = params.shape[1]
+        n_cols = self.params_training.shape[1]
 
         # Add model layers, note choice of activation function (relu)
         # We will use 3 hidden layers and an output layer
@@ -78,10 +79,7 @@ def __init__(self, params, logL, split=0.8):
         # Need to choose training optimiser, and the loss function
         model.compile(optimizer='adam', loss='mean_squared_error')
 
-        self._history = model.fit(params_training, logL_training,
-                                  validation_data=(params_test, logL_test),
-                                  epochs=100)
-        self._model = model
+        return model
 
     def __call__(self, x):
         """Calculate proxy loglikelihood.
@@ -97,5 +95,11 @@ def __call__(self, x):
 
         """
         x_ = numpy.atleast_2d(x)
-        y = self._model.predict(x_)
-        return numpy.squeeze(y)
+        y = self.model.predict(x_)
+        return float(numpy.squeeze(y))
+
+    def uncertainty(self):
+        """Uncertainty value for the trained keras model."""
+        test_loss = numpy.sqrt(self.history.history['val_loss'])
+
+        return numpy.squeeze(test_loss.min())
diff --git a/pybambi/neuralnetworks/nearestneighbor.py b/pybambi/neuralnetworks/nearestneighbour.py
similarity index 54%
rename from pybambi/neuralnetworks/nearestneighbor.py
rename to pybambi/neuralnetworks/nearestneighbour.py
index 9ddc145..6dc2aea 100644
--- a/pybambi/neuralnetworks/nearestneighbor.py
+++ b/pybambi/neuralnetworks/nearestneighbour.py
@@ -1,18 +1,18 @@
-"""Nearest neighbor interpolation predictor.
+"""Nearest neighbour interpolation predictor.
 
 Author: Will Handley (wh260@cam.ac.uk)
 Date: November 2018
 
 
-This implements a nearest neighbor interpolation, and is designed as a
+This implements a nearest neighbour interpolation, and is designed as a
 placeholder predictor, rather than an actual neural network
 """
 import numpy
 from pybambi.neuralnetworks.base import Predictor
 
 
-class NearestNeighborInterpolation(Predictor):
-    """Nearest Neighbor interpolation.
+class NearestNeighbourInterpolation(Predictor):
+    """Nearest Neighbour interpolation.
 
     Returns the loglikelihood of the training point closest in parameter space
 
@@ -28,11 +28,13 @@ class NearestNeighborInterpolation(Predictor):
 
     """
 
-    def __init__(self, params, logL):
+    def __init__(self, params, logL, split=0.8):
         """Construct predictor from training data."""
-        super(NearestNeighborInterpolation, self).__init__(params, logL)
-        self._params = params[:]
-        self._logL = logL[:]
+        super(NearestNeighbourInterpolation,
+              self).__init__(params, logL, split)
+
+        self.std = numpy.std([self(par)-logL for par, logl in
+                              zip(self.params_testing, self.logL_testing)])
 
     def __call__(self, x):
         """Calculate proxy loglikelihood.
@@ -47,6 +49,10 @@ def __call__(self, x):
         proxy loglikelihood value(s)
 
         """
-        distances = numpy.linalg.norm(self._params - x, axis=1)
+        distances = numpy.linalg.norm(self.params_training - x, axis=1)
         i = numpy.argmin(distances)
-        return self._logL[i]
+        return self.logL_training[i]
+
+    def uncertainty(self):
+        """Rough uncertainty for the nearest neighbour model."""
+        return self.std
diff --git a/pybambi/polychord.py b/pybambi/polychord.py
index eac429b..5ddec70 100644
--- a/pybambi/polychord.py
+++ b/pybambi/polychord.py
@@ -4,9 +4,10 @@
 Date: November 2018
 """
 import os
+import numpy
 
 
-def run_polychord(loglikelihood, prior, dumper, nDims, nlive, root,
+def run_polychord(loglikelihood, prior, dumper, nDims, nlive, root, ndump,
                   num_repeats):
     """Run PolyChord.
 
@@ -47,6 +48,9 @@ def run_polychord(loglikelihood, prior, dumper, nDims, nlive, root,
     root: str
         base name for output files
 
+    ndump: int
+        How many iterations between dumper function calls
+
     num_repeats: int
         Length of chain to generate new live points
 
@@ -62,12 +66,14 @@ def run_polychord(loglikelihood, prior, dumper, nDims, nlive, root,
     settings.num_repeats = num_repeats
     settings.do_clustering = True
     settings.read_resume = False
+    settings.compression_factor = numpy.exp(-float(ndump)/nlive)
+    settings.precision_criterion = 0.01
 
     def polychord_loglikelihood(theta):
         return loglikelihood(theta), []
 
     def polychord_dumper(live, dead, logweights, logZ, logZerr):
-        dumper(live[:, :-1])
+        dumper(live[:, :-2], live[:, -1], dead[:, :-2], dead[:, -1])
 
     pypolychord.run_polychord(polychord_loglikelihood, nDims, nDerived,
                               settings, prior, polychord_dumper)
diff --git a/test.py b/test.py
index ad19b9d..9f0d4ea 100644
--- a/test.py
+++ b/test.py
@@ -9,8 +9,8 @@ def loglikelihood(theta):
     """ Spherical Gaussian Likelihood """
     sigma = 0.1
     nDims = len(theta)
-    logL = -log(2*pi*sigma*sigma)*nDims/2.0
-    logL -= sum((theta/sigma)**2) / 2 + log(2) * nDims
+    logL = -sum((theta/sigma)**2) / 2
+    logL += -log(2*pi*sigma*sigma)*nDims/2.0 + log(2) * nDims
     return logL
 
 
@@ -19,4 +19,7 @@ def prior(cube):
     return -1 + 2 * cube
 
 
-run_pyBAMBI(loglikelihood, prior, nDims, nested_sampler='polychord', nlive=500)
+run_pyBAMBI(loglikelihood, prior, nDims,
+            nested_sampler='multinest',
+            nlive=500,
+            learner='keras')
diff --git a/tests/neuralnetworks/test_kerasnet.py b/tests/neuralnetworks/test_kerasnet.py
index 59253ee..02e51fa 100644
--- a/tests/neuralnetworks/test_kerasnet.py
+++ b/tests/neuralnetworks/test_kerasnet.py
@@ -19,8 +19,6 @@ def test_KerasNet():
     # Initialise the keras net example, which includes training
     p = KerasNetInterpolation(params, logL)
 
-    # Test 1D input has better than 5% accuracy
-    assert numpy.abs((p(params[0]) - logL[0]) / logL[0]) < 0.05
-
-    # Test 2D input has better than 5% accuracy
-    assert numpy.abs((p(params) - logL) / logL).max() < 0.05
+    # Test input has better than 5% accuracy
+    for i, l in enumerate(logL):
+        assert abs(p(params[i]) - l) < 1.
diff --git a/tests/neuralnetworks/test_nearestneighbor.py b/tests/neuralnetworks/test_nearestneighbour.py
similarity index 57%
rename from tests/neuralnetworks/test_nearestneighbor.py
rename to tests/neuralnetworks/test_nearestneighbour.py
index 174bde9..1381a2c 100644
--- a/tests/neuralnetworks/test_nearestneighbor.py
+++ b/tests/neuralnetworks/test_nearestneighbour.py
@@ -1,11 +1,15 @@
-from pybambi.neuralnetworks.nearestneighbor import NearestNeighborInterpolation
 import numpy
+from pybambi.neuralnetworks.nearestneighbour \
+        import NearestNeighbourInterpolation
 
 
-def test_NearestNeighborInterpolation():
+def test_NearestNeighbourInterpolation():
+    numpy.random.seed(0)
     params = numpy.array([[0, 0], [0, 1], [1, 1], [1, 0]])
+    params = numpy.repeat(params, 10, axis=0)
     logL = numpy.array([1, 2, 3, 4])
-    p = NearestNeighborInterpolation(params, logL)
+    logL = numpy.repeat(logL, 10)
+    p = NearestNeighbourInterpolation(params, logL)
     for t, l in zip(params, logL):
         assert p(t) == l