diff --git a/CHANGELOG.md b/CHANGELOG.md
index 46e2f1831f9ea38f7340440a7b0cabec28957749..0dba60a7e1d94b6f93bfc100099916541d209e83 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,7 @@ Changes currently on master, but not under a tag.
 - Added method to result to get injection recovery credible levels
 - Added function to generate a pp-plot from many results to core/result.py
 - Fixed a bug which caused `Interferometer.detector_tensor` not to update when `latitude`, `longitude`, `xarm_azimuth`, `yarm_azimuth`, `xarm_tilt`, `yarm_tilt` were updated.
+- Extracted time and frequency series behaviour from `WaveformGenerator` and `InterferometerStrainData` and moved it to `series.gw.CoupledTimeAndFrequencySeries`
 
 ### Changes
 - Switch the ordering the key-word arguments in `result.read_in_result` to put
diff --git a/bilby/core/likelihood.py b/bilby/core/likelihood.py
index 3698852d262addcd1525343b8bbf876ae956f256..910f6eb3226c14ea4f05e2124a4a024c6411a781 100644
--- a/bilby/core/likelihood.py
+++ b/bilby/core/likelihood.py
@@ -1,8 +1,9 @@
 from __future__ import division, print_function
-
 import copy
+
 import numpy as np
 from scipy.special import gammaln
+
 from .utils import infer_parameters_from_function
 
 
diff --git a/bilby/core/prior.py b/bilby/core/prior.py
index 82b9c59d8ef1aa2bcc757836e84ffe706e95719b..fac02ac7928e3b422f3627cb8fa7c77b5922d620 100644
--- a/bilby/core/prior.py
+++ b/bilby/core/prior.py
@@ -1,17 +1,18 @@
 from __future__ import division
 
-import numpy as np
-from scipy.interpolate import interp1d
-from scipy.integrate import cumtrapz
-from scipy.special import erf, erfinv
-import scipy.stats
 import os
 from collections import OrderedDict
 from future.utils import iteritems
 
-from .utils import logger, infer_args_from_method
-from . import utils
+import numpy as np
+import scipy.stats
+from scipy.integrate import cumtrapz
+from scipy.interpolate import interp1d
+from scipy.special import erf, erfinv
+
+# Keep import bilby statement, it is necessary for some eval() statements
 import bilby  # noqa
+from .utils import logger, infer_args_from_method, check_directory_exists_and_if_not_mkdir
 
 
 class PriorDict(OrderedDict):
@@ -50,7 +51,7 @@ class PriorDict(OrderedDict):
             Output file naming scheme
         """
 
-        utils.check_directory_exists_and_if_not_mkdir(outdir)
+        check_directory_exists_and_if_not_mkdir(outdir)
         prior_file = os.path.join(outdir, "{}.prior".format(label))
         logger.debug("Writing priors to {}".format(prior_file))
         with open(prior_file, "w") as outfile:
diff --git a/bilby/core/result.py b/bilby/core/result.py
index c8e36a83a6333aa8033b261ea466cd07f4e15346..5a9963cf91dac47bf3b8db47653edd692bf1887a 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -2,13 +2,15 @@ from __future__ import division
 
 import os
 from distutils.version import LooseVersion
+from collections import OrderedDict, namedtuple
+
 import numpy as np
 import deepdish
 import pandas as pd
 import corner
+import scipy.stats
 import matplotlib
 import matplotlib.pyplot as plt
-from collections import OrderedDict, namedtuple
 
 from . import utils
 from .utils import (logger, infer_parameters_from_function,
@@ -927,6 +929,44 @@ class Result(object):
                     return np.all(A == B)
         return False
 
+    @property
+    def kde(self):
+        """ Kernel density estimate built from the stored posterior
+
+        Uses `scipy.stats.gaussian_kde` to generate the kernel density
+        """
+        try:
+            return self._kde
+        except AttributeError:
+            self._kde = scipy.stats.gaussian_kde(
+                self.posterior[self.search_parameter_keys].values.T)
+            return self._kde
+
+    def posterior_probability(self, sample):
+        """ Calculate the posterior probabily for a new sample
+
+        This queries a Kernel Density Estimate of the posterior to calculate
+        the posterior probability density for the new sample.
+
+        Parameters
+        ----------
+        sample: dict, or list of dictionaries
+            A dictionary containing all the keys from
+            self.search_parameter_keys and corresponding values at which to
+            calculate the posterior probability
+
+        Returns
+        -------
+        p: array-like,
+            The posterior probability of the sample
+
+        """
+        if isinstance(sample, dict):
+            sample = [sample]
+        ordered_sample = [[s[key] for key in self.search_parameter_keys]
+                          for s in sample]
+        return self.kde(ordered_sample)
+
 
 def plot_multiple(results, filename=None, labels=None, colours=None,
                   save=True, evidences=False, **kwargs):
diff --git a/bilby/core/sampler/base_sampler.py b/bilby/core/sampler/base_sampler.py
index b088610732aaa2e9f08d86011a0f9cbed1dbb695..f47308b4f8ad20a4c2305e27d019deb015c78cbb 100644
--- a/bilby/core/sampler/base_sampler.py
+++ b/bilby/core/sampler/base_sampler.py
@@ -1,7 +1,9 @@
 from __future__ import absolute_import
 import datetime
 import numpy as np
+
 from pandas import DataFrame
+
 from ..utils import logger, command_line_args
 from ..prior import Prior, PriorDict
 from ..result import Result, read_in_result
@@ -143,9 +145,9 @@ class Sampler(object):
         external_sampler_name = self.__class__.__name__.lower()
         try:
             self.external_sampler = __import__(external_sampler_name)
-        except ImportError:
-            raise ImportError(
-                "Sampler {} not installed on this system".format(external_sampler_name))
+        except (ImportError, SystemExit, ModuleNotFoundError):
+            raise SamplerNotInstalledError(
+                "Sampler {} is not installed on this system".format(external_sampler_name))
 
     def _verify_kwargs_against_default_kwargs(self):
         """
@@ -473,3 +475,15 @@ class MCMCSampler(Sampler):
         except emcee.autocorr.AutocorrError as e:
             self.result.max_autocorrelation_time = None
             logger.info("Unable to calculate autocorr time: {}".format(e))
+
+
+class Error(Exception):
+    """ Base class for all exceptions raised by this module """
+
+
+class SamplerError(Error):
+    """ Base class for Error related to samplers in this module """
+
+
+class SamplerNotInstalledError(SamplerError):
+    """ Base class for Error raised by not installed samplers """
diff --git a/bilby/core/sampler/cpnest.py b/bilby/core/sampler/cpnest.py
index 6642393030741fb873085e59784d3e6ff710c457..951bff7883a607a1e59c23c3d8e6b75d186cedc5 100644
--- a/bilby/core/sampler/cpnest.py
+++ b/bilby/core/sampler/cpnest.py
@@ -1,8 +1,10 @@
 from __future__ import absolute_import
+
 import numpy as np
 from pandas import DataFrame
-from ..utils import logger, check_directory_exists_and_if_not_mkdir
+
 from .base_sampler import NestedSampler
+from ..utils import logger, check_directory_exists_and_if_not_mkdir
 
 
 class Cpnest(NestedSampler):
diff --git a/bilby/core/sampler/dynesty.py b/bilby/core/sampler/dynesty.py
index 1175ab62165f3d863f8b6516912652b88bd80747..3e5209c6e1be19c68b3985f57c8834017314d572 100644
--- a/bilby/core/sampler/dynesty.py
+++ b/bilby/core/sampler/dynesty.py
@@ -2,9 +2,11 @@ from __future__ import absolute_import
 
 import os
 import sys
+
 import numpy as np
 from pandas import DataFrame
 from deepdish.io import load, save
+
 from ..utils import logger, check_directory_exists_and_if_not_mkdir
 from .base_sampler import Sampler, NestedSampler
 
diff --git a/bilby/core/sampler/emcee.py b/bilby/core/sampler/emcee.py
index 5385fe4410a13db36661e2d5d337f37b46d3cbeb..932462884755676ff52d721ef458b83f27c6e5cf 100644
--- a/bilby/core/sampler/emcee.py
+++ b/bilby/core/sampler/emcee.py
@@ -1,6 +1,8 @@
 from __future__ import absolute_import, print_function
+
 import numpy as np
 from pandas import DataFrame
+
 from ..utils import logger, get_progress_bar
 from .base_sampler import MCMCSampler
 
diff --git a/bilby/core/sampler/nestle.py b/bilby/core/sampler/nestle.py
index ffcc14e0f62935cbd9617b4fefa91b820082eddf..652dc0d3205115ef6040203f94671a11ea9f6416 100644
--- a/bilby/core/sampler/nestle.py
+++ b/bilby/core/sampler/nestle.py
@@ -2,6 +2,7 @@ from __future__ import absolute_import
 
 import numpy as np
 from pandas import DataFrame
+
 from .base_sampler import NestedSampler
 
 
diff --git a/bilby/core/sampler/ptemcee.py b/bilby/core/sampler/ptemcee.py
index f73afdb53de2868ce74d565000f99631e7c00715..ca252471061a4ca03d9291e05dd1734f58d56206 100644
--- a/bilby/core/sampler/ptemcee.py
+++ b/bilby/core/sampler/ptemcee.py
@@ -1,6 +1,7 @@
 from __future__ import absolute_import, division, print_function
 
 import numpy as np
+
 from ..utils import get_progress_bar, logger
 from . import Emcee
 
@@ -57,7 +58,6 @@ class Ptemcee(Emcee):
     def run_sampler(self):
         import ptemcee
         tqdm = get_progress_bar()
-
         sampler = ptemcee.Sampler(dim=self.ndim, logl=self.log_likelihood,
                                   logp=self.log_prior, **self.sampler_init_kwargs)
         self.pos0 = [[self.get_random_draw_from_prior()
diff --git a/bilby/core/sampler/pymc3.py b/bilby/core/sampler/pymc3.py
index c35d5bd8becd5a696365b3b990dc7b01b141af29..bdfecffa86ade96880f3265b23813bfca9cf8ce3 100644
--- a/bilby/core/sampler/pymc3.py
+++ b/bilby/core/sampler/pymc3.py
@@ -1,12 +1,16 @@
 from __future__ import absolute_import, print_function
 
 from collections import OrderedDict
+
 import numpy as np
 
 from ..utils import derivatives, logger, infer_args_from_method
-from ..prior import Prior
+from ..prior import Prior, DeltaFunction, Sine, Cosine, PowerLaw
 from ..result import Result
 from .base_sampler import Sampler, MCMCSampler
+from ..likelihood import GaussianLikelihood, PoissonLikelihood, ExponentialLikelihood, \
+    StudentTLikelihood
+from ...gw.likelihood import BasicGravitationalWaveTransient, GravitationalWaveTransient
 
 
 class Pymc3(MCMCSampler):
@@ -64,6 +68,20 @@ class Pymc3(MCMCSampler):
         self.draws = draws
         self.chains = self.__kwargs['chains']
 
+    @staticmethod
+    def _import_external_sampler():
+        import pymc3
+        from pymc3.sampling import STEP_METHODS
+        from pymc3.theanof import floatX
+        return pymc3, STEP_METHODS, floatX
+
+    @staticmethod
+    def _import_theano():
+        import theano  # noqa
+        import theano.tensor as tt
+        from theano.compile.ops import as_op  # noqa
+        return theano, tt, as_op
+
     def _verify_parameters(self):
         """
         Change `_verify_parameters()` to just pass, i.e., don't try and
@@ -217,8 +235,6 @@ class Pymc3(MCMCSampler):
         Map the bilby delta function prior to a single value for PyMC3.
         """
 
-        from ..prior import DeltaFunction
-
         # check prior is a DeltaFunction
         if isinstance(self.priors[key], DeltaFunction):
             return self.priors[key].peak
@@ -230,17 +246,10 @@ class Pymc3(MCMCSampler):
         Map the bilby Sine prior to a PyMC3 style function
         """
 
-        from ..prior import Sine
-
         # check prior is a Sine
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
         if isinstance(self.priors[key], Sine):
-            import pymc3
-
-            try:
-                import theano.tensor as tt
-                from pymc3.theanof import floatX
-            except ImportError:
-                raise ImportError("You must have Theano installed to use PyMC3")
 
             class Pymc3Sine(pymc3.Continuous):
                 def __init__(self, lower=0., upper=np.pi):
@@ -277,18 +286,10 @@ class Pymc3(MCMCSampler):
         Map the bilby Cosine prior to a PyMC3 style function
         """
 
-        from ..prior import Cosine
-
         # check prior is a Cosine
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
         if isinstance(self.priors[key], Cosine):
-            import pymc3
-
-            # import theano
-            try:
-                import theano.tensor as tt
-                from pymc3.theanof import floatX
-            except ImportError:
-                raise ImportError("You must have Theano installed to use PyMC3")
 
             class Pymc3Cosine(pymc3.Continuous):
                 def __init__(self, lower=-np.pi / 2., upper=np.pi / 2.):
@@ -324,23 +325,15 @@ class Pymc3(MCMCSampler):
         Map the bilby PowerLaw prior to a PyMC3 style function
         """
 
-        from ..prior import PowerLaw
-
         # check prior is a PowerLaw
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
         if isinstance(self.priors[key], PowerLaw):
-            import pymc3
 
             # check power law is set
             if not hasattr(self.priors[key], 'alpha'):
                 raise AttributeError("No 'alpha' attribute set for PowerLaw prior")
 
-            # import theano
-            try:
-                import theano.tensor as tt
-                from pymc3.theanof import floatX
-            except ImportError:
-                raise ImportError("You must have Theano installed to use PyMC3")
-
             if self.priors[key].alpha < -1.:
                 # use Pareto distribution
                 palpha = -(1. + self.priors[key].alpha)
@@ -385,11 +378,8 @@ class Pymc3(MCMCSampler):
             raise ValueError("Prior for '{}' is not a Power Law".format(key))
 
     def run_sampler(self):
-        import pymc3
         # set the step method
-
-        from pymc3.sampling import STEP_METHODS
-
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
         step_methods = {m.__name__.lower(): m.__name__ for m in STEP_METHODS}
         if 'step' in self.__kwargs:
             self.step_method = self.__kwargs.pop('step')
@@ -470,8 +460,7 @@ class Pymc3(MCMCSampler):
         self.setup_prior_mapping()
 
         self.pymc3_priors = OrderedDict()
-
-        import pymc3
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
 
         # set the parameter prior distributions (in the model context manager)
         with self.pymc3_model:
@@ -534,18 +523,10 @@ class Pymc3(MCMCSampler):
         Convert any bilby likelihoods to PyMC3 distributions.
         """
 
-        try:
-            import theano  # noqa
-            import theano.tensor as tt
-            from theano.compile.ops import as_op  # noqa
-        except ImportError:
-            raise ImportError("Could not import theano")
-
-        from ..likelihood import GaussianLikelihood, PoissonLikelihood, ExponentialLikelihood, \
-            StudentTLikelihood
-        from ...gw.likelihood import BasicGravitationalWaveTransient, GravitationalWaveTransient
-
         # create theano Op for the log likelihood if not using a predefined model
+        pymc3, STEP_METHODS, floatX = self._import_external_sampler()
+        theano, tt, as_op = self._import_theano()
+
         class LogLike(tt.Op):
 
             itypes = [tt.dvector]
@@ -605,8 +586,6 @@ class Pymc3(MCMCSampler):
 
                 outputs[0][0] = grads
 
-        import pymc3
-
         with self.pymc3_model:
             #  check if it is a predefined likelhood function
             if isinstance(self.likelihood, GaussianLikelihood):
diff --git a/bilby/core/sampler/pymultinest.py b/bilby/core/sampler/pymultinest.py
index 0bfdb0eb96b33a5193ef59ddf06eaf60361c0aa8..d566740ad4d70e2463eaed256123733c63c3bf73 100644
--- a/bilby/core/sampler/pymultinest.py
+++ b/bilby/core/sampler/pymultinest.py
@@ -1,8 +1,8 @@
 from __future__ import absolute_import
 
 from ..utils import check_directory_exists_and_if_not_mkdir
-from .base_sampler import NestedSampler
 from ..utils import logger
+from .base_sampler import NestedSampler
 
 
 class Pymultinest(NestedSampler):
@@ -56,7 +56,7 @@ class Pymultinest(NestedSampler):
         https://github.com/JohannesBuchner/PyMultiNest/issues/115
         """
         if not self.kwargs['outputfiles_basename']:
-            self.kwargs['outputfiles_basename'] =\
+            self.kwargs['outputfiles_basename'] = \
                 '{}/pm_{}/'.format(self.outdir, self.label)
         if self.kwargs['outputfiles_basename'].endswith('/') is False:
             self.kwargs['outputfiles_basename'] = '{}/'.format(
diff --git a/bilby/core/utils.py b/bilby/core/utils.py
index 6d68f35275791df90c7e94a74adc49a122eb9b70..56eafe0ec8741823a88a0e80554ec36741c82c88 100644
--- a/bilby/core/utils.py
+++ b/bilby/core/utils.py
@@ -1,13 +1,15 @@
 from __future__ import division
+
 import logging
 import os
-import numpy as np
 from math import fmod
 import argparse
 import traceback
 import inspect
 import subprocess
 
+import numpy as np
+
 logger = logging.getLogger('bilby')
 
 # Constants
diff --git a/bilby/gw/calibration.py b/bilby/gw/calibration.py
index 1ed90538a15d683d698ed6358197d16ffedfd067..d95a12f002a391188ec5efcc3a9dbc44594b0e39 100644
--- a/bilby/gw/calibration.py
+++ b/bilby/gw/calibration.py
@@ -49,6 +49,9 @@ class Recalibrate(object):
         self.params.update({key[len(self.prefix):]: params[key] for key in params
                             if self.prefix in key})
 
+    def __eq__(self, other):
+        return self.__dict__ == other.__dict__
+
 
 class CubicSpline(Recalibrate):
 
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index d33226ef511a732d4ec7441858c4481c02057ad3..bc1d0ec417f4939b1f17391341ad2cee3d41d7a5 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -1,4 +1,5 @@
 from __future__ import division
+
 import numpy as np
 from pandas import DataFrame
 
diff --git a/bilby/gw/detector.py b/bilby/gw/detector.py
index a35a7b64e8455fa16fd6ee43162622a0db497485..22845a6d283bdb3a2c34bbb295ee16058f3f0c26 100644
--- a/bilby/gw/detector.py
+++ b/bilby/gw/detector.py
@@ -1,11 +1,13 @@
 from __future__ import division, print_function, absolute_import
 
 import os
+import sys
 
 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.signal.windows import tukey
 from scipy.interpolate import interp1d
+import deepdish as dd
 
 from . import utils as gwutils
 from ..core import utils
@@ -38,12 +40,12 @@ class InterferometerList(list):
 
         list.__init__(self)
         if type(interferometers) == str:
-            raise ValueError("Input must not be a string")
+            raise TypeError("Input must not be a string")
         for ifo in interferometers:
             if type(ifo) == str:
                 ifo = get_empty_interferometer(ifo)
             if type(ifo) not in [Interferometer, TriangularInterferometer]:
-                raise ValueError("Input list of interferometers are not all Interferometer objects")
+                raise TypeError("Input list of interferometers are not all Interferometer objects")
             else:
                 self.append(ifo)
         self._check_interferometers()
@@ -108,7 +110,7 @@ class InterferometerList(list):
         """
         if injection_polarizations is None:
             if waveform_generator is not None:
-                injection_polarizations =\
+                injection_polarizations = \
                     waveform_generator.frequency_domain_strain(parameters)
             else:
                 raise ValueError(
@@ -213,6 +215,48 @@ class InterferometerList(list):
         return {interferometer.name: interferometer.meta_data
                 for interferometer in self}
 
+    @staticmethod
+    def _hdf5_filename_from_outdir_label(outdir, label):
+        return os.path.join(outdir, label + '.h5')
+
+    def to_hdf5(self, outdir='outdir', label='ifo_list'):
+        """ Saves the object to a hdf5 file
+
+        Parameters
+        ----------
+        outdir: str, optional
+            Output directory name of the file
+        label: str, optional
+            Output file name, is 'ifo_list' if not given otherwise. A list of
+            the included interferometers will be appended.
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of InterferometerList is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+        label = label + '_' + ''.join(ifo.name for ifo in self)
+        utils.check_directory_exists_and_if_not_mkdir(outdir)
+        dd.io.save(self._hdf5_filename_from_outdir_label(outdir, label), self)
+
+    @classmethod
+    def from_hdf5(cls, filename=None):
+        """ Loads in an InterferometerList object from an hdf5 file
+
+        Parameters
+        ----------
+        filename: str
+            If given, try to load from this filename
+
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of InterferometerList is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+        res = dd.io.load(filename)
+        if res.__class__ == list:
+            res = cls(res)
+        if res.__class__ != cls:
+            raise TypeError('The loaded object is not a InterferometerList')
+        return res
+
 
 class InterferometerStrainData(object):
     """ Strain data for an interferometer """
@@ -248,6 +292,21 @@ class InterferometerStrainData(object):
         self._time_domain_strain = None
         self._time_array = None
 
+    def __eq__(self, other):
+        if self.minimum_frequency == other.minimum_frequency \
+                and self.maximum_frequency == other.maximum_frequency \
+                and self.roll_off == other.roll_off \
+                and self.window_factor == other.window_factor \
+                and self.sampling_frequency == other.sampling_frequency \
+                and self.duration == other.duration \
+                and self.start_time == other.start_time \
+                and np.array_equal(self.time_array, other.time_array) \
+                and np.array_equal(self.frequency_array, other.frequency_array) \
+                and np.array_equal(self.frequency_domain_strain, other.frequency_domain_strain) \
+                and np.array_equal(self.time_domain_strain, other.time_domain_strain):
+            return True
+        return False
+
     def time_within_data(self, time):
         """ Check if time is within the data span
 
@@ -815,6 +874,22 @@ class Interferometer(object):
             maximum_frequency=maximum_frequency)
         self.meta_data = dict()
 
+    def __eq__(self, other):
+        if self.name == other.name and \
+                self.length == other.length and \
+                self.latitude == other.latitude and \
+                self.longitude == other.longitude and \
+                self.elevation == other.elevation and \
+                self.xarm_azimuth == other.xarm_azimuth and \
+                self.xarm_tilt == other.xarm_tilt and \
+                self.yarm_azimuth == other.yarm_azimuth and \
+                self.yarm_tilt == other.yarm_tilt and \
+                self.power_spectral_density.__eq__(other.power_spectral_density) and \
+                self.calibration_model == other.calibration_model and \
+                self.strain_data == other.strain_data:
+            return True
+        return False
+
     def __repr__(self):
         return self.__class__.__name__ + '(name=\'{}\', power_spectral_density={}, minimum_frequency={}, ' \
                                          'maximum_frequency={}, length={}, latitude={}, longitude={}, elevation={}, ' \
@@ -1249,7 +1324,7 @@ class Interferometer(object):
 
         if injection_polarizations is None:
             if waveform_generator is not None:
-                injection_polarizations =\
+                injection_polarizations = \
                     waveform_generator.frequency_domain_strain(parameters)
             else:
                 raise ValueError(
@@ -1580,6 +1655,48 @@ class Interferometer(object):
             fig.savefig(
                 '{}/{}_{}_time_domain_data.png'.format(outdir, self.name, label))
 
+    @staticmethod
+    def _hdf5_filename_from_outdir_label(outdir, label):
+        return os.path.join(outdir, label + '.h5')
+
+    def to_hdf5(self, outdir='outdir', label=None):
+        """ Save the object to a hdf5 file
+
+        Attributes
+        ----------
+        outdir: str, optional
+            Output directory name of the file, defaults to 'outdir'.
+        label: str, optional
+            Output file name, is self.name if not given otherwise.
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of Interferometer is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+        if label is None:
+            label = self.name
+        utils.check_directory_exists_and_if_not_mkdir('outdir')
+        filename = self._hdf5_filename_from_outdir_label(outdir, label)
+        dd.io.save(filename, self)
+
+    @classmethod
+    def from_hdf5(cls, filename=None):
+        """ Loads in an Interferometer object from an hdf5 file
+
+        Parameters
+        ----------
+        filename: str
+            If given, try to load from this filename
+
+        """
+        if sys.version_info[0] < 3:
+            raise NotImplementedError('Pickling of Interferometer is not supported in Python 2.'
+                                      'Use Python 3 instead.')
+
+        res = dd.io.load(filename)
+        if res.__class__ != cls:
+            raise TypeError('The loaded object is not an Interferometer')
+        return res
+
 
 class TriangularInterferometer(InterferometerList):
 
@@ -1648,6 +1765,15 @@ class PowerSpectralDensity(object):
         self.psd_file = psd_file
         self.asd_file = asd_file
 
+    def __eq__(self, other):
+        if self.psd_file == other.psd_file \
+                and self.asd_file == other.asd_file \
+                and np.array_equal(self.frequency_array, other.frequency_array) \
+                and np.array_equal(self.psd_array, other.psd_array) \
+                and np.array_equal(self.asd_array, other.asd_array):
+            return True
+        return False
+
     def __repr__(self):
         if self.asd_file is not None or self.psd_file is not None:
             return self.__class__.__name__ + '(psd_file=\'{}\', asd_file=\'{}\')' \
@@ -1769,12 +1895,12 @@ class PowerSpectralDensity(object):
 
     @property
     def asd_file(self):
-        return self.__asd_file
+        return self._asd_file
 
     @asd_file.setter
     def asd_file(self, asd_file):
         asd_file = self.__validate_file_name(file=asd_file)
-        self.__asd_file = asd_file
+        self._asd_file = asd_file
         if asd_file is not None:
             self.__import_amplitude_spectral_density()
             self.__check_file_was_asd_file()
@@ -1788,12 +1914,12 @@ class PowerSpectralDensity(object):
 
     @property
     def psd_file(self):
-        return self.__psd_file
+        return self._psd_file
 
     @psd_file.setter
     def psd_file(self, psd_file):
         psd_file = self.__validate_file_name(file=psd_file)
-        self.__psd_file = psd_file
+        self._psd_file = psd_file
         if psd_file is not None:
             self.__import_power_spectral_density()
             self.__check_file_was_psd_file()
@@ -2155,16 +2281,16 @@ def load_data_from_cache_file(
             frame_duration = float(frame_duration)
             if frame_name[:4] == 'file':
                 frame_name = frame_name[16:]
-            if not data_set & (frame_start < segment_start) &\
+            if not data_set & (frame_start < segment_start) & \
                     (segment_start < frame_start + frame_duration):
                 ifo.set_strain_data_from_frame_file(
                     frame_name, 4096, segment_duration,
                     start_time=segment_start,
                     channel=channel_name, buffer_time=0)
                 data_set = True
-            if not psd_set & (frame_start < psd_start) &\
+            if not psd_set & (frame_start < psd_start) & \
                     (psd_start + psd_duration < frame_start + frame_duration):
-                ifo.power_spectral_density =\
+                ifo.power_spectral_density = \
                     PowerSpectralDensity.from_frame_file(
                         frame_name, psd_start_time=psd_start,
                         psd_duration=psd_duration,
diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index eafa9da755b7beb497a6e243df902727d776e9ff..b56a377c603764b4e7853e9a34e438f39161f037 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -285,7 +285,7 @@ class GravitationalWaveTransient(likelihood.Likelihood):
             self.interferometers.start_time + np.linspace(
                 0, self.interferometers.duration,
                 int(self.interferometers.duration / 2 *
-                    self.waveform_generator.sampling_frequency))[1:]
+                    self.waveform_generator.sampling_frequency + 1))[1:]
         self.time_prior_array =\
             self.priors['geocent_time'].prob(times) * delta_tc
 
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 46f62708c94e5b57488943e6bdfe37a418df52a0..d9eb726240e5d421e35290b7ad056c5334954f91 100644
--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
@@ -1,6 +1,8 @@
 import os
+
 import numpy as np
 from scipy.interpolate import UnivariateSpline
+
 from ..core.prior import (PriorDict, Uniform, FromFile, Prior, DeltaFunction,
                           Gaussian, Interped)
 from ..core.utils import logger
diff --git a/bilby/gw/series.py b/bilby/gw/series.py
index aab94dee94d96db575af6ae6ee42c7b719250f7b..a83a0a658bce8dcc1f2351f9dcfe0f25a82f4c01 100644
--- a/bilby/gw/series.py
+++ b/bilby/gw/series.py
@@ -20,6 +20,12 @@ class CoupledTimeAndFrequencySeries(object):
         self.start_time = start_time
         self._frequency_array_updated = False
         self._time_array_updated = False
+        self._frequency_array = None
+        self._time_array = None
+
+    def __repr__(self):
+        return self.__class__.__name__ + '(duration={}, sampling_frequency={}, start_time={})'\
+            .format(self.duration, self.sampling_frequency, self.start_time)
 
     @property
     def frequency_array(self):
@@ -31,7 +37,7 @@ class CoupledTimeAndFrequencySeries(object):
         """
         if self._frequency_array_updated is False:
             if self.sampling_frequency and self.duration:
-                self.frequency_array = utils.create_frequency_series(
+                self._frequency_array = utils.create_frequency_series(
                     sampling_frequency=self.sampling_frequency,
                     duration=self.duration)
             else:
diff --git a/bilby/gw/waveform_generator.py b/bilby/gw/waveform_generator.py
index 6e9135731291f0f80172fea753583d5d5d22c2f0..6e9a5f6a79a698ec963062529b9bb02177df4874 100644
--- a/bilby/gw/waveform_generator.py
+++ b/bilby/gw/waveform_generator.py
@@ -1,4 +1,5 @@
 import numpy as np
+
 from ..core import utils
 from ..gw.series import CoupledTimeAndFrequencySeries
 
@@ -167,7 +168,7 @@ class WaveformGenerator(object):
         model_strain = dict()
         for key in transformed_model_strain:
             if transformation_function == utils.nfft:
-                model_strain[key], self.frequency_array = \
+                model_strain[key], _ = \
                     transformation_function(transformed_model_strain[key], self.sampling_frequency)
             else:
                 model_strain[key] = transformation_function(transformed_model_strain[key], self.sampling_frequency)
diff --git a/bilby/hyper/likelihood.py b/bilby/hyper/likelihood.py
index 200b2e69fd573d24a1006ae2582daeac32e10a4f..3325bca17c82e5d2629f52332db9b474e4d2025d 100644
--- a/bilby/hyper/likelihood.py
+++ b/bilby/hyper/likelihood.py
@@ -1,7 +1,9 @@
 from __future__ import division, print_function
 
 import logging
+
 import numpy as np
+
 from ..core.likelihood import Likelihood
 from .model import Model
 
diff --git a/bilby/hyper/model.py b/bilby/hyper/model.py
index 6fe300cf9a83d427fbce6a5993a2c07060460a03..5eb1fe9d9d73aad72aa5f3cc247ad17116fb7433 100644
--- a/bilby/hyper/model.py
+++ b/bilby/hyper/model.py
@@ -24,11 +24,9 @@ class Model(object):
                 self.parameters[key] = None
 
     def prob(self, data):
+        probability = 1.0
         for ii, function in enumerate(self.models):
-            if ii == 0:
-                probability = function(data, **self._get_function_parameters(function))
-            else:
-                probability *= function(data, **self._get_function_parameters(function))
+            probability *= function(data, **self._get_function_parameters(function))
         return probability
 
     def _get_function_parameters(self, func):
diff --git a/test/detector_test.py b/test/detector_test.py
index 7ca202e6de1305c3e54e0194343cf06e5472cdd5..028e1d33f6e16eb69c5965663a6a2e817a619fc7 100644
--- a/test/detector_test.py
+++ b/test/detector_test.py
@@ -7,16 +7,18 @@ from mock import MagicMock
 from mock import patch
 import numpy as np
 import scipy.signal.windows
-import gwpy
 import os
+import sys
+from shutil import rmtree
 import logging
+import deepdish as dd
 
 
-class TestDetector(unittest.TestCase):
+class TestInterferometer(unittest.TestCase):
 
     def setUp(self):
         self.name = 'name'
-        self.power_spectral_density = MagicMock()
+        self.power_spectral_density = bilby.gw.detector.PowerSpectralDensity.from_aligo()
         self.minimum_frequency = 10
         self.maximum_frequency = 20
         self.length = 30
@@ -37,6 +39,7 @@ class TestDetector(unittest.TestCase):
                                                     xarm_tilt=self.xarm_tilt, yarm_tilt=self.yarm_tilt)
         self.ifo.strain_data.set_from_frequency_domain_strain(
             np.linspace(0, 4096, 4097), sampling_frequency=4096, duration=2)
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir')
 
     def tearDown(self):
         del self.name
@@ -52,6 +55,7 @@ class TestDetector(unittest.TestCase):
         del self.xarm_tilt
         del self.yarm_tilt
         del self.ifo
+        rmtree('outdir')
 
     def test_name_setting(self):
         self.assertEqual(self.ifo.name, self.name)
@@ -323,7 +327,10 @@ class TestDetector(unittest.TestCase):
             signal = 1
             expected = [signal, signal, self.ifo.power_spectral_density_array, self.ifo.strain_data.duration]
             actual = self.ifo.optimal_snr_squared(signal=signal)
-            self.assertListEqual(expected, actual)
+            self.assertEqual(expected[0], actual[0])
+            self.assertEqual(expected[1], actual[1])
+            self.assertTrue(np.array_equal(expected[2], actual[2]))
+            self.assertEqual(expected[3], actual[3])
 
     def test_matched_filter_snr_squared(self):
         """ Merely checks parameters are given in the right order """
@@ -334,7 +341,9 @@ class TestDetector(unittest.TestCase):
                                                            self.ifo.strain_data.duration]]
             actual = self.ifo.matched_filter_snr_squared(signal=signal)
             self.assertTrue(np.array_equal(expected[0], actual[0]))  # array-like element has to be evaluated separately
-            self.assertListEqual(expected[1], actual[1])
+            self.assertEqual(expected[1][0], actual[1][0])
+            self.assertTrue(np.array_equal(expected[1][1], actual[1][1]))
+            self.assertEqual(expected[1][2], actual[1][2])
 
     def test_repr(self):
         expected = 'Interferometer(name=\'{}\', power_spectral_density={}, minimum_frequency={}, ' \
@@ -346,6 +355,145 @@ class TestDetector(unittest.TestCase):
                     float(self.yarm_tilt))
         self.assertEqual(expected, repr(self.ifo))
 
+    def test_to_and_from_hdf5_loading(self):
+        if sys.version_info[0] < 3:
+            with self.assertRaises(NotImplementedError):
+                self.ifo.to_hdf5(outdir='outdir', label='test')
+        else:
+            self.ifo.to_hdf5(outdir='outdir', label='test')
+            filename = self.ifo._hdf5_filename_from_outdir_label(outdir='outdir', label='test')
+            recovered_ifo = bilby.gw.detector.Interferometer.from_hdf5(filename)
+            self.assertEqual(self.ifo, recovered_ifo)
+
+    def test_to_and_from_hdf5_wrong_class(self):
+        if sys.version_info[0] < 3:
+            pass
+        else:
+            bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir')
+            dd.io.save('./outdir/psd.h5', self.power_spectral_density)
+            filename = self.ifo._hdf5_filename_from_outdir_label(outdir='outdir', label='psd')
+            with self.assertRaises(TypeError):
+                bilby.gw.detector.Interferometer.from_hdf5(filename)
+
+
+class TestInterferometerEquals(unittest.TestCase):
+
+    def setUp(self):
+        self.name = 'name'
+        self.power_spectral_density_1 = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        self.power_spectral_density_2 = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        self.minimum_frequency = 10
+        self.maximum_frequency = 20
+        self.length = 30
+        self.latitude = 1
+        self.longitude = 2
+        self.elevation = 3
+        self.xarm_azimuth = 4
+        self.yarm_azimuth = 5
+        self.xarm_tilt = 0.
+        self.yarm_tilt = 0.
+        # noinspection PyTypeChecker
+        self.duration = 1
+        self.sampling_frequency = 200
+        self.frequency_array = bilby.utils.create_frequency_series(sampling_frequency=self.sampling_frequency,
+                                                                   duration=self.duration)
+        self.strain = self.frequency_array
+        self.ifo_1 = bilby.gw.detector.Interferometer(name=self.name,
+                                                      power_spectral_density=self.power_spectral_density_1,
+                                                      minimum_frequency=self.minimum_frequency,
+                                                      maximum_frequency=self.maximum_frequency, length=self.length,
+                                                      latitude=self.latitude, longitude=self.longitude,
+                                                      elevation=self.elevation,
+                                                      xarm_azimuth=self.xarm_azimuth, yarm_azimuth=self.yarm_azimuth,
+                                                      xarm_tilt=self.xarm_tilt, yarm_tilt=self.yarm_tilt)
+        self.ifo_2 = bilby.gw.detector.Interferometer(name=self.name,
+                                                      power_spectral_density=self.power_spectral_density_2,
+                                                      minimum_frequency=self.minimum_frequency,
+                                                      maximum_frequency=self.maximum_frequency, length=self.length,
+                                                      latitude=self.latitude, longitude=self.longitude,
+                                                      elevation=self.elevation,
+                                                      xarm_azimuth=self.xarm_azimuth, yarm_azimuth=self.yarm_azimuth,
+                                                      xarm_tilt=self.xarm_tilt, yarm_tilt=self.yarm_tilt)
+        self.ifo_1.set_strain_data_from_frequency_domain_strain(frequency_array=self.frequency_array,
+                                                                frequency_domain_strain=self.strain)
+        self.ifo_2.set_strain_data_from_frequency_domain_strain(frequency_array=self.frequency_array,
+                                                                frequency_domain_strain=self.strain)
+
+    def tearDown(self):
+        del self.name
+        del self.power_spectral_density_1
+        del self.power_spectral_density_2
+        del self.minimum_frequency
+        del self.maximum_frequency
+        del self.length
+        del self.latitude
+        del self.longitude
+        del self.elevation
+        del self.xarm_azimuth
+        del self.yarm_azimuth
+        del self.xarm_tilt
+        del self.yarm_tilt
+        del self.ifo_1
+        del self.ifo_2
+        del self.sampling_frequency
+        del self.duration
+        del self.frequency_array
+        del self.strain
+
+    def test_eq_true(self):
+        self.assertEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_psd(self):
+        self.ifo_1.power_spectral_density.psd_array[0] = 1234
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_minimum_frequency(self):
+        self.ifo_1.minimum_frequency -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_maximum_frequency(self):
+        self.ifo_1.minimum_frequency -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_length(self):
+        self.ifo_1.length -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_latitude(self):
+        self.ifo_1.latitude -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_longitude(self):
+        self.ifo_1.longitude -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_elevation(self):
+        self.ifo_1.elevation -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_xarm_azimuth(self):
+        self.ifo_1.xarm_azimuth -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_xarmtilt(self):
+        self.ifo_1.xarm_tilt -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_yarm_azimuth(self):
+        self.ifo_1.yarm_azimuth -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_yarm_tilt(self):
+        self.ifo_1.yarm_tilt -= 1
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
+    def test_eq_false_different_ifo_strain_data(self):
+        self.strain = bilby.utils.create_frequency_series(sampling_frequency=self.sampling_frequency/2,
+                                                          duration=self.duration*2)
+        self.ifo_1.set_strain_data_from_frequency_domain_strain(frequency_array=self.frequency_array,
+                                                                frequency_domain_strain=self.strain)
+        self.assertNotEqual(self.ifo_1, self.ifo_2)
+
 
 class TestInterferometerStrainData(unittest.TestCase):
 
@@ -370,42 +518,6 @@ class TestInterferometerStrainData(unittest.TestCase):
                 frequency_domain_strain=np.array([0, 1, 2]), frequency_array=np.array([5, 15, 25]))
             self.assertTrue(np.array_equal(self.ifosd.frequency_mask, [False, True, False]))
 
-    def test_frequency_array_setting_direct(self):
-        with mock.patch('bilby.core.utils.create_frequency_series') as m:
-            m.return_value = np.array([5, 15, 25])
-            self.ifosd.set_from_frequency_domain_strain(
-                frequency_domain_strain=np.array([0, 1, 2]), frequency_array=np.array([5, 15, 25]))
-            self.assertTrue(np.array_equal(self.ifosd.frequency_array, np.array(np.array([5, 15, 25]))))
-
-    def test_duration_setting(self):
-        with mock.patch('bilby.core.utils.create_frequency_series') as m:
-            m.return_value = np.array([0, 1, 2])
-            self.ifosd.set_from_frequency_domain_strain(
-                frequency_domain_strain=np.array([0, 1, 2]), frequency_array=np.array([0, 1, 2]))
-            self.assertAlmostEqual(self.ifosd.duration, 1)
-
-    def test_sampling_frequency_setting(self):
-        with mock.patch('bilby.core.utils.create_frequency_series') as n:
-            with mock.patch('bilby.core.utils.get_sampling_frequency_and_duration_from_frequency_array') as m:
-                m.return_value = 8, 456
-                n.return_value = np.array([1, 2, 3])
-                self.ifosd.set_from_frequency_domain_strain(
-                    frequency_domain_strain=np.array([0, 1, 2]), frequency_array=np.array([0, 1, 2]))
-                self.assertEqual(8, self.ifosd.sampling_frequency)
-
-    def test_frequency_array_setting(self):
-        duration = 3
-        sampling_frequency = 1
-        with mock.patch('bilby.core.utils.create_frequency_series') as m:
-            m.return_value = [1, 2, 3]
-            self.ifosd.set_from_frequency_domain_strain(
-                frequency_domain_strain=np.array([0, 1, 2]), duration=duration,
-                sampling_frequency=sampling_frequency)
-            self.assertTrue(np.array_equal(
-                self.ifosd.frequency_array,
-                bilby.core.utils.create_frequency_series(duration=duration,
-                                                         sampling_frequency=sampling_frequency)))
-
     def test_set_data_fails(self):
         with mock.patch('bilby.core.utils.create_frequency_series') as m:
             m.return_value = [1, 2, 3]
@@ -458,60 +570,6 @@ class TestInterferometerStrainData(unittest.TestCase):
         self.assertTrue(np.all(
             self.ifosd.frequency_domain_strain == frequency_domain_strain * self.ifosd.frequency_mask))
 
-    def test_time_array_when_set(self):
-        test_array = np.array([1, 2, 3])
-        self.ifosd.time_array = test_array
-        self.assertTrue(np.array_equal(test_array, self.ifosd.time_array))
-
-    def test_time_array_when_not_set(self):
-        with mock.patch('bilby.core.utils.create_time_series') as m:
-            self.ifosd.start_time = 3
-            self.ifosd.sampling_frequency = 1000
-            self.ifosd.duration = 5
-            m.return_value = 4
-            self.assertEqual(m.return_value, self.ifosd.time_array)
-            m.assert_called_with(sampling_frequency=1000,
-                                 duration=5,
-                                 starting_time=3)
-
-    def test_time_array_without_sampling_frequency(self):
-        self.ifosd.sampling_frequency = None
-        self.ifosd.duration = 4
-        with self.assertRaises(ValueError):
-            test = self.ifosd.time_array
-
-    def test_time_array_without_duration(self):
-        self.ifosd.sampling_frequency = 4096
-        self.ifosd.duration = None
-        with self.assertRaises(ValueError):
-            test = self.ifosd.time_array
-
-    def test_frequency_array_when_set(self):
-        test_array = np.array([1, 2, 3])
-        self.ifosd.frequency_array = test_array
-        self.assertTrue(np.array_equal(test_array, self.ifosd.frequency_array))
-
-    def test_frequency_array_when_not_set(self):
-        with mock.patch('bilby.core.utils.create_frequency_series') as m:
-            m.return_value = [1, 2, 3]
-            self.ifosd.sampling_frequency = 1000
-            self.ifosd.duration = 5
-            self.assertListEqual(m.return_value, self.ifosd.frequency_array)
-            m.assert_called_with(sampling_frequency=1000,
-                                 duration=5)
-
-    def test_frequency_array_without_sampling_frequency(self):
-        self.ifosd.sampling_frequency = None
-        self.ifosd.duration = 4
-        with self.assertRaises(ValueError):
-            test = self.ifosd.frequency_array
-
-    def test_frequency_array_without_duration(self):
-        self.ifosd.sampling_frequency = 4096
-        self.ifosd.duration = None
-        with self.assertRaises(ValueError):
-            test = self.ifosd.frequency_array
-
     def test_time_within_data_before(self):
         self.ifosd.start_time = 3
         self.ifosd.duration = 2
@@ -613,18 +671,103 @@ class TestInterferometerStrainData(unittest.TestCase):
             self.ifosd.frequency_domain_strain = np.array([1])
 
 
+class TestInterferometerStrainDataEquals(unittest.TestCase):
+
+    def setUp(self):
+        self.minimum_frequency = 10
+        self.maximum_frequency = 20
+        self.roll_off = 0.2
+        self.sampling_frequency = 100
+        self.duration = 2
+        self.frequency_array = bilby.utils.create_frequency_series(sampling_frequency=self.sampling_frequency,
+                                                                   duration=self.duration)
+        self.strain = self.frequency_array
+        self.ifosd_1 = bilby.gw.detector.InterferometerStrainData(minimum_frequency=self.minimum_frequency,
+                                                                  maximum_frequency=self.maximum_frequency,
+                                                                  roll_off=self.roll_off)
+        self.ifosd_2 = bilby.gw.detector.InterferometerStrainData(minimum_frequency=self.minimum_frequency,
+                                                                  maximum_frequency=self.maximum_frequency,
+                                                                  roll_off=self.roll_off)
+        self.ifosd_1.set_from_frequency_domain_strain(frequency_domain_strain=self.strain,
+                                                      frequency_array=self.frequency_array)
+        self.ifosd_2.set_from_frequency_domain_strain(frequency_domain_strain=self.strain,
+                                                      frequency_array=self.frequency_array)
+
+    def tearDown(self):
+        del self.minimum_frequency
+        del self.maximum_frequency
+        del self.roll_off
+        del self.sampling_frequency
+        del self.duration
+        del self.frequency_array
+        del self.strain
+        del self.ifosd_1
+        del self.ifosd_2
+
+    def test_eq_true(self):
+        self.assertEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_minimum_frequency(self):
+        self.ifosd_1.minimum_frequency -= 1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_maximum_frequency(self):
+        self.ifosd_1.maximum_frequency -= 1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_roll_off(self):
+        self.ifosd_1.roll_off -= 0.1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_window_factor(self):
+        self.ifosd_1.roll_off -= 0.1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_sampling_frequency(self):
+        self.ifosd_1.sampling_frequency -= 0.1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_sampling_duration(self):
+        self.ifosd_1.duration -= 0.1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_start_time(self):
+        self.ifosd_1.start_time -= 0.1
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_frequency_array(self):
+        new_frequency_array = bilby.utils.create_frequency_series(sampling_frequency=self.sampling_frequency/2,
+                                                                  duration=self.duration*2)
+        self.ifosd_1.frequency_array = new_frequency_array
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_frequency_domain_strain(self):
+        new_strain = bilby.utils.create_frequency_series(sampling_frequency=self.sampling_frequency/2,
+                                                         duration=self.duration*2)
+        self.ifosd_1._frequency_domain_strain = new_strain
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_time_array(self):
+        new_time_array = bilby.utils.create_time_series(sampling_frequency=self.sampling_frequency/2,
+                                                        duration=self.duration*2)
+        self.ifosd_1.time_array = new_time_array
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+    def test_eq_different_time_domain_strain(self):
+        new_strain = bilby.utils.create_time_series(sampling_frequency=self.sampling_frequency/2,
+                                                    duration=self.duration*2)
+        self.ifosd_1._time_domain_strain= new_strain
+        self.assertNotEqual(self.ifosd_1, self.ifosd_2)
+
+
 class TestInterferometerList(unittest.TestCase):
 
     def setUp(self):
         self.frequency_arrays = np.linspace(0, 4096, 4097)
         self.name1 = 'name1'
         self.name2 = 'name2'
-        self.power_spectral_density1 = MagicMock()
-        self.power_spectral_density1.get_noise_realisation = MagicMock(return_value=(self.frequency_arrays,
-                                                                                     self.frequency_arrays))
-        self.power_spectral_density2 = MagicMock()
-        self.power_spectral_density2.get_noise_realisation = MagicMock(return_value=(self.frequency_arrays,
-                                                                                     self.frequency_arrays))
+        self.power_spectral_density1 = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        self.power_spectral_density2 = bilby.gw.detector.PowerSpectralDensity.from_aligo()
         self.minimum_frequency1 = 10
         self.minimum_frequency2 = 10
         self.maximum_frequency1 = 20
@@ -667,6 +810,7 @@ class TestInterferometerList(unittest.TestCase):
         self.ifo2.strain_data.set_from_frequency_domain_strain(
             self.frequency_arrays, sampling_frequency=4096, duration=2)
         self.ifo_list = bilby.gw.detector.InterferometerList([self.ifo1, self.ifo2])
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir')
 
     def tearDown(self):
         del self.frequency_arrays
@@ -697,20 +841,21 @@ class TestInterferometerList(unittest.TestCase):
         del self.ifo1
         del self.ifo2
         del self.ifo_list
+        rmtree('outdir')
 
     def test_init_with_string(self):
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             bilby.gw.detector.InterferometerList("string")
 
     def test_init_with_string_list(self):
         """ Merely checks if this ends up in the right bracket """
         with mock.patch('bilby.gw.detector.get_empty_interferometer') as m:
-            m.side_effect = ValueError
-            with self.assertRaises(ValueError):
+            m.side_effect = TypeError
+            with self.assertRaises(TypeError):
                 bilby.gw.detector.InterferometerList(['string'])
 
     def test_init_with_other_object(self):
-        with self.assertRaises(ValueError):
+        with self.assertRaises(TypeError):
             bilby.gw.detector.InterferometerList([object()])
 
     def test_init_with_actual_ifos(self):
@@ -833,6 +978,26 @@ class TestInterferometerList(unittest.TestCase):
         names = [ifo.name for ifo in self.ifo_list]
         self.assertListEqual([self.ifo1.name, new_ifo.name, self.ifo2.name], names)
 
+    def test_to_and_from_hdf5_loading(self):
+        if sys.version_info[0] < 3:
+            with self.assertRaises(NotImplementedError):
+                self.ifo_list.to_hdf5(outdir='outdir', label='test')
+        else:
+            self.ifo_list.to_hdf5(outdir='outdir', label='test')
+            filename = 'outdir/test_name1name2.h5'
+            recovered_ifo = bilby.gw.detector.InterferometerList.from_hdf5(filename)
+            self.assertListEqual(self.ifo_list, recovered_ifo)
+
+    def test_to_and_from_hdf5_wrong_class(self):
+        if sys.version_info[0] < 3:
+            pass
+        else:
+            dd.io.save('./outdir/psd.h5', self.ifo_list[0].power_spectral_density)
+            filename = self.ifo_list._hdf5_filename_from_outdir_label(
+                outdir='outdir', label='psd')
+            with self.assertRaises(TypeError):
+                bilby.gw.detector.InterferometerList.from_hdf5(filename)
+
 
 class TestPowerSpectralDensityWithoutFiles(unittest.TestCase):
 
@@ -1030,5 +1195,61 @@ class TestPowerSpectralDensityWithFiles(unittest.TestCase):
         self.assertEqual(expected, repr(psd))
 
 
+class TestPowerSpectralDensityEquals(unittest.TestCase):
+
+    def setUp(self):
+        self.psd_from_file_1 = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        self.psd_from_file_2 = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        self.frequency_array = np.linspace(1, 100)
+        self.psd_array = np.linspace(1, 100)
+        self.psd_from_array_1 = bilby.gw.detector.PowerSpectralDensity. \
+            from_power_spectral_density_array(frequency_array=self.frequency_array, psd_array= self.psd_array)
+        self.psd_from_array_2 = bilby.gw.detector.PowerSpectralDensity. \
+            from_power_spectral_density_array(frequency_array=self.frequency_array, psd_array= self.psd_array)
+
+    def tearDown(self):
+        del self.psd_from_file_1
+        del self.psd_from_file_2
+        del self.frequency_array
+        del self.psd_array
+        del self.psd_from_array_1
+        del self.psd_from_array_2
+
+    def test_eq_true_from_array(self):
+        self.assertEqual(self.psd_from_array_1, self.psd_from_array_2)
+
+    def test_eq_true_from_file(self):
+        self.assertEqual(self.psd_from_file_1, self.psd_from_file_2)
+
+    def test_eq_false_different_psd_file_name(self):
+        self.psd_from_file_1._psd_file = 'some_other_name'
+        self.assertNotEqual(self.psd_from_file_1, self.psd_from_file_2)
+
+    def test_eq_false_different_asd_file_name(self):
+        self.psd_from_file_1._psd_file = None
+        self.psd_from_file_2._psd_file = None
+        self.psd_from_file_1._asd_file = 'some_name'
+        self.psd_from_file_2._asd_file = 'some_other_name'
+        self.assertNotEqual(self.psd_from_file_1, self.psd_from_file_2)
+
+    def test_eq_false_different_frequency_array(self):
+        self.psd_from_file_1.frequency_array[0] = 0.5
+        self.psd_from_array_1.frequency_array[0] = 0.5
+        self.assertNotEqual(self.psd_from_file_1, self.psd_from_file_2)
+        self.assertNotEqual(self.psd_from_array_1, self.psd_from_array_2)
+
+    def test_eq_false_different_psd(self):
+        self.psd_from_file_1.psd_array[0] = 0.53544321
+        self.psd_from_array_1.psd_array[0] = 0.53544321
+        self.assertNotEqual(self.psd_from_file_1, self.psd_from_file_2)
+        self.assertNotEqual(self.psd_from_array_1, self.psd_from_array_2)
+
+    def test_eq_false_different_asd(self):
+        self.psd_from_file_1.asd_array[0] = 0.53544321
+        self.psd_from_array_1.asd_array[0] = 0.53544321
+        self.assertNotEqual(self.psd_from_file_1, self.psd_from_file_2)
+        self.assertNotEqual(self.psd_from_array_1, self.psd_from_array_2)
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/test/result_test.py b/test/result_test.py
index 3cc23eb2adc2753b4743fc95518e56ea6a7bd54d..2308bceac466372aa4e58d5fb23d361f1a99093b 100644
--- a/test/result_test.py
+++ b/test/result_test.py
@@ -240,6 +240,30 @@ class TestResult(unittest.TestCase):
         with self.assertRaises(TypeError):
             self.result.get_all_injection_credible_levels()
 
+    def test_kde(self):
+        kde = self.result.kde
+        import scipy.stats
+        self.assertEqual(type(kde), scipy.stats.kde.gaussian_kde)
+        self.assertEqual(kde.d, 2)
+
+    def test_posterior_probability(self):
+        sample = dict(x=0, y=0.1)
+        self.assertTrue(
+            isinstance(self.result.posterior_probability(sample), np.ndarray))
+        self.assertTrue(
+            len(self.result.posterior_probability(sample)), 1)
+        self.assertEqual(
+            self.result.posterior_probability(sample)[0],
+            self.result.kde([0, 0.1]))
+
+    def test_multiple_posterior_probability(self):
+        sample = [dict(x=0, y=0.1), dict(x=0.8, y=0)]
+        self.assertTrue(
+            isinstance(self.result.posterior_probability(sample), np.ndarray))
+        self.assertTrue(
+            all(self.result.posterior_probability(sample)
+                == self.result.kde([[0, 0.1], [0.8, 0]])))
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/test/series_test.py b/test/series_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..685164abf6f11192f77a089582bfb07571fc8d55
--- /dev/null
+++ b/test/series_test.py
@@ -0,0 +1,102 @@
+from __future__ import absolute_import
+import unittest
+import bilby
+import numpy as np
+
+
+class TestCoupledTimeAndFrequencySeries(unittest.TestCase):
+
+    def setUp(self):
+        self.duration = 2
+        self.sampling_frequency = 4096
+        self.start_time = -1
+        self.series = bilby.gw.series.CoupledTimeAndFrequencySeries(duration=self.duration,
+                                                                    sampling_frequency=self.sampling_frequency,
+                                                                    start_time=self.start_time)
+
+    def tearDown(self):
+        del self.duration
+        del self.sampling_frequency
+        del self.start_time
+        del self.series
+
+    def test_repr(self):
+        expected = 'CoupledTimeAndFrequencySeries(duration={}, sampling_frequency={}, start_time={})'\
+            .format(self.series.duration,
+                    self.series.sampling_frequency,
+                    self.series.start_time)
+        self.assertEqual(expected, repr(self.series))
+
+    def test_duration_from_init(self):
+        self.assertEqual(self.duration, self.series.duration)
+
+    def test_sampling_from_init(self):
+        self.assertEqual(self.sampling_frequency, self.series.sampling_frequency)
+
+    def test_start_time_from_init(self):
+        self.assertEqual(self.start_time, self.series.start_time)
+
+    def test_frequency_array_type(self):
+        self.assertIsInstance(self.series.frequency_array, np.ndarray)
+
+    def test_time_array_type(self):
+        self.assertIsInstance(self.series.time_array, np.ndarray)
+
+    def test_frequency_array_from_init(self):
+        expected = bilby.core.utils.create_frequency_series(sampling_frequency=self.sampling_frequency,
+                                                            duration=self.duration)
+        self.assertTrue(np.array_equal(expected, self.series.frequency_array))
+
+    def test_time_array_from_init(self):
+        expected = bilby.core.utils.create_time_series(sampling_frequency=self.sampling_frequency,
+                                                       duration=self.duration,
+                                                       starting_time=self.start_time)
+        self.assertTrue(np.array_equal(expected, self.series.time_array))
+
+    def test_frequency_array_setter(self):
+        new_sampling_frequency = 100
+        new_duration = 3
+        new_frequency_array = bilby.core.utils.create_frequency_series(sampling_frequency=new_sampling_frequency,
+                                                                       duration=new_duration)
+        self.series.frequency_array = new_frequency_array
+        self.assertTrue(np.array_equal(new_frequency_array, self.series.frequency_array))
+        self.assertLessEqual(np.abs(new_sampling_frequency - self.series.sampling_frequency), 1)
+        self.assertAlmostEqual(new_duration, self.series.duration)
+        self.assertAlmostEqual(self.start_time, self.series.start_time)
+
+    def test_time_array_setter(self):
+        new_sampling_frequency = 100
+        new_duration = 3
+        new_start_time = 4
+        new_time_array = bilby.core.utils.create_time_series(sampling_frequency=new_sampling_frequency,
+                                                             duration=new_duration,
+                                                             starting_time=new_start_time)
+        self.series.time_array = new_time_array
+        self.assertTrue(np.array_equal(new_time_array, self.series.time_array))
+        self.assertAlmostEqual(new_sampling_frequency, self.series.sampling_frequency, places=1)
+        self.assertAlmostEqual(new_duration, self.series.duration, places=1)
+        self.assertAlmostEqual(new_start_time, self.series.start_time, places=1)
+
+    def test_time_array_without_sampling_frequency(self):
+        self.series.sampling_frequency = None
+        self.series.duration = 4
+        with self.assertRaises(ValueError):
+            test = self.series.time_array
+
+    def test_time_array_without_duration(self):
+        self.series.sampling_frequency = 4096
+        self.series.duration = None
+        with self.assertRaises(ValueError):
+            test = self.series.time_array
+
+    def test_frequency_array_without_sampling_frequency(self):
+        self.series.sampling_frequency = None
+        self.series.duration = 4
+        with self.assertRaises(ValueError):
+            test = self.series.frequency_array
+
+    def test_frequency_array_without_duration(self):
+        self.series.sampling_frequency = 4096
+        self.series.duration = None
+        with self.assertRaises(ValueError):
+            test = self.series.frequency_array