Repository: StingraySoftware/stingray
Pull request #162 (2017-02-13 22:45:05)

[WIP] Map fitting astropy models
This pull request implements maximum likelihood and maximum a posteriori fitting of power spectra (and light curves, too).

It uses `astropy.modeling` for parametric models, `scipy.optimize` for optimization and `emcee` for MCMC.

It's currently work in progress, until I have time to finish the unit tests.

<!-- Reviewable:start -->
---
This change is [<img src="https://reviewable.io/review_button.svg" height="34" align="absmiddle" alt="Reviewable"/>](https://reviewable.io/reviews/stingraysoftware/stingray/162)
<!-- Reviewable:end -->



Patch:
diff --git a/.travis.yml b/.travis.yml
index 7804c50d..eede1553 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -6,6 +6,7 @@ python:
     - 3.3
     - 3.4
     - 3.5
+    - 3.6
 
 # Setting sudo to false opts in to Travis-CI container-based builds.
 sudo: false
@@ -27,12 +28,9 @@ env:
         # The following versions are the 'default' for tests, unless
         # overidden underneath. They are defined here in order to save having
         # to repeat them for all configurations.
-        - NUMPY_VERSION=1.9
-        - SCIPY_VERSION=0.17.0
         - ASTROPY_VERSION=stable
-        - H5PY_VERSION=2.5.0
         - SETUP_CMD='test'
-        - PIP_DEPENDENCIES=''
+        - PIP_DEPENDENCIES='emcee statsmodels corner'
         # For this package-template, we include examples of Cython modules,
         # so Cython is required for testing. If your package does not include
         # Cython code, you can set CONDA_DEPENDENCIES=''
@@ -50,13 +48,21 @@ matrix:
         - python: 2.7
           env: CONDA_DEPENDENCIES='Cython numpy scipy nose matplotlib'
 
-        - python: 3.4
+        - python: 3.5
           env: CONDA_DEPENDENCIES='Cython numpy scipy nose matplotlib'
 
+        # Try without using corner
+        - python: 3.5
+          env: PIP_DEPENDENCIES='emcee statsmodels' SETUP_CMD='test --coverage'
+
         # Do a coverage test in Python 2.
         - python: 2.7
           env: SETUP_CMD='test --coverage'
 
+        # Do a coverage test in Python 3.
+        - python: 3.5
+          env: SETUP_CMD='test --coverage'
+
         # Check for sphinx doc build warnings - we do this first because it
         # may run for a long time
         - python: 2.7
@@ -66,7 +72,7 @@ matrix:
         - python: 2.7
           env: ASTROPY_VERSION=development
 
-        - python: 3.4
+        - python: 3.5
           env: ASTROPY_VERSION=development
 
     allow_failures:
diff --git a/appveyor.yml b/appveyor.yml
index 1161fdb5..352b9aef 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -5,27 +5,32 @@ environment:
 
     global:
         PYTHON: "C:\\conda"
-        MINICONDA_VERSION: "latest"
+        MINICONDA_VERSION: ""
         CMD_IN_ENV: "cmd /E:ON /V:ON /C .\\ci-helpers\\appveyor\\windows_sdk.cmd"
         PYTHON_ARCH: "64" # needs to be set for CMD_IN_ENV to succeed. If a mix
                           # of 32 bit and 64 bit builds are needed, move this
                           # to the matrix section.
         CONDA_DEPENDENCIES: "scipy numpy nose h5py astropy matplotlib"
-        PIP_DEPENDENCIES: ""
+        PIP_DEPENDENCIES: "emcee statsmodels"
 
     matrix:
 
-        - PYTHON_VERSION: "2.7"
-          NUMPY_VERSION: "stable"
         - PYTHON_VERSION: "3.4"
-          NUMPY_VERSION: "stable"
+          NUMPY_VERSION: "1.11"
+          CONDA_PY: "34"
         - PYTHON_VERSION: "3.5"
+          NUMPY_VERSION: "1.11"
+          CONDA_PY: "35"
+        - PYTHON_VERSION: "2.7"
           NUMPY_VERSION: "stable"
-
+          CONDA_PY: "27"
+        - PYTHON_VERSION: "3.6"
+          NUMPY_VERSION: "stable"
+          CONDA_PY: "36"
 matrix:
     fast_finish: true
     allow_failures:
-        - PYTHON_VERSION: "3.5"
+        - PYTHON_VERSION: "3.6"
           NUMPY_VERSION: "stable"
 
 platform:
diff --git a/requirements.txt b/requirements.txt
index da5f9157..474e4d16 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
+six
 numpy
 scipy
 matplotlib
diff --git a/setup.cfg b/setup.cfg
index bc738a5e..8b52cc62 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -7,7 +7,7 @@ all_files = 1
 upload-dir = docs/_build/html
 show-response = 1
 
-[pytest]
+[tool:pytest]
 minversion = 2.2
 norecursedirs = build docs/_build
 doctest_plus = enabled
diff --git a/stingray/__init__.py b/stingray/__init__.py
index bf181896..63d16947 100644
--- a/stingray/__init__.py
+++ b/stingray/__init__.py
@@ -16,5 +16,6 @@ if not _ASTROPY_SETUP_:
     from stingray.utils import *
     from stingray.powerspectrum import *
     from stingray.crossspectrum import *
+    from stingray import *
     from stingray.exceptions import *
     from stingray.covariancespectrum import *
diff --git a/stingray/gti.py b/stingray/gti.py
index 2d4a8b38..f2e67ca9 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -467,11 +467,11 @@ def join_gtis(gti0, gti1):
 
     g0 = gti0.flatten()
     # Opening GTI: type = 1; Closing: type = -1
-    g0_type = np.asarray(list(zip(-np.ones(len(g0) / 2),
-                                  np.ones(len(g0) / 2))))
+    g0_type = np.asarray(list(zip(-np.ones(int(len(g0) / 2)),
+                                  np.ones(int(len(g0) / 2)))))
     g1 = gti1.flatten()
-    g1_type = np.asarray(list(zip(-np.ones(len(g1) / 2),
-                                  np.ones(len(g1) / 2))))
+    g1_type = np.asarray(list(zip(-np.ones(int(len(g1) / 2)),
+                                  np.ones(int(len(g1) / 2)))))
 
     g_all = np.append(g0, g1)
     g_type_all = np.append(g0_type, g1_type)
diff --git a/stingray/io.py b/stingray/io.py
index 5c97b0a9..54aea5b4 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -436,7 +436,7 @@ def _retrieve_hdf5_object(filename):
         attr_keys = hf.attrs.keys()
         data = {}
 
-        dset_copy = dset_keys[:]
+        dset_copy = list(dset_keys)[:]
         for key in dset_keys:
 
             # Make sure key hasn't been removed
@@ -453,7 +453,7 @@ def _retrieve_hdf5_object(filename):
                 else:
                     data[key] = hf[key].value
         
-        attr_copy = attr_keys[:]
+        attr_copy = list(attr_keys)[:]
         for key in attr_keys:
             
             # Make sure key hasn't been removed
diff --git a/stingray/modeling/__init__.py b/stingray/modeling/__init__.py
new file mode 100644
index 00000000..ff7f34cc
--- /dev/null
+++ b/stingray/modeling/__init__.py
@@ -0,0 +1,17 @@
+# Licensed under MIT license - see LICENSE.rst
+
+"""
+Library of Time Series Methods For Astronomical X-ray Data.
+"""
+
+# Affiliated packages may add whatever they like to this file, but
+# should keep this content at the top.
+# ----------------------------------------------------------------------------
+from .._astropy_init import *
+# ----------------------------------------------------------------------------
+
+# For egg_info test builds to pass, put package imports here.
+if not _ASTROPY_SETUP_:
+     from stingray.modeling.posterior import *
+     from stingray.modeling.parameterestimation import *
+     from stingray.modeling.scripts import *
\ No newline at end of file
diff --git a/stingray/modeling/parameterestimation.py b/stingray/modeling/parameterestimation.py
new file mode 100644
index 00000000..f8e6271d
--- /dev/null
+++ b/stingray/modeling/parameterestimation.py
@@ -0,0 +1,796 @@
+
+__all__ = ["OptimizationResults", "ParameterEstimation", "PSDParEst",
+           "SamplingResults"]
+
+
+# check whether matplotlib is installed for easy plotting
+try:
+    import matplotlib.pyplot as plt
+    from matplotlib.ticker import MaxNLocator
+    can_plot = True
+except ImportError:
+    can_plot = False
+
+
+# check whether emcee is installed for sampling
+try:
+    import emcee
+    can_sample = True
+except ImportError:
+    can_sample = False
+
+try:
+    import corner
+    use_corner = True
+except ImportError:
+    use_corner = False
+
+import numpy as np
+import scipy
+import scipy.optimize
+import scipy.stats
+import scipy.signal
+import copy
+
+try:
+    from statsmodels.tools.numdiff import approx_hess
+    comp_hessian = True
+except ImportError:
+    comp_hessian = False
+
+from stingray.modeling.posterior import Posterior, PSDPosterior, LogLikelihood
+from astropy.modeling.fitting import _fitter_to_model_params
+
+
+class OptimizationResults(object):
+
+    def __init__(self, lpost, res, neg=True):
+        """
+        Helper class that will contain the results of the regression.
+        Less fiddly than a dictionary.
+
+        Parameters
+        ----------
+        lpost: instance of Posterior or one of its subclasses
+            The object containing the function that is being optimized
+            in the regression
+
+
+        res: instance of scipy's OptimizeResult class
+            The object containing the results from a optimization run
+
+        """
+        self.neg = neg
+        self.result = res.fun
+        self.p_opt = res.x
+        self.model = lpost.model
+
+        self._compute_covariance(lpost, res)
+        self._compute_model(lpost)
+        self._compute_criteria(lpost)
+        self._compute_statistics(lpost)
+
+    def _compute_covariance(self, lpost, res):
+
+        if hasattr(res, "hess_inv"):
+            if not isinstance(res.hess_inv, np.ndarray):
+                self.cov = np.asarray(res.hess_inv.todense())
+            else:
+                self.cov = res.hess_inv
+
+            self.err = np.sqrt(np.diag(self.cov))
+        else:
+            if comp_hessian:
+                # calculate Hessian approximating with finite differences
+                print("Approximating Hessian with finite differences ...")
+
+                phess = approx_hess(self.p_opt, lpost)
+
+                self.cov = np.linalg.inv(phess)
+                self.err = np.sqrt(np.diag(np.abs(self.cov)))
+
+            else:
+                self.cov = None
+                self.err = None
+
+    def _compute_model(self, lpost):
+
+        _fitter_to_model_params(lpost.model, self.p_opt)
+
+        self.mfit = lpost.model(lpost.x)
+
+    def _compute_criteria(self, lpost):
+
+        if isinstance(lpost, Posterior):
+            self.deviance = -2.0*lpost.loglikelihood(self.p_opt, neg=False)
+        elif isinstance(lpost, LogLikelihood):
+            self.deviance = 2.0*self.result
+
+        # Akaike Information Criterion
+        self.aic = self.result+2.0*self.p_opt.shape[0]
+
+        # Bayesian Information Criterion
+        self.bic = self.result + self.p_opt.shape[0]*np.log(lpost.x.shape[0])
+
+        # Deviance Information Criterion
+        # TODO: Add Deviance Information Criterion
+
+    def _compute_statistics(self, lpost):
+        try:
+            self.mfit
+        except AttributeError:
+            self._compute_model(lpost)
+
+        self.merit = np.sum(((lpost.y-self.mfit)/self.mfit)**2.0)
+        self.dof = lpost.y.shape[0] - float(self.p_opt.shape[0])
+        self.sexp = 2.0*len(lpost.x)*len(self.p_opt)
+        self.ssd = np.sqrt(2.0*self.sexp)
+        self.sobs = np.sum(lpost.y-self.mfit)
+
+    def print_summary(self, lpost):
+
+        print("The best-fit model parameters plus errors are:")
+
+        fixed = [lpost.model.fixed[n] for n in lpost.model.param_names]
+        parnames = [n for n, f in zip(lpost.model.param_names, fixed) \
+                    if f is False]
+
+        for i, (x, y, p) in enumerate(zip(self.p_opt, self.err, parnames)):
+            print("%i) Parameter %s: %f.5 +/- %f.5f"%(i, p, x, y))
+
+        print("\n")
+
+        print("Fitting statistics: ")
+        print(" -- number of data points: %i"%(len(lpost.x)))
+
+        try:
+            self.deviance
+        except AttributeError:
+            self._compute_criteria(lpost)
+
+        print(" -- Deviance [-2 log L] D = %f.3"%self.deviance)
+        print(" -- The Akaike Information Criterion of the model is: " +
+              str(self.aic) + ".")
+
+        print(" -- The Bayesian Information Criterion of the model is: " +
+              str(self.bic) + ".")
+
+        try:
+            self.merit
+        except AttributeError:
+            self._compute_statistics(lpost)
+
+        print(" -- The figure-of-merit function for this model " +
+              " is: %f.5f"%self.merit +
+              " and the fit for %i dof is %f.3f"%(self.dof,
+                                                  self.merit/self.dof))
+
+        print(" -- Summed Residuals S = %f.5f"%self.sobs)
+        print(" -- Expected S ~ %f.5 +/- %f.5"%(self.sexp, self.ssd))
+        print(" -- merit function (SSE) M = %f.5f \n\n"%self.merit)
+
+        return
+
+
+class ParameterEstimation(object):
+
+    def __init__(self, fitmethod='BFGS', max_post=True):
+        """
+        Parameter estimation of two-dimensional data, either via
+        optimization or MCMC.
+        Note: optimization with bounds is not supported. If something like
+        this is required, define (uniform) priors in the ParametricModel
+        instances to be used below.
+
+        Parameters:
+        -----------
+        fitmethod: string, optional, default "L-BFGS-B"
+            Any of the strings allowed in scipy.optimize.minimize in
+            the method keyword. Sets the fit method to be used.
+
+        max_post: bool, optional, default True
+            If True, then compute the Maximum-A-Posteriori estimate. If False,
+            compute a Maximum Likelihood estimate.
+        """
+
+        self.fitmethod = fitmethod
+
+        self.max_post = max_post
+
+    def fit(self, lpost, t0, neg=True, scipy_optimize_options=None):
+        """
+        Do either a Maximum A Posteriori or Maximum Likelihood
+        fit to the data.
+
+        Parameters:
+        -----------
+        lpost: Posterior (or subclass) instance
+            and instance of class Posterior or one of its subclasses
+            that defines the function to be minized (either in loglikelihood
+            or logposterior)
+
+        t0 : {list | numpy.ndarray}
+            List/array with set of initial parameters
+
+        neg : bool, optional, default True
+            Boolean to be passed to `lpost`, setting whether to use the
+            *negative* posterior or the *negative* log-likelihood. Since
+            `Posterior` and `LogLikelihood` objects are generally defined in
+
+        scipy_optimize_options : dict, optional, default None
+            A dictionary with options for `scipy.optimize.minimize`,
+            directly passed on as keyword arguments.
+
+        Returns:
+        --------
+        fitparams: dict
+            A dictionary with the fit results
+            TODO: Add description of keywords in the class!
+        """
+
+        if not isinstance(lpost, Posterior) and not isinstance(lpost,
+                                                               LogLikelihood):
+            raise TypeError("lpost must be a subclass of "
+                            "Posterior or LogLikelihoood.")
+
+        if not len(t0) == lpost.npar:
+            raise ValueError("Parameter set t0 must be of right "
+                             "length for model in lpost.")
+
+        if scipy.__version__ < "0.10.0":
+            args = [neg]
+        else:
+            args = (neg,)
+
+        if not scipy_optimize_options:
+            scipy_optimize_options = {}
+
+        # different commands for different fitting methods,
+        # at least until scipy 0.11 is out
+        funcval = 100.0
+        i = 0
+        while funcval == 100 or funcval == 200 or \
+                funcval == 0.0 or not np.isfinite(funcval):
+
+            if i > 20:
+                raise Exception("Fitting unsuccessful!")
+            # perturb parameters slightly
+            t0_p = np.random.multivariate_normal(t0, np.diag(np.abs(t0)/100.))
+
+            # if max_post is True, do the Maximum-A-Posteriori Fit
+            if self.max_post:
+                opt = scipy.optimize.minimize(lpost, t0_p,
+                                              method=self.fitmethod,
+                                              args=args, tol=1.e-10,
+                                              **scipy_optimize_options)
+
+            # if max_post is False, then do a Maximum Likelihood Fit
+            else:
+                if isinstance(lpost, Posterior):
+                    # This could be a `Posterior` object
+                    opt = scipy.optimize.minimize(lpost.loglikelihood, t0_p,
+                                                  method=self.fitmethod,
+                                                  args=args, tol=1.e-10,
+                                                 **scipy_optimize_options)
+
+                elif isinstance(lpost, LogLikelihood):
+                    # Except this could be a `LogLikelihood object
+                    # In which case, use the evaluate function
+                    # if it's not either, give up and break!
+                    opt = scipy.optimize.minimize(lpost.evaluate, t0_p,
+                                                  method=self.fitmethod,
+                                                  args=args, tol=1.e-10,
+                                                  **scipy_optimize_options)
+                else:
+                    raise TypeError("lpost must be a Posterior or LogLikelihood"
+                                    " object.")
+
+            funcval = opt.fun
+            i += 1
+
+        res = OptimizationResults(lpost, opt, neg=neg)
+
+        return res
+
+    def compute_lrt(self, lpost1, t1, lpost2, t2, neg=True, max_post=False):
+        """
+        This function computes the Likelihood Ratio Test between two
+        nested models.
+
+        Parameters
+        ----------
+        lpost1 : object of a subclass of Posterior
+            The posterior object for model 1
+
+        t1 : iterable
+            The starting parameters for model 1
+
+        lpost2 : object of a subclass of Posterior
+            The posterior object for model 2
+
+        t2 : iterable
+            The starting parameters for model 2
+
+        neg : bool, optional, default True
+            Boolean flag to decide whether to use the negative log-likelihood
+            or log-posterior
+
+        max_post: bool, optional, default False
+            If True, set the internal state to do the optimization with the
+            log-likelihood rather than the log-posterior.
+
+        """
+
+        self.max_post = max_post
+
+        # fit data with both models
+        res1 = self.fit(lpost1, t1, neg=neg)
+        res2 = self.fit(lpost2, t2, neg=neg)
+
+        # compute log likelihood ratio as difference between the deviances
+        lrt = res1.deviance - res2.deviance
+
+        return lrt
+
+    def sample(self, lpost, t0, cov=None,
+               nwalkers=500, niter=100, burnin=100, threads=1,
+               print_results=True, plot=False, namestr="test"):
+        """
+        Sample the posterior distribution defined in `lpost` using MCMC.
+        Here we use the `emcee` package, but other implementations could
+        in principle be used.
+
+        Parameters
+        ----------
+        t0 : iterable
+            list or array containing the starting parameters. Its length
+            must match `lpost.model.npar`.
+
+
+        nwalkers : int
+            The number of walkers (chains) to use during the MCMC procedure.
+            The more walkers are used, the slower the estimation will be, but
+            the better the final distribution is likely to be.
+
+        niter : int
+            The number of iterations to run the MCMC chains for. The larger this
+            number, the longer the estimation will take, but the higher the
+            chance that the walkers have actually converged on the true
+            posterior distribution.
+
+        burnin : int
+            The number of iterations to run the walkers before convergence is
+            assumed to have occurred. This part of the chain will be discarded
+            before sampling from what is then assumed to be the posterior
+            distribution desired.
+
+        threads : int
+            The number of threads for parallelization.
+            Default is 1, i.e. no parallelization
+
+        print_results : bool
+            Boolean flag setting whether the results of the MCMC run should
+            be printed to standard output. Default: True
+
+        plot : bool
+            Boolean flag setting whether summary plots of the MCMC chains
+            should be produced. Default: False
+
+        namestr : str
+            Optional string for output file names for the plotting.
+
+        Returns
+        -------
+
+        res : SamplingResults object
+
+        """
+        if not can_sample:
+            raise ImportError("emcee not installed! Can't sample!")
+
+        ndim = len(t0)
+
+        if cov is None:
+            # do a MAP fitting step to find good starting positions for
+            # the sampler
+            res = self.fit(lpost, t0, neg=True)
+            cov = res.cov
+        # sample random starting positions for each walker from
+        # a multivariate Gaussian
+        p0 = np.array([np.random.multivariate_normal(t0, cov) for
+                       i in range(nwalkers)])
+
+        # initialize the sampler
+        sampler = emcee.EnsembleSampler(nwalkers, ndim, lpost, args=[False],
+                                        threads=threads)
+
+        # run the burn-in
+        pos, prob, state = sampler.run_mcmc(p0, burnin)
+
+        sampler.reset()
+
+        # do the actual MCMC run
+        _, _, _ = sampler.run_mcmc(pos, niter, rstate0=state)
+
+        res = SamplingResults(sampler)
+
+        if print_results:
+            res.print_results()
+
+        if plot:
+            res.plot_results(namestr + "_corner.pdf")
+
+        return res
+
+
+class SamplingResults(object):
+
+    def __init__(self, sampler, ci_min=0.05, ci_max=0.95):
+        """
+        Helper class that will contain the results of the sampling
+        in a handly format.
+        Less fiddly than a dictionary.
+
+        Parameters
+        ----------
+        sampler: emcee.EnsembleSampler object
+            The object containing the sampler that's done all the work.
+
+        ci_min: float out of [0,1]
+            The lower bound percentile for printing confidence intervals
+            on the parameters
+
+        ci_max: float out of [0,1]
+            The upper bound percentile for printing confidence intervals
+            on the parameters
+
+        """
+
+        # store all the samples
+        self.samples = sampler.flatchain
+
+        self.nwalkers = np.float(sampler.chain.shape[0])
+        self.niter = np.float(sampler.iterations)
+
+        # store number of dimensions
+        self.ndim = sampler.dim
+
+        # compute and store acceptance fraction
+        self.acceptance = np.nanmean(sampler.acceptance_fraction)
+        self.L = self.acceptance*self.samples.shape[0]
+
+        self._check_convergence(sampler)
+        self._infer(ci_min, ci_max)
+
+    def _check_convergence(self, sampler):
+
+        # compute and store autocorrelation time
+        try:
+            self.acor = sampler.acor
+        except emcee.autocorr.AutocorrError:
+            print("Chains too short to compute autocorrelation lengths.")
+
+        self.rhat = self._compute_rhat(sampler)
+
+    def _compute_rhat(self, sampler):
+
+        # between-sequence variance
+        mean_samples_iter = np.nanmean(sampler.chain, axis=1)
+
+        # mean over the means over iterations: (self.ndim)
+        mean_samples = np.nanmean(sampler.chain, axis=(0,1))
+
+        # now compute between-sequence variance
+        bb = (self.niter / (self.nwalkers - 1)) * np.sum((mean_samples_iter -
+                                                          mean_samples)**2.,
+                                                         axis=0)
+
+        # compute variance of each chain
+        var_samples = np.nanvar(sampler.chain, axis=1)
+
+        # compute mean of variance
+        ww = np.nanmean(var_samples, axis=0)
+
+        # compute weighted average of ww and bb:
+        rhat = ((self.niter - 1) / self.niter) * ww + (1 / self.niter) * bb
+
+        return rhat
+
+    def _infer(self, ci_min=0.05, ci_max=0.95):
+        self.mean = np.mean(self.samples, axis=0)
+        self.std = np.std(self.samples, axis=0)
+        self.ci = np.percentile(self.samples, [ci_min, ci_max], axis=0)
+
+    def print_results(self):
+        """
+        Print results of the MCMC run.
+
+        """
+
+        print("-- The acceptance fraction is: %f.5"%self.acceptance)
+        try:
+            print("-- The autocorrelation time is: %f.5"%self.acor)
+        except AttributeError:
+            pass
+        print("R_hat for the parameters is: " + str(self.rhat))
+
+        print("-- Posterior Summary of Parameters: \n")
+        print("parameter \t mean \t\t sd \t\t 5% \t\t 95% \n")
+        print("---------------------------------------------\n")
+        for i in range(self.ndim):
+            print("theta[" + str(i) + "] \t " +
+                  str(self.mean[i]) + "\t" + str(self.std[i]) + "\t" +
+                  str(self.ci[0, i]) + "\t" + str(self.ci[1, i]) + "\n")
+
+    def plot_results(self, filename, nsamples=1000):
+        """
+        Plot some results in a triangle plot.
+        If installed, will use `corner` for the plotting
+        (available here https://github.com/dfm/corner.py or
+        through pip), if not, uses its own code to make a triangle
+        plot.
+
+        Parameters
+        ----------
+
+        filename: str
+            Name of the output file with the figure
+        nsamples: int
+            The maximum number of samples used for plotting.
+
+        """
+        assert can_plot, "Need to have matplotlib installed for plotting"
+        if use_corner:
+            corner.corner(self.samples, labels=None,
+                          quantiles=[0.16, 0.5, 0.84],
+                          show_titles=True, title_args={"fontsize": 12})
+
+        else:
+            fig = plt.figure(figsize=(15, 15))
+            plt.subplots_adjust(top=0.925, bottom=0.025,
+                                left=0.025, right=0.975,
+                                wspace=0.2, hspace=0.2)
+
+            ind_all = np.random.choice(np.arange(self.samples.shape[0]),
+                                       size=nsamples)
+            samples = self.samples[ind_all]
+            for i in range(self.ndim):
+                for j in range(self.ndim):
+                    xmin, xmax = samples[:, j].min(), samples[:, j].max()
+                    ymin, ymax = samples[:, i].min(), samples[:, i].max()
+                    ax = fig.add_subplot(self.ndim, self.ndim, i*self.ndim+j+1)
+
+                    ax.xaxis.set_major_locator(MaxNLocator(5))
+                    ax.ticklabel_format(style="sci", scilimits=(-2, 2))
+
+                    if i == j:
+                        ntemp, binstemp, patchestemp = \
+                            ax.hist(samples[:, i], 30, normed=True,
+                                    histtype='stepfilled')
+                        ax.axis([ymin, ymax, 0, max(ntemp)*1.2])
+
+                    else:
+
+                        ax.axis([xmin, xmax, ymin, ymax])
+
+                        # make a scatter plot first
+                        ax.scatter(samples[:, j], samples[:, i], s=7)
+                        # then add contours
+                        xmin, xmax = samples[:, j].min(), samples[:, j].max()
+                        ymin, ymax = samples[:, i].min(), samples[:, i].max()
+
+                        # Perform Kernel density estimate on data
+                        try:
+                            xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
+                            positions = np.vstack([xx.ravel(), yy.ravel()])
+                            values = np.vstack([samples[:, j], samples[:, i]])
+                            kernel = scipy.stats.gaussian_kde(values)
+                            zz = np.reshape(kernel(positions).T, xx.shape)
+
+                            ax.contour(xx, yy, zz, 7)
+                        except ValueError:
+                            print("Not making contours.")
+
+        plt.savefig(filename, format='pdf')
+        plt.close()
+        return
+
+
+class PSDParEst(ParameterEstimation):
+
+    def __init__(self, ps, fitmethod='BFGS', max_post=True):
+
+        self.ps = ps
+        ParameterEstimation.__init__(self, fitmethod=fitmethod,
+                                     max_post=max_post)
+
+    def fit(self, lpost, t0, neg=True):
+
+        self.lpost = lpost
+
+        res = ParameterEstimation.fit(self, self.lpost, t0, neg=neg)
+
+        res.maxpow, res.maxfreq, res.maxind = \
+            self._compute_highest_outlier(self.lpost, res)
+
+        return res
+
+    def sample(self, lpost, t0, cov=None,
+               nwalkers=500, niter=100, burnin=100, threads=1,
+               print_results=True, plot=False, namestr="test"):
+
+        self.lpost = lpost
+
+        fit_res = ParameterEstimation.fit(self, self.lpost, t0, neg=True)
+
+        res = ParameterEstimation.sample(self, self.lpost, fit_res.p_opt,
+                                         cov=fit_res.cov,
+                                         nwalkers=nwalkers,
+                                         niter=niter, burnin=burnin,
+                                         threads=threads,
+                                         print_results=print_results, plot=plot,
+                                         namestr=namestr)
+
+        return res
+
+    def _compute_highest_outlier(self, lpost, res, nmax=1):
+
+        residuals = 2.0 * lpost.y/ res.mfit
+
+        ratio_sort = copy.copy(residuals)
+        ratio_sort.sort()
+        max_y = ratio_sort[-nmax:]
+
+        max_x = np.zeros(max_y.shape[0])
+        max_ind = np.zeros(max_y.shape[0])
+
+        for i, my in enumerate(max_y):
+            max_x[i], max_ind[i] = self._find_outlier(lpost.x, residuals, my)
+
+        return max_y, max_x, max_ind
+
+    @staticmethod
+    def _find_outlier(xdata, ratio, max_y):
+        max_ind = np.where(ratio == max_y)[0][0]+1
+        max_x = xdata[max_ind]
+
+        return max_x, max_ind
+
+    def plotfits(self, res1, res2=None, save_plot=False,
+                 namestr='test', log=False):
+
+        if not can_plot:
+            print("No matplotlib imported. Can't plot!")
+        else:
+            # make a figure
+            f = plt.figure(figsize=(12, 10))
+            # adjust subplots such that the space between the top and bottom
+            # of each are zero
+            plt.subplots_adjust(hspace=0.0, wspace=0.4)
+
+            # first subplot of the grid, twice as high as the other two
+            # This is the periodogram with the two fitted models overplotted
+            s1 = plt.subplot2grid((4, 1), (0, 0), rowspan=2)
+
+            if log:
+                logx = np.log10(self.ps.freq)
+                logy = np.log10(self.ps.power)
+                logpar1 = np.log10(res1.mfit)
+
+                p1, = s1.plot(logx, logy, color='black', linestyle='steps-mid')
+                p2, = s1.plot(logx, logpar1, color='blue', lw=2)
+                s1.set_xlim([min(logx), max(logx)])
+                s1.set_ylim([min(logy)-1.0, max(logy)+1])
+                if self.ps.norm == "leahy":
+                    s1.set_ylabel('log(Leahy-Normalized Power)', fontsize=18)
+                elif self.ps.norm == "rms":
+                    s1.set_ylabel('log(RMS-Normalized Power)', fontsize=18)
+                else:
+                    s1.set_ylabel("log(Power)", fontsize=18)
+
+            else:
+                p1, = s1.plot(self.ps.freq, self.ps.power,
+                              color='black', linestyle='steps-mid')
+                p2, = s1.plot(self.ps.freq, res1.mfit,
+                              color='blue', lw=2)
+
+                s1.set_xscale("log")
+                s1.set_yscale("log")
+
+                s1.set_xlim([min(self.ps.freq), max(self.ps.freq)])
+                s1.set_ylim([min(self.ps.freq)/10.0,
+                             max(self.ps.power)*10.0])
+
+                if self.ps.norm == "leahy":
+                    s1.set_ylabel('Leahy-Normalized Power', fontsize=18)
+                elif self.ps.norm == "rms":
+                    s1.set_ylabel('RMS-Normalized Power', fontsize=18)
+                else:
+                    s1.set_ylabel("Power", fontsize=18)
+
+            if res2 is not None:
+                if log:
+                    logpar2 = np.log10(res2.mfit)
+                    p3, = s1.plot(logx, logpar2, color='red', lw=2)
+                else:
+                    p3, = s1.plot(self.ps.freq, res2.mfit,
+                                  color='red', lw=2)
+                s1.legend([p1, p2, p3], ["data", "model 1 fit", "model 2 fit"])
+            else:
+                s1.legend([p1, p2], ["data", "model fit"])
+
+            s1.set_title("Periodogram and fits for data set " + namestr,
+                         fontsize=18)
+
+            # second subplot: power/model for Power law and straight line
+            s2 = plt.subplot2grid((4, 1), (2, 0), rowspan=1)
+            pldif = self.ps.power / res1.mfit
+
+            s2.set_ylabel("Residuals, \n first model",
+                          fontsize=18)
+
+            if log:
+                s2.plot(logx, pldif, color='black', linestyle='steps-mid')
+                s2.plot(logx, np.ones(self.ps.freq.shape[0]),
+                        color='blue', lw=2)
+                s2.set_xlim([min(logx), max(logx)])
+                s2.set_ylim([min(pldif), max(pldif)])
+
+            else:
+                s2.plot(self.ps.freq, pldif, color='black',
+                        linestyle='steps-mid')
+                s2.plot(self.ps.power, np.ones(self.ps.power.shape[0]),
+                        color='blue', lw=2)
+
+                s2.set_xscale("log")
+                s2.set_yscale("log")
+                s2.set_xlim([min(self.ps.freq), max(self.ps.freq)])
+                s2.set_ylim([min(pldif), max(pldif)])
+
+            if res2 is not None:
+                bpldif = self.ps.power/res2.mfit
+
+            # third subplot: power/model for bent power law and straight line
+                s3 = plt.subplot2grid((4, 1), (3, 0), rowspan=1)
+
+                if log:
+                    s3.plot(logx, bpldif, color='black', linestyle='steps-mid')
+                    s3.plot(logx, np.ones(len(self.ps.freq)),
+                            color='red', lw=2)
+                    s3.axis([min(logx), max(logx), min(bpldif), max(bpldif)])
+                    s3.set_xlabel("log(Frequency) [Hz]", fontsize=18)
+
+                else:
+                    s3.plot(self.ps.freq, bpldif,
+                            color='black', linestyle='steps-mid')
+                    s3.plot(self.ps.freq, np.ones(len(self.ps.freq)),
+                            color='red', lw=2)
+                    s3.set_xscale("log")
+                    s3.set_yscale("log")
+                    s3.set_xlim([min(self.ps.freq), max(self.ps.freq)])
+                    s3.set_ylim([min(bpldif), max(bpldif)])
+                    s3.set_xlabel("Frequency [Hz]", fontsize=18)
+
+                s3.set_ylabel("Residuals, \n second model",
+                              fontsize=18)
+
+            else:
+                if log:
+                    s2.set_xlabel("log(Frequency) [Hz]", fontsize=18)
+                else:
+                    s2.set_xlabel("Frequency [Hz]", fontsize=18)
+
+            ax = plt.gca()
+
+            for label in ax.get_xticklabels() + ax.get_yticklabels():
+                label.set_fontsize(14)
+
+            # make sure xticks are taken from first plots, but don't
+            # appear there
+            plt.setp(s1.get_xticklabels(), visible=False)
+
+            if save_plot:
+                # save figure in png file and close plot device
+                plt.savefig(namestr + '_ps_fit.png', format='png')
+
+        return
diff --git a/stingray/modeling/posterior.py b/stingray/modeling/posterior.py
new file mode 100644
index 00000000..4986e4d3
--- /dev/null
+++ b/stingray/modeling/posterior.py
@@ -0,0 +1,580 @@
+from __future__ import division
+
+import abc
+
+import numpy as np
+from scipy.special import gamma as scipy_gamma
+from scipy.special import gammaln as scipy_gammaln
+from astropy.modeling.fitting import _fitter_to_model_params
+from astropy.modeling import models
+
+from stingray import Lightcurve, Powerspectrum
+
+# TODO: Add checks and balances to code
+
+#from stingray.modeling.parametricmodels import logmin
+
+__all__ = ["set_logprior", "Posterior", "PSDPosterior", "LogLikelihood",
+           "PSDLogLikelihood", "GaussianLogLikelihood",
+           "PoissonPosterior", "GaussianPosterior",
+           "PriorUndefinedError", "LikelihoodUndefinedError"]
+
+logmin = -10000000000000000.0
+
+class PriorUndefinedError(Exception):
+    pass
+
+class LikelihoodUndefinedError(Exception):
+    pass
+
+class IncorrectParameterError(Exception):
+    pass
+
+def set_logprior(lpost, priors):
+    """
+    This function constructs the `logprior` method required to successfully
+    use a `Posterior` object.
+
+    All instances of lass `Posterior` and its subclasses require to implement a
+    `logprior` methods. However, priors are strongly problem-dependent and
+    therefore usually user-defined.
+
+    This function allows for setting the `logprior` method on any instance
+    of class `Posterior` efficiently by allowing the user to pass a
+    dictionary of priors and an instance of class `Posterior`.
+
+    Parameters
+    ----------
+    lpost : Posterior object
+        An instance of class Posterior or any of its subclasses
+
+    priors : dictionary
+        A dictionary containing the prior definitions. Keys are parameter
+        names as defined by the model used in `lpost`. Items are functions
+        that take a parameter as input and return the log-prior probability
+         of that parameter.
+
+    Returns
+    -------
+        logprior : function
+            The function definition for the prior
+
+    Example
+    -------
+    Make a light curve and power spectrum
+
+    >>> photon_arrivals = np.sort(np.random.uniform(0,1000, size=10000))
+    >>> lc = Lightcurve.make_lightcurve(photon_arrivals, dt=1.0)
+    >>> ps = Powerspectrum(lc, norm="frac")
+
+    Define the model
+
+    >>> pl = models.PowerLaw1D()
+    >>> pl.x_0.fixed = True
+
+    Instantiate the posterior:
+
+    >>> lpost = PSDPosterior(ps, pl)
+
+    Define the priors:
+
+    >>> p_alpha = lambda alpha: ((-1. <= alpha) & (alpha <= 5.))
+    >>> p_amplitude = lambda amplitude: ((-10 <= np.log(amplitude)) &
+    ...                                 ((np.log(amplitude) <= 10.0)))
+    >>> priors = {"alpha":p_alpha, "amplitude":p_amplitude}
+
+    Set the logprior method in the lpost object:
+
+    >>> lpost.logprior = set_logprior(lpost, priors)
+    """
+
+    # get the number of free parameters in the model
+    #free_params = [p for p in lpost.model.param_names if not
+    #                getattr(lpost.model, p).fixed]
+    free_params = [key for key, l in lpost.model.fixed.items() if not l]
+
+    # define the logprior
+    def logprior(t0, neg=False):
+        """
+        The logarithm of the prior distribution for the
+        model defined in self.model.
+
+        Parameters:
+        ------------
+        t0: {list | numpy.ndarray}
+            The list with parameters for the model
+
+        Returns:
+        --------
+        logp: float
+            The logarithm of the prior distribution for the model and
+            parameters given.
+        """
+
+        if len(t0) != len(free_params):
+            raise IncorrectParameterError("The number of parameters passed into "
+                                          "the prior does not match the number "
+                                          "of parameters in the model.")
+
+        logp = 0.0 # initialize log-prior
+        ii = 0 # counter for the variable parameter
+
+        # loop through all parameter names, but only compute
+        # prior for those that are not fixed
+        # Note: need to do it this way to preserve order of parameters
+        # correctly!
+        for pname in lpost.model.param_names:
+            if not lpost.model.fixed[pname]:
+                logp += np.log(priors[pname](t0[ii]))
+                ii += 1
+
+        if not np.isfinite(logp):
+            logp = logmin
+
+        if neg:
+            return -logp
+        else:
+            return logp
+
+    return logprior
+
+
+class LogLikelihood(object):
+    __metaclass__ = abc.ABCMeta
+
+    def __init__(self, x, y, model):
+        """
+        x : iterable
+            x-coordinate of the data. Could be multi-dimensional.
+
+        y : iterable
+            y-coordinate of the data. Could be multi-dimensional.
+
+        model: probably astropy.modeling.FittableModel instance
+            Your model
+        """
+        self.x = x
+        self.y = y
+
+        self.model = model
+
+    @abc.abstractmethod
+    def evaluate(self, parameters):
+        """
+        This is where you define your log-likelihood. Do this!
+
+        """
+        pass
+
+    def __call__(self, parameters, neg=False):
+        return self.evaluate(parameters, neg)
+
+
+class GaussianLogLikelihood(LogLikelihood):
+
+    def __init__(self, x, y, yerr, model):
+        """
+        A Gaussian likelihood.
+
+        Parameters
+        ----------
+        x : iterable
+            x-coordinate of the data
+
+        y : iterable
+            y-coordinte of the data
+
+        yerr: iterable
+            the uncertainty on the data, as standard deviation
+
+        model: an Astropy Model instance
+            The model to use in the likelihood.
+
+        """
+
+        self.x = x
+        self.y = y
+        self.yerr = yerr
+        self.model = model
+
+        self.npar = 0
+        for pname in self.model.param_names:
+            if not self.model.fixed[pname]:
+                self.npar += 1
+
+
+    def evaluate(self, pars, neg=False):
+        if np.size(pars) != self.npar:
+            raise IncorrectParameterError("Input parameters must" +
+                                          " match model parameters!")
+
+        _fitter_to_model_params(self.model, pars)
+
+        mean_model = self.model(self.x)
+
+        loglike = np.sum(-0.5*np.log(2.*np.pi) - np.log(self.yerr) -
+                         (self.y-mean_model)**2/(2.*self.yerr**2))
+
+        if not np.isfinite(loglike):
+            loglike = logmin
+
+        if neg:
+            return -loglike
+        else:
+            return loglike
+
+
+class PoissonLogLikelihood(LogLikelihood):
+
+    def __init__(self, x, y, model):
+        """
+        A Gaussian likelihood.
+
+        Parameters
+        ----------
+        x : iterable
+            x-coordinate of the data
+
+        y : iterable
+            y-coordinte of the data
+
+        model: an Astropy Model instance
+            The model to use in the likelihood.
+
+        """
+
+        self.x = x
+        self.y = y
+        self.model = model
+        self.npar = 0
+        for pname in self.model.param_names:
+            if not self.model.fixed[pname]:
+                self.npar += 1
+
+    def evaluate(self, pars, neg=False):
+
+        if np.size(pars) != self.npar:
+            raise IncorrectParameterError("Input parameters must" +
+                                          " match model parameters!")
+
+        _fitter_to_model_params(self.model, pars)
+
+        mean_model = self.model(self.x)
+
+        loglike = np.sum(-mean_model + self.y*np.log(mean_model) \
+               - scipy_gammaln(self.y + 1.))
+
+        if not np.isfinite(loglike):
+            loglike = logmin
+
+        if neg:
+            return -loglike
+        else:
+            return loglike
+
+
+class PSDLogLikelihood(LogLikelihood):
+
+    def __init__(self, x, y, model, m=1):
+        """
+        A Gaussian likelihood.
+
+        Parameters
+        ----------
+        x : iterable
+            x-coordinate of the data
+
+        y : iterable
+            y-coordinte of the data
+
+        model: an Astropy Model instance
+            The model to use in the likelihood.
+
+        m : int
+            1/2 of the degrees of freedom
+
+        """
+
+        self.x = x
+        self.y = y
+        self.model = model
+        self.m = m
+        self.npar = 0
+        for pname in self.model.param_names:
+            if not self.model.fixed[pname]:
+                self.npar += 1
+
+    def evaluate(self, pars, neg=False):
+
+        if np.size(pars) != self.npar:
+            raise IncorrectParameterError("Input parameters must" +
+                                          " match model parameters!")
+
+        _fitter_to_model_params(self.model, pars)
+
+        mean_model = self.model(self.x)
+
+        if self.m == 1:
+            loglike = -np.sum(np.log(mean_model)) - \
+                      np.sum(self.y/mean_model)
+
+        else:
+            loglike = -2.0*self.m*(np.sum(np.log(mean_model)) +
+                               np.sum(self.y/mean_model) +
+                               np.sum((2.0 / (2. * self.m) - 1.0) *
+                                      np.log(self.y)))
+
+
+        if not np.isfinite(loglike):
+            loglike = logmin
+
+        if neg:
+            return -loglike
+        else:
+            return loglike
+
+
+class Posterior(object):
+
+    def __init__(self, x, y, model):
+        """
+        Define a posterior object.
+
+        The posterior describes the Bayesian probability distribution of
+        a set of parameters $\theta$ given some observed data $D$ and
+        some prior assumptions $I$.
+
+        It is defined as
+
+            $p(\theta | D, I) = p(D | \theta, I) p(\theta | I)/p(D| I)
+
+        where $p(D | \theta, I)$ describes the likelihood, i.e. the
+        sampling distribution of the data and the (parametric) model, and
+        $p(\theta | I)$ describes the prior distribution, i.e. our information
+        about the parameters $\theta$ before we gathered the data.
+        The marginal likelihood $p(D| I)$ describes the probability of
+        observing the data given the model assumptions, integrated over the
+        space of all parameters.
+
+        Parameters
+        ----------
+        x : iterable
+            The abscissa or independent variable of the data. This could
+            in principle be a multi-dimensional array.
+
+        y : iterable
+            The ordinate or dependent variable of the data.
+
+        model: astropy.modeling.models class instance
+            The parametric model supposed to represent the data. For details
+            see the astropy.modeling documentation
+
+
+        References
+        ----------
+
+        * Sivia, D. S., and J. Skilling. "Data Analysis:
+            A Bayesian Tutorial. 2006."
+        * Gelman, Andrew, et al. Bayesian data analysis. Vol. 2. Boca Raton,
+            FL, USA: Chapman & Hall/CRC, 2014.
+        * von Toussaint, Udo. "Bayesian inference in physics."
+            Reviews of Modern Physics 83.3 (2011): 943.
+        *Hogg, David W. "Probability Calculus for inference".
+            arxiv: 1205.4446
+
+        """
+        self.x = x
+        self.y = y
+
+        self.model = model
+
+        self.npar = 0
+        for pname in self.model.param_names:
+            if not self.model.fixed[pname]:
+                self.npar += 1
+
+    def logposterior(self, t0, neg=False):
+
+        if not hasattr(self, "logprior"):
+            raise PriorUndefinedError("There is no prior implemented. " +
+                                      "Cannot calculate posterior!")
+
+        if not hasattr(self, "loglikelihood"):
+            raise LikelihoodUndefinedError("There is no likelihood implemented. " +
+                                           "Cannot calculate posterior!")
+
+        lpost = self.loglikelihood(t0) + self.logprior(t0)
+
+        if neg is True:
+            return -lpost
+        else:
+            return lpost
+
+    def __call__(self, t0, neg=False):
+        return self.logposterior(t0, neg=neg)
+
+
+class PSDPosterior(Posterior):
+
+    def __init__(self, ps, model, priors=None):
+        """
+        Posterior distribution for power spectra.
+        Uses an exponential distribution for the errors in the likelihood,
+        or a $\chi^2$ distribution with $2M$ degrees of freedom, where $M$ is
+        the number of frequency bins or power spectra averaged in each bin.
+
+
+        Parameters
+        ----------
+        ps: {Powerspectrum | AveragedPowerspectrum} instance
+            the Powerspectrum object containing the data
+
+        model: instance of any subclass of parameterclass.ParametricModel
+            The model for the power spectrum. Note that in order to define
+            the posterior properly, the ParametricModel subclass must be
+            instantiated with the hyperpars parameter set, or there won't
+            be a prior to be calculated! If all this object is used
+            for a maximum likelihood-style analysis, no prior is required.
+
+        priors : dict of form {"parameter name": function}, optional
+            A dictionary with the definitions for the prior probabilities.
+            For each parameter in `model`, there must be a prior defined with
+            a key of the exact same name as stored in `model.param_names`.
+            The item for each key is a function definition defining the prior
+            (e.g. a lambda function or a `scipy.stats.distribution.pdf`.
+            If `priors = None`, then no prior is set. This means priors need
+            to be added by hand using the `set_logprior` function defined in
+            this module. Note that it is impossible to call the posterior object
+            itself or the `self.logposterior` method without defining a prior.
+
+        Attributes
+        ----------
+        ps: {Powerspectrum | AveragedPowerspectrum} instance
+            the Powerspectrum object containing the data
+
+        m: int
+            The number of averaged periodograms or frequency bins in ps.
+            Useful for binned/averaged periodograms, since the value of
+            m will change the likelihood function!
+
+        x: numpy.ndarray
+            The independent variable (list of frequencies) stored in ps.freq
+
+        y: numpy.ndarray
+            The dependent variable (list of powers) stored in ps.power
+
+        model: instance of any subclass of parameterclass.ParametricModel
+               The model for the power spectrum. Note that in order to define
+               the posterior properly, the ParametricModel subclass must be
+               instantiated with the hyperpars parameter set, or there won't
+               be a prior to be calculated! If all this object is used
+               for a maximum likelihood-style analysis, no prior is required.
+
+        """
+        self.loglikelihood = PSDLogLikelihood(ps.freq,
+                                              ps.power,
+                                              model, m=ps.m)
+
+        self.m = ps.m
+        Posterior.__init__(self, ps.freq, ps.power, model)
+
+        if not priors is None:
+            self.logprior = set_logprior(self, priors)
+
+
+class PoissonPosterior(Posterior):
+
+    def __init__(self, x, y, model, priors=None):
+        """
+        Posterior for Poisson lightcurve data. Primary intended use is for
+        modelling X-ray light curves, but alternative uses are conceivable.
+
+        TODO: Include astropy.modeling models
+
+        Parameters
+        ----------
+        x : numpy.ndarray
+            The independent variable (e.g. time stamps of a light curve)
+
+        y : numpy.ndarray
+            The dependent variable (e.g. counts per bin of a light curve)
+
+        model: instance of any subclass of parameterclass.ParametricModel
+            The model for the power spectrum. Note that in order to define
+            the posterior properly, the ParametricModel subclass must be
+            instantiated with the hyperpars parameter set, or there won't
+            be a prior to be calculated! If all this object is used
+            for a maximum likelihood-style analysis, no prior is required.
+
+        priors : dict of form {"parameter name": function}, optional
+            A dictionary with the definitions for the prior probabilities.
+            For each parameter in `model`, there must be a prior defined with
+            a key of the exact same name as stored in `model.param_names`.
+            The item for each key is a function definition defining the prior
+            (e.g. a lambda function or a `scipy.stats.distribution.pdf`.
+            If `priors = None`, then no prior is set. This means priors need
+            to be added by hand using the `set_logprior` function defined in
+            this module. Note that it is impossible to call the posterior object
+            itself or the `self.logposterior` method without defining a prior.
+
+        Attributes
+        ----------
+        x: numpy.ndarray
+            The independent variable (list of frequencies) stored in ps.freq
+
+        y: numpy.ndarray
+            The dependent variable (list of powers) stored in ps.power
+
+        model: instance of any subclass of parameterclass.ParametricModel
+               The model for the power spectrum. Note that in order to define
+               the posterior properly, the ParametricModel subclass must be
+               instantiated with the hyperpars parameter set, or there won't
+               be a prior to be calculated! If all this object is used
+               for a maximum likelihood-style analysis, no prior is required.
+
+        """
+        self.x = x
+        self.y = y
+
+        self.loglikelihood = PoissonLogLikelihood(self.x, self.y, model)
+
+        Posterior.__init__(self, self.x, self.y, model)
+
+        if not priors is None:
+            self.logprior = set_logprior(self, priors)
+
+
+class GaussianPosterior(Posterior):
+
+    def __init__(self, x, y, yerr, model, priors=None):
+        """
+        A general class for two-dimensional data following a Gaussian
+        sampling distribution.
+
+        Parameters
+        ----------
+        x: numpy.ndarray
+            independent variable
+
+        y: numpy.ndarray
+            dependent variable
+
+        yerr: numpy.ndarray
+            measurement uncertainties for y
+
+        model: instance of any subclass of parameterclass.ParametricModel
+            The model for the power spectrum. Note that in order to define
+            the posterior properly, the ParametricModel subclass must be
+            instantiated with the hyperpars parameter set, or there won't
+            be a prior to be calculated! If all this object is used
+            for a maximum likelihood-style analysis, no prior is required.
+
+        """
+        self.loglikelihood = GaussianLogLikelihood(x, y, yerr, model)
+
+        Posterior.__init__(self, x, y, model)
+
+        self.yerr = yerr
+
+        if not priors is None:
+            self.logprior = set_logprior(self, priors)
diff --git a/stingray/modeling/scripts.py b/stingray/modeling/scripts.py
new file mode 100644
index 00000000..df91c9be
--- /dev/null
+++ b/stingray/modeling/scripts.py
@@ -0,0 +1,244 @@
+
+import numpy as np
+from astropy.modeling import models
+
+from stingray.modeling import PSDParEst, PSDPosterior, PSDLogLikelihood
+from stingray import Powerspectrum
+
+__all__ = ["fit_powerspectrum", "fit_lorentzians"]
+
+
+def fit_powerspectrum(ps, model, starting_pars, max_post=False, priors=None,
+              fitmethod="L-BFGS-B"):
+    """
+    Fit a number of Lorentzians to a power spectrum, possibly including white
+    noise. Each Lorentzian has three parameters (amplitude, centroid position,
+    full-width at half maximum), plus one extra parameter if the white noise
+    level should be fit as well. Priors for each parameter can be included in
+    case `max_post = True`, in which case the function will attempt a
+    Maximum-A-Posteriori fit. Priors must be specified as a dictionary with one
+    entry for each parameter.
+    The parameter names are `(amplitude_i, x_0_i, fwhm_i)` for each `i` out of
+    a total of `N` Lorentzians. The white noise level has a parameter
+    `amplitude_(N+1)`. For example, a model with two Lorentzians and a
+    white noise level would have parameters:
+    [amplitude_0, x_0_0, fwhm_0, amplitude_1, x_0_1, fwhm_1, amplitude_2].
+
+    Parameters
+    ----------
+    ps : Powerspectrum
+        A Powerspectrum object with the data to be fit
+
+    nlor : int
+        The number of Lorentzians to fit
+
+    starting_pars : iterable
+        The list of starting guesses for the optimizer. See explanation above
+        for ordering of parameters in this list.
+
+    fit_whitenoise : bool, optional, default True
+        If True, the code will attempt to fit a white noise level along with
+        the Lorentzians. Be sure to include a starting parameter for the
+        optimizer in `starting_pars`!
+
+    max_post : bool, optional, default False
+        If True, perform a Maximum-A-Posteriori fit of the data rather than a
+        Maximum Likelihood fit. Note that this requires priors to be specified,
+        otherwise this will cause an exception!
+
+    priors : {dict | None}, optional, default None
+        Dictionary with priors for the MAP fit. This should be of the form
+        {"parameter name": probability distribution, ...}
+
+    fitmethod : string, optional, default "L-BFGS-B"
+        Specifies an optimization algorithm to use. Supply any valid option for
+        `scipy.optimize.minimize`.
+
+    Returns
+    -------
+    parest : PSDParEst object
+        A PSDParEst object for further analysis
+
+    res : OptimizationResults object
+        The OptimizationResults object storing useful results and quantities
+        relating to the fit
+
+    Example
+    -------
+
+    We start by making an example power spectrum with three Lorentzians
+    >>> m = 1
+    >>> nfreq = 100000
+    >>> freq = np.linspace(1, 1000, nfreq)
+
+    >>> np.random.seed(100)  # set the seed for the random number generator
+    >>> noise = np.random.exponential(size=nfreq)
+
+    >>> model = models.PowerLaw1D() + models.Const1D()
+    >>> model.x_0_0.fixed = True
+
+    >>> alpha_0 = 2.0
+    >>> amplitude_0 = 100.0
+    >>> amplitude_1 = 2.0
+
+    >>> model.alpha_0 = alpha_0
+    >>> model.amplitude_0 = amplitude_0
+    >>> model.amplitude_1 = amplitude_1
+
+    >>> p = model(freq)
+    >>> power = noise * p
+
+    >>> ps = Powerspectrum()
+    >>> ps.freq = freq
+    >>> ps.power = power
+    >>> ps.m = m
+    >>> ps.df = freq[1] - freq[0]
+    >>> ps.norm = "leahy"
+
+    Now we have to guess starting parameters. For each Lorentzian, we have
+    amplitude, centroid position and fwhm, and this pattern repeats for each
+    Lorentzian in the fit. The white noise level is the last parameter.
+    >>> t0 = [80, 1.5, 2.5]
+
+    Let's also make a model to test:
+    >>> model_to_test = models.PowerLaw1D() + models.Const1D()
+    >>> model_to_test.x_0_0.fixed = True
+
+    We're ready for doing the fit:
+    >>> parest, res = fit_powerspectrum(ps, model_to_test, t0)
+
+    `res` contains a whole array of useful information about the fit, for
+    example the parameters at the optimum:
+    >>> p_opt = res.p_opt
+
+    """
+    if priors:
+        lpost = PSDPosterior(ps, model, priors=priors)
+    else:
+        lpost = PSDLogLikelihood(ps.freq, ps.power, model, m=ps.m)
+
+    parest = PSDParEst(ps, fitmethod=fitmethod, max_post=max_post)
+    res = parest.fit(lpost, starting_pars, neg=True)
+
+    return parest, res
+
+def fit_lorentzians(ps, nlor, starting_pars, fit_whitenoise=True,
+                    max_post=False, priors=None,
+                    fitmethod="L-BFGS-B"):
+    """
+    Fit a number of Lorentzians to a power spectrum, possibly including white
+    noise. Each Lorentzian has three parameters (amplitude, centroid position,
+    full-width at half maximum), plus one extra parameter if the white noise
+    level should be fit as well. Priors for each parameter can be included in
+    case `max_post = True`, in which case the function will attempt a
+    Maximum-A-Posteriori fit. Priors must be specified as a dictionary with one
+    entry for each parameter.
+    The parameter names are `(amplitude_i, x_0_i, fwhm_i)` for each `i` out of
+    a total of `N` Lorentzians. The white noise level has a parameter
+    `amplitude_(N+1)`. For example, a model with two Lorentzians and a
+    white noise level would have parameters:
+    [amplitude_0, x_0_0, fwhm_0, amplitude_1, x_0_1, fwhm_1, amplitude_2].
+
+    Parameters
+    ----------
+    ps : Powerspectrum
+        A Powerspectrum object with the data to be fit
+
+    nlor : int
+        The number of Lorentzians to fit
+
+    starting_pars : iterable
+        The list of starting guesses for the optimizer. See explanation above
+        for ordering of parameters in this list.
+
+    fit_whitenoise : bool, optional, default True
+        If True, the code will attempt to fit a white noise level along with
+        the Lorentzians. Be sure to include a starting parameter for the
+        optimizer in `starting_pars`!
+
+    max_post : bool, optional, default False
+        If True, perform a Maximum-A-Posteriori fit of the data rather than a
+        Maximum Likelihood fit. Note that this requires priors to be specified,
+        otherwise this will cause an exception!
+
+    priors : {dict | None}, optional, default None
+        Dictionary with priors for the MAP fit. This should be of the form
+        {"parameter name": probability distribution, ...}
+
+    fitmethod : string, optional, default "L-BFGS-B"
+        Specifies an optimization algorithm to use. Supply any valid option for
+        `scipy.optimize.minimize`.
+
+    Returns
+    -------
+    parest : PSDParEst object
+        A PSDParEst object for further analysis
+
+    res : OptimizationResults object
+        The OptimizationResults object storing useful results and quantities
+        relating to the fit
+
+    Example
+    -------
+
+    We start by making an example power spectrum with three Lorentzians
+    >>> np.random.seed(400)
+    >>> nlor = 3
+
+    >>> x_0_0 = 0.5
+    >>> x_0_1 = 2.0
+    >>> x_0_2 = 7.5
+
+    >>> amplitude_0 = 150.0
+    >>> amplitude_1 = 50.0
+    >>> amplitude_2 = 15.0
+
+    >>> fwhm_0 = 0.1
+    >>> fwhm_1 = 1.0
+    >>> fwhm_2 = 0.5
+
+    We will also include a white noise level:
+    >>> whitenoise = 2.0
+
+    >>> model = models.Lorentz1D(amplitude_0, x_0_0, fwhm_0) + \\
+    ...         models.Lorentz1D(amplitude_1, x_0_1, fwhm_1) + \\
+    ...         models.Lorentz1D(amplitude_2, x_0_2, fwhm_2) + \\
+    ...         models.Const1D(whitenoise)
+
+    >>> freq = np.linspace(0.01, 10.0, 10.0/0.01)
+    >>> p = model(freq)
+    >>> noise = np.random.exponential(size=len(freq))
+
+    >>> power = p*noise
+    >>> ps = Powerspectrum()
+    >>> ps.freq = freq
+    >>> ps.power = power
+    >>> ps.df = ps.freq[1] - ps.freq[0]
+    >>> ps.m = 1
+
+    Now we have to guess starting parameters. For each Lorentzian, we have
+    amplitude, centroid position and fwhm, and this pattern repeats for each
+    Lorentzian in the fit. The white noise level is the last parameter.
+    >>> t0 = [150, 0.4, 0.2, 50, 2.3, 0.6, 20, 8.0, 0.4, 2.1]
+
+    We're ready for doing the fit:
+    >>> parest, res = fit_lorentzians(ps, nlor, t0)
+
+    `res` contains a whole array of useful information about the fit, for
+    example the parameters at the optimum:
+    >>> p_opt = res.p_opt
+
+    """
+
+    model = models.Lorentz1D()
+
+    if nlor > 1:
+        for i in range( nlor -1):
+            model += models.Lorentz1D()
+
+    if fit_whitenoise:
+        model += models.Const1D()
+
+    return fit_powerspectrum(ps, model, starting_pars, max_post=max_post,
+                             priors=priors, fitmethod=fitmethod)
+
diff --git a/stingray/scraps.py b/stingray/scraps.py
deleted file mode 100644
index 1603042b..00000000
--- a/stingray/scraps.py
+++ /dev/null
@@ -1,1 +0,0 @@
-__author__ = 'danielahuppenkothen'
diff --git a/stingray/simulator/simulator.py b/stingray/simulator/simulator.py
index ae7f0258..3264180f 100644
--- a/stingray/simulator/simulator.py
+++ b/stingray/simulator/simulator.py
@@ -255,10 +255,10 @@ class Simulator(object):
         """
 
         # Fill in 0 entries until the start time
-        h_zeros = np.zeros(start/self.dt)
+        h_zeros = np.zeros(int(start/self.dt))
 
         # Define constant impulse response
-        h_ones = np.ones(width/self.dt) * intensity
+        h_ones = np.ones(int(width/self.dt)) * intensity
 
         return np.append(h_zeros, h_ones)
 
diff --git a/stingray/simulator/transfer.py b/stingray/simulator/transfer.py
index f8cc3f80..5b543049 100644
--- a/stingray/simulator/transfer.py
+++ b/stingray/simulator/transfer.py
@@ -256,10 +256,10 @@ def simple_ir(dt=0.125, start=0, width=1000, intensity=1):
     """
 
     # Fill in 0 entries until the start time
-    h_zeros = np.zeros(start/dt)
+    h_zeros = np.zeros(int(start/dt))
 
     # Define constant impulse response
-    h_ones = np.ones(width/dt) * intensity
+    h_ones = np.ones(int(width/dt)) * intensity
 
     return np.append(h_zeros, h_ones)
 


Test patch:
diff --git a/stingray/modeling/tests/__init__.py b/stingray/modeling/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/stingray/modeling/tests/test_parameterestimation.py b/stingray/modeling/tests/test_parameterestimation.py
new file mode 100644
index 00000000..7717ac52
--- /dev/null
+++ b/stingray/modeling/tests/test_parameterestimation.py
@@ -0,0 +1,768 @@
+import numpy as np
+import scipy.stats
+import os
+
+from astropy.tests.helper import pytest
+from astropy.modeling import models
+from astropy.modeling.fitting import _fitter_to_model_params
+
+from stingray import Powerspectrum
+
+from stingray.modeling import ParameterEstimation, PSDParEst, \
+    OptimizationResults, SamplingResults
+from stingray.modeling import PSDPosterior, set_logprior, PSDLogLikelihood
+
+try:
+    from statsmodels.tools.numdiff import approx_hess
+    comp_hessian = True
+except ImportError:
+    comp_hessian = False
+
+try:
+    import emcee
+    can_sample = True
+except ImportError:
+    can_sample = False
+
+
+class TestParameterEstimation(object):
+
+    @classmethod
+    def setup_class(cls):
+        m = 1
+        nfreq = 100000
+        freq = np.arange(nfreq)
+        noise = np.random.exponential(size=nfreq)
+        power = noise * 2.0
+
+        ps = Powerspectrum()
+        ps.freq = freq
+        ps.power = power
+        ps.m = m
+        ps.df = freq[1] - freq[0]
+        ps.norm = "leahy"
+
+        cls.ps = ps
+        cls.a_mean, cls.a_var = 2.0, 1.0
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(amplitude)
+
+        cls.priors = {"amplitude": p_amplitude}
+        cls.lpost = PSDPosterior(cls.ps, cls.model)
+        cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
+
+    def test_par_est_initializes(self):
+        pe = ParameterEstimation()
+
+    def test_parest_stores_max_post_correctly(self):
+        """
+        Make sure the keyword for Maximum A Posteriori fits is stored correctly
+        as a default.
+        """
+        pe = ParameterEstimation()
+
+        assert pe.max_post is True, "max_post should be set to True as a default."
+
+    def test_object_works_with_loglikelihood_object(self):
+        llike = PSDLogLikelihood(self.ps.freq, self.ps.power, self.model)
+        pe = ParameterEstimation()
+        res = pe.fit(llike, [2.0])
+
+        pass
+
+    def test_fit_fails_when_object_is_not_posterior_or_likelihood(self):
+        x = np.ones(10)
+        y = np.ones(10)
+        pe = ParameterEstimation()
+        with pytest.raises(TypeError):
+            res = pe.fit(x, y)
+
+    def test_fit_fails_without_lpost_or_t0(self):
+        pe = ParameterEstimation()
+        with pytest.raises(TypeError):
+            res = pe.fit()
+
+    def test_fit_fails_without_t0(self):
+        pe = ParameterEstimation()
+        with pytest.raises(TypeError):
+            res = pe.fit(np.ones(10))
+
+    def test_fit_fails_with_incorrect_number_of_parameters(self):
+        pe = ParameterEstimation()
+        t0 = [1, 2]
+        with pytest.raises(ValueError):
+            res = pe.fit(self.lpost, t0)
+
+    def test_fit_method_works_with_correct_parameter(self):
+        pe = ParameterEstimation()
+        t0 = [2.0]
+        res = pe.fit(self.lpost, t0)
+
+    def test_fit_method_returns_optimization_results_object(self):
+        pe = ParameterEstimation()
+        t0 = [2.0]
+        res = pe.fit(self.lpost, t0)
+        assert isinstance(res,
+                          OptimizationResults), "res must be of type OptimizationResults"
+
+    def test_compute_lrt_fails_when_garbage_goes_in(self):
+        pe = ParameterEstimation()
+        t0 = [2.0]
+
+        with pytest.raises(TypeError):
+            pe.compute_lrt(self.lpost, t0, None, t0)
+
+        with pytest.raises(ValueError):
+            pe.compute_lrt(self.lpost, t0[:-1], self.lpost, t0)
+
+    def test_compute_lrt_sets_max_post_to_false(self):
+        t0 = [2.0]
+        pe = ParameterEstimation(max_post=True)
+
+        assert pe.max_post is True
+        delta_deviance = pe.compute_lrt(self.lpost, t0, self.lpost, t0)
+
+        assert pe.max_post is False
+
+    def test_compute_lrt_computes_deviance_correctly(self):
+        t0 = [2.0]
+        pe = ParameterEstimation()
+
+        delta_deviance = pe.compute_lrt(self.lpost, t0, self.lpost, t0)
+
+        assert delta_deviance < 1e-7
+
+    def test_sampler_runs(self):
+        pe = ParameterEstimation()
+        if os.path.exists("test_corner.pdf"):
+            os.unlink("test_corner.pdf")
+        sample_res = pe.sample(self.lpost, [2.0], nwalkers=100, niter=10,
+                               burnin=50, print_results=True, plot=True)
+
+        assert os.path.exists("test_corner.pdf")
+        assert sample_res.acceptance > 0.25
+        assert isinstance(sample_res, SamplingResults)
+
+
+class TestOptimizationResults(object):
+
+    @classmethod
+    def setup_class(cls):
+        m = 1
+        nfreq = 100000
+        freq = np.arange(nfreq)
+        noise = np.random.exponential(size=nfreq)
+        power = noise * 2.0
+
+        ps = Powerspectrum()
+        ps.freq = freq
+        ps.power = power
+        ps.m = m
+        ps.df = freq[1] - freq[0]
+        ps.norm = "leahy"
+
+        cls.ps = ps
+        cls.a_mean, cls.a_var = 2.0, 1.0
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(amplitude)
+
+        cls.priors = {"amplitude": p_amplitude}
+        cls.lpost = PSDPosterior(cls.ps, cls.model)
+        cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
+
+        cls.fitmethod = "BFGS"
+        cls.max_post = True
+        cls.t0 = [2.0]
+        cls.neg = True
+        cls.opt = scipy.optimize.minimize(cls.lpost, cls.t0,
+                                          method=cls.fitmethod,
+                                          args=cls.neg, tol=1.e-10)
+
+    def test_object_initializes(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+
+    def test_object_has_right_attributes(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+
+        assert hasattr(res, "p_opt")
+        assert hasattr(res, "result")
+        assert hasattr(res, "deviance")
+        assert hasattr(res, "aic")
+        assert hasattr(res, "bic")
+        assert hasattr(res, "model")
+        assert isinstance(res.model, models.Const1D)
+
+    def test_p_opt_is_correct(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+        assert res.p_opt == self.opt.x, "res.p_opt must be the same as opt.x!"
+        assert np.isclose(res.p_opt[0], 2.0, atol=0.1, rtol=0.1)
+
+    def test_model_is_same_as_in_lpost(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+        assert res.model == self.lpost.model
+
+    def test_result_is_same_as_in_opt(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+        assert res.result == self.opt.fun
+
+    def test_compute_model_works_correctly(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+        mean_model = np.ones_like(self.lpost.x) * self.opt.x[0]
+        assert np.all(
+            res.mfit == mean_model), "res.model should be exactly " \
+                                     "the model for the data."
+
+    def test_compute_criteria_computes_criteria_correctly(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+
+        test_aic = 169440.83719024697
+        test_bic = 169245.62163088709
+        test_deviance = 337950.6823459795
+
+        assert np.isclose(res.aic, test_aic, atol=0.1, rtol=0.1)
+        assert np.isclose(res.bic, test_bic, atol=0.1, rtol=0.1)
+        assert np.isclose(res.deviance, test_deviance, atol=0.1, rtol=0.1)
+
+    def test_merit_calculated_correctly(self):
+        res = OptimizationResults(self.lpost, self.opt, neg=self.neg)
+
+        test_merit = 98770.654981073574
+        assert np.isclose(res.merit, test_merit, atol=0.1, rtol=0.1)
+
+    def test_res_is_of_correct_type(self):
+        pe = ParameterEstimation()
+        t0 = [2.0]
+        res = pe.fit(self.lpost, t0)
+
+        assert isinstance(res, OptimizationResults)
+
+
+class OptimizationResultsSubclassDummy(OptimizationResults):
+
+    def __init__(self, lpost, res, neg):
+        self.neg = neg
+        self.result = res.fun
+        self.p_opt = res.x
+        self.model = lpost.model
+
+
+class TestOptimizationResultInternalFunctions(object):
+
+    @classmethod
+    def setup_class(cls):
+        m = 1
+        nfreq = 100000
+        freq = np.linspace(1, 1000, nfreq)
+
+        np.random.seed(100)  # set the seed for the random number generator
+        noise = np.random.exponential(size=nfreq)
+
+        cls.model = models.PowerLaw1D() + models.Const1D()
+        cls.model.x_0_0.fixed = True
+
+        cls.alpha_0 = 2.0
+        cls.amplitude_0 = 100.0
+        cls.amplitude_1 = 2.0
+
+        cls.model.alpha_0 = cls.alpha_0
+        cls.model.amplitude_0 = cls.amplitude_0
+        cls.model.amplitude_1 = cls.amplitude_1
+
+        p = cls.model(freq)
+        power = noise * p
+
+        ps = Powerspectrum()
+        ps.freq = freq
+        ps.power = power
+        ps.m = m
+        ps.df = freq[1] - freq[0]
+        ps.norm = "leahy"
+
+        cls.ps = ps
+        cls.a_mean, cls.a_var = 2.0, 1.0
+        cls.a2_mean, cls.a2_var = 100.0, 10.0
+
+        p_amplitude_1 = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(amplitude)
+
+        p_alpha_0 = lambda alpha: \
+            scipy.stats.uniform(0.0, 5.0).pdf(alpha)
+
+        p_amplitude_0 = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a2_mean, scale=cls.a2_var).pdf(
+                amplitude)
+
+        cls.priors = {"amplitude_1": p_amplitude_1,
+                      "amplitude_0": p_amplitude_0,
+                      "alpha_0": p_alpha_0}
+
+        cls.lpost = PSDPosterior(cls.ps, cls.model)
+        cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
+
+        cls.fitmethod = "BFGS"
+        cls.max_post = True
+        cls.t0 = [cls.amplitude_0, cls.alpha_0, cls.amplitude_1]
+        cls.neg = True
+        cls.opt = scipy.optimize.minimize(cls.lpost, cls.t0,
+                                          method=cls.fitmethod,
+                                          args=cls.neg, tol=1.e-5)
+
+        cls.optres = OptimizationResultsSubclassDummy(cls.lpost, cls.opt,
+                                                      neg=True)
+
+    def test_output(self, capsys):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_covariance(self.lpost, self.opt)
+        optres.print_summary(self.lpost)
+        out, err = capsys.readouterr()
+        assert "The best-fit model parameter" in out
+        assert "merit function (SSE)" in out
+
+    def test_compute_model(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_model(self.lpost)
+
+        assert hasattr(optres,
+                       "mfit"), "OptimizationResult object should have mfit " \
+                                "attribute at this point!"
+
+        _fitter_to_model_params(self.model, self.opt.x)
+        mfit_test = self.model(self.lpost.x)
+
+        assert np.all(optres.mfit == mfit_test)
+
+    def test_compute_statistics_computes_mfit(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        assert hasattr(optres, "mfit") is False
+        optres._compute_statistics(self.lpost)
+
+        assert hasattr(optres, "mfit")
+
+    def test_compute_statistics_computes_all_statistics(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_statistics(self.lpost)
+
+        assert hasattr(optres, "merit")
+        assert hasattr(optres, "dof")
+        assert hasattr(optres, "sexp")
+        assert hasattr(optres, "ssd")
+        assert hasattr(optres, "sobs")
+
+    def test_compute_statistics_returns_correct_values(self):
+        test_merit = 99765.718448514497
+        test_dof = 99997.0
+        test_sexp = 600000.0
+        test_ssd = 1095.4451150103323
+        test_sobs = -154.901207861497
+
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_statistics(self.lpost)
+
+        assert np.isclose(test_merit, optres.merit, atol=0.01, rtol=0.01)
+        assert test_dof == optres.dof
+        assert test_sexp == optres.sexp
+        assert test_ssd == optres.ssd
+        assert np.isclose(test_sobs, optres.sobs, atol=0.01, rtol=0.01)
+
+    def test_compute_criteria_returns_correct_attributes(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_criteria(self.lpost)
+
+        assert hasattr(optres, "aic")
+        assert hasattr(optres, "bic")
+        assert hasattr(optres, "deviance")
+
+    def test_compute_criteria_returns_correct_values(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_criteria(self.lpost)
+
+        test_aic = 170988.9174964963
+        test_bic = 171017.4562728912
+        test_deviance = 341954.16168676887
+
+        assert np.isclose(test_aic, optres.aic)
+        assert np.isclose(test_bic, optres.bic)
+        assert np.isclose(test_deviance, optres.deviance)
+
+    def test_compute_covariance_with_hess_inverse(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
+                                                  neg=True)
+
+        optres._compute_covariance(self.lpost, self.opt)
+
+        assert np.all(optres.cov == np.asarray(self.opt.hess_inv))
+        assert np.all(optres.err == np.sqrt(np.diag(self.opt.hess_inv)))
+
+    def test_compute_covariance_without_hess_inverse(self):
+        fitmethod = "powell"
+        opt = scipy.optimize.minimize(self.lpost, self.t0,
+                                      method=fitmethod,
+                                      args=self.neg, tol=1.e-10)
+
+        optres = OptimizationResultsSubclassDummy(self.lpost, opt,
+                                                  neg=True)
+
+        optres._compute_covariance(self.lpost, opt)
+
+        if comp_hessian:
+            phess = approx_hess(opt.x, self.lpost)
+            hess_inv = np.linalg.inv(phess)
+
+            assert np.all(optres.cov == hess_inv)
+            assert np.all(optres.err == np.sqrt(np.diag(np.abs(hess_inv))))
+
+
+if can_sample:
+    class SamplingResultsDummy(SamplingResults):
+        def __init__(self, sampler, ci_min=0.05, ci_max=0.95):
+            # store all the samples
+            self.samples = sampler.flatchain
+
+            self.nwalkers = np.float(sampler.chain.shape[0])
+            self.niter = np.float(sampler.iterations)
+
+            # store number of dimensions
+            self.ndim = sampler.dim
+
+            # compute and store acceptance fraction
+            self.acceptance = np.nanmean(sampler.acceptance_fraction)
+            self.L = self.acceptance * self.samples.shape[0]
+
+
+    class TestSamplingResults(object):
+        @classmethod
+        def setup_class(cls):
+            m = 1
+            nfreq = 100000
+            freq = np.arange(nfreq)
+            noise = np.random.exponential(size=nfreq)
+            power = noise * 2.0
+
+            ps = Powerspectrum()
+            ps.freq = freq
+            ps.power = power
+            ps.m = m
+            ps.df = freq[1] - freq[0]
+            ps.norm = "leahy"
+
+            cls.ps = ps
+            cls.a_mean, cls.a_var = 2.0, 1.0
+
+            cls.model = models.Const1D()
+
+            p_amplitude = lambda amplitude: \
+                scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(
+                    amplitude)
+
+            cls.priors = {"amplitude": p_amplitude}
+            cls.lpost = PSDPosterior(cls.ps, cls.model)
+            cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
+
+            cls.fitmethod = "BFGS"
+            cls.max_post = True
+            cls.t0 = [2.0]
+            cls.neg = True
+
+            pe = ParameterEstimation()
+            res = pe.fit(cls.lpost, cls.t0)
+
+            cls.nwalkers = 100
+            cls.niter = 200
+
+            np.random.seed(200)
+            p0 = np.array(
+                [np.random.multivariate_normal(res.p_opt, res.cov) for
+                 i in range(cls.nwalkers)])
+
+            cls.sampler = emcee.EnsembleSampler(cls.nwalkers,
+                                                len(res.p_opt), cls.lpost,
+                                                args=[False], threads=1)
+
+            _, _, _ = cls.sampler.run_mcmc(p0, cls.niter)
+
+        def test_sample_results_object_initializes(self):
+            SamplingResults(self.sampler)
+
+        def test_sample_results_produces_attributes(self):
+            s = SamplingResults(self.sampler)
+
+            assert s.samples.shape[0] == self.nwalkers * self.niter
+
+        def test_sampling_results_acceptance_ratio(self):
+            s = SamplingResults(self.sampler)
+
+            assert s.acceptance > 0.25
+            assert np.isclose(s.L,
+                              s.acceptance * self.nwalkers * self.niter)
+
+        def test_check_convergence_works(self):
+            s = SamplingResultsDummy(self.sampler)
+            s._check_convergence(self.sampler)
+
+            assert hasattr(s, "rhat")
+
+        def test_rhat_computes_correct_answer(self):
+            s = SamplingResults(self.sampler)
+
+            rhat_test = 3.81886815e-06
+
+            assert np.isclose(rhat_test, s.rhat[0], atol=0.001, rtol=0.001)
+
+        def test_infer_works(self):
+            s = SamplingResultsDummy(self.sampler)
+            s._infer()
+
+            assert hasattr(s, "mean")
+            assert hasattr(s, "std")
+            assert hasattr(s, "ci")
+
+        def test_infer_computes_correct_values(self):
+            s = SamplingResults(self.sampler)
+
+            test_mean = 2.00190793
+            test_std = 0.00195719
+            test_ci = [[1.99435539], [1.9971502]]
+
+            assert np.isclose(test_mean, s.mean[0], atol=0.01, rtol=0.01)
+            assert np.isclose(test_std, s.std[0], atol=0.01, rtol=0.01)
+            assert np.all(np.isclose(test_ci, s.ci, atol=0.01, rtol=0.01))
+
+
+class TestPSDParEst(object):
+
+    @classmethod
+    def setup_class(cls):
+
+        m = 1
+        nfreq = 100000
+        freq = np.linspace(1, 10.0, nfreq)
+
+        np.random.seed(100) # set the seed for the random number generator
+        noise = np.random.exponential(size=nfreq)
+
+        cls.model = models.Lorentz1D() + models.Const1D()
+
+        cls.x_0_0 = 2.0
+        cls.fwhm_0 = 0.1
+        cls.amplitude_0 = 100.0
+
+        cls.amplitude_1 = 2.0
+
+        cls.model.x_0_0 = cls.x_0_0
+        cls.model.fwhm_0 = cls.fwhm_0
+        cls.model.amplitude_0 = cls.amplitude_0
+        cls.model.amplitude_1 = cls.amplitude_1
+
+        p = cls.model(freq)
+
+        np.random.seed(400)
+        power = noise*p
+
+        ps = Powerspectrum()
+        ps.freq = freq
+        ps.power = power
+        ps.m = m
+        ps.df = freq[1]-freq[0]
+        ps.norm = "leahy"
+
+        cls.ps = ps
+        cls.a_mean, cls.a_var = 2.0, 1.0
+        cls.a2_mean, cls.a2_var = 100.0, 10.0
+
+        p_amplitude_1 = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(amplitude)
+
+        p_x_0_0 = lambda alpha: \
+            scipy.stats.uniform(0.0, 5.0).pdf(alpha)
+
+        p_fwhm_0 = lambda alpha: \
+            scipy.stats.uniform(0.0, 0.5).pdf(alpha)
+
+        p_amplitude_0 = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a2_mean, scale=cls.a2_var).pdf(amplitude)
+
+        cls.priors = {"amplitude_1": p_amplitude_1,
+                      "amplitude_0": p_amplitude_0,
+                      "x_0_0": p_x_0_0,
+                      "fwhm_0": p_fwhm_0}
+
+        cls.lpost = PSDPosterior(cls.ps, cls.model)
+        cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
+
+        cls.fitmethod = "BFGS"
+        cls.max_post = True
+        cls.t0 = [cls.x_0_0, cls.fwhm_0, cls.amplitude_0, cls.amplitude_1]
+        cls.neg = True
+
+
+    def test_par_est_initializes(self):
+        pe = PSDParEst(self.ps)
+
+    def test_parest_stores_max_post_correctly(self):
+        """
+        Make sure the keyword for Maximum A Posteriori fits is stored correctly
+        as a default.
+        """
+        pe = PSDParEst(self.ps)
+
+        assert pe.max_post is True, "max_post should be set to True as a default."
+
+    def test_fit_fails_when_object_is_not_posterior_or_likelihood(self):
+        x = np.ones(10)
+        y = np.ones(10)
+        pe = PSDParEst(self.ps)
+        with pytest.raises(TypeError):
+            res = pe.fit(x, y)
+
+    def test_fit_fails_without_lpost_or_t0(self):
+        pe = PSDParEst(self.ps)
+        with pytest.raises(TypeError):
+            res = pe.fit()
+
+    def test_fit_fails_without_t0(self):
+        pe = PSDParEst(self.ps)
+        with pytest.raises(TypeError):
+            res = pe.fit(np.ones(10))
+
+    def test_fit_fails_with_incorrect_number_of_parameters(self):
+        pe = PSDParEst(self.ps)
+        t0 = [1,2]
+        with pytest.raises(ValueError):
+            res = pe.fit(self.lpost, t0)
+
+    def test_fit_method_works_with_correct_parameter(self):
+        pe = PSDParEst(self.ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+        t0 = [2.0, 1, 1, 1]
+        res = pe.fit(lpost, t0)
+
+    def test_fit_method_returns_optimization_results_object(self):
+        pe = PSDParEst(self.ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+        t0 = [2.0, 1, 1, 1]
+        res = pe.fit(self.lpost, t0)
+        assert isinstance(res, OptimizationResults), "res must be of type " \
+                                                     "OptimizationResults"
+
+    def test_plotfits(self):
+        pe = PSDParEst(self.ps)
+        t0 = [2.0, 1, 1, 1]
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        res = pe.fit(lpost, t0)
+
+        pe.plotfits(res, save_plot=True)
+
+        assert os.path.exists("test_ps_fit.png")
+        os.unlink("test_ps_fit.png")
+
+    def test_plotfits_log(self):
+        pe = PSDParEst(self.ps)
+        t0 = [2.0, 1, 1, 1]
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        res = pe.fit(self.lpost, t0)
+
+        pe.plotfits(res, save_plot=True, log=True)
+
+        assert os.path.exists("test_ps_fit.png")
+        os.unlink("test_ps_fit.png")
+
+    def test_plotfits2(self):
+        t0 = [2.0, 1, 1, 1]
+        ps = Powerspectrum()
+        ps.freq = self.ps.freq
+        ps.power = self.ps.power
+        ps.m = self.ps.m
+        ps.df = self.ps.df
+        ps.norm = "rms"
+        pe = PSDParEst(ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        res = pe.fit(self.lpost, t0)
+
+        pe.plotfits(res, res2=res, save_plot=True)
+
+        assert os.path.exists("test_ps_fit.png")
+        os.unlink("test_ps_fit.png")
+
+    def test_plotfits_log2(self):
+        ps = Powerspectrum()
+        ps.freq = self.ps.freq
+        ps.power = self.ps.power
+        ps.m = self.ps.m
+        ps.df = self.ps.df
+        ps.norm = "rms"
+        pe = PSDParEst(ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        t0 = [2.0, 1, 1, 1]
+        res = pe.fit(self.lpost, t0)
+
+        pe.plotfits(res, res2=res, save_plot=True, log=True)
+
+        assert os.path.exists("test_ps_fit.png")
+        os.unlink("test_ps_fit.png")
+
+    def test_compute_lrt_fails_when_garbage_goes_in(self):
+        pe = PSDParEst(self.ps)
+        t0 = [2.0, 1, 1, 1]
+
+        with pytest.raises(TypeError):
+            pe.compute_lrt(self.lpost, t0, None, t0)
+
+        with pytest.raises(ValueError):
+            pe.compute_lrt(self.lpost, t0[:-1], self.lpost, t0)
+
+    def test_compute_lrt_sets_max_post_to_false(self):
+        t0 = [2.0, 1, 1, 1]
+        pe = PSDParEst(self.ps, max_post=True)
+
+        assert pe.max_post is True
+        delta_deviance = pe.compute_lrt(self.lpost, t0, self.lpost, t0)
+
+        assert pe.max_post is False
+
+    def test_compute_lrt_computes_deviance_correctly(self):
+
+        t0 = [2.0, 1, 1, 1]
+        pe = PSDParEst(self.ps, max_post=True)
+
+        # MB: This is a little too random
+        delta_deviance = pe.compute_lrt(self.lpost, t0,
+                                        self.lpost, t0)
+
+        assert np.absolute(delta_deviance) < 1.5e-4
+
+    def test_sampler_runs(self):
+
+        pe = PSDParEst(self.ps)
+        lpost = PSDPosterior(self.ps, self.model, self.priors)
+
+        sample_res = pe.sample(lpost, [2.0, 0.1, 100, 2.0], nwalkers=50,
+                               niter=10, burnin=15, print_results=True,
+                               plot=True)
+        assert os.path.exists("test_corner.pdf")
+        os.unlink("test_corner.pdf")
+        assert sample_res.acceptance > 0.25
+        assert isinstance(sample_res, SamplingResults)
diff --git a/stingray/modeling/tests/test_posterior.py b/stingray/modeling/tests/test_posterior.py
new file mode 100644
index 00000000..8ce2830b
--- /dev/null
+++ b/stingray/modeling/tests/test_posterior.py
@@ -0,0 +1,613 @@
+import numpy as np
+import scipy.stats
+import copy
+
+from astropy.tests.helper import pytest
+from astropy.modeling import models
+from scipy.special import gammaln as scipy_gammaln
+
+from stingray import Lightcurve, Powerspectrum
+from stingray.modeling import Posterior, PSDPosterior, PoissonPosterior, GaussianPosterior
+from stingray.modeling import set_logprior
+from stingray.modeling.posterior import logmin
+from stingray.modeling.posterior import IncorrectParameterError
+
+np.random.seed(20150907)
+
+class TestSetPrior(object):
+
+    @classmethod
+    def setup_class(cls):
+        photon_arrivals = np.sort(np.random.uniform(0,1000, size=10000))
+        cls.lc = Lightcurve.make_lightcurve(photon_arrivals, dt=1.0)
+        cls.ps = Powerspectrum(cls.lc, norm="frac")
+        pl = models.PowerLaw1D()
+        pl.x_0.fixed = True
+
+        cls.lpost = PSDPosterior(cls.ps, pl)
+
+    def test_set_prior_runs(self):
+        p_alpha = lambda alpha: ((-1. <= alpha) & (alpha <= 5.))/6.0
+        p_amplitude = lambda amplitude: ((-10 <= np.log(amplitude)) &
+                                         ((np.log(amplitude) <= 10.0)))/20.0
+
+        priors = {"alpha":p_alpha, "amplitude":p_amplitude}
+        self.lpost.logprior = set_logprior(self.lpost, priors)
+
+    def test_prior_executes_correctly(self):
+        p_alpha = lambda alpha: ((-1. <= alpha) & (alpha <= 5.))/6.0
+        p_amplitude = lambda amplitude: ((-10 <= np.log(amplitude)) &
+                                         ((np.log(amplitude) <= 10.0)))/20.0
+
+        priors = {"alpha":p_alpha, "amplitude":p_amplitude}
+        self.lpost.logprior = set_logprior(self.lpost, priors)
+        true_logprior = np.log(1./6.) + np.log(1./20.0)
+        assert self.lpost.logprior([np.exp(0.0), np.exp(0.0)]) ==  true_logprior
+
+    def test_prior_returns_logmin_outside_prior_range(self):
+        p_alpha = lambda alpha: ((-1. <= alpha) & (alpha <= 5.))/6.0
+        p_amplitude = lambda amplitude: ((-10 <= np.log(amplitude)) &
+                                         ((np.log(amplitude) <= 10.0)))/20.0
+
+        priors = {"alpha":p_alpha, "amplitude":p_amplitude}
+        self.lpost.logprior = set_logprior(self.lpost, priors)
+        assert self.lpost.logprior([-2.0, np.exp(11.0)]) ==  logmin
+
+
+class PosteriorClassDummy(Posterior):
+    """
+    This is a test class that tests the basic functionality of the
+    Posterior superclass.
+    """
+    def __init__(self, x, y, model):
+        Posterior.__init__(self, x, y, model)
+
+    def  loglikelihood(self, t0, neg=False):
+        loglike = 1.0
+        return loglike
+
+    def logprior(self, t0):
+        lp = 2.0
+        return lp
+
+class TestPosterior(object):
+
+    @classmethod
+    def setup_class(cls):
+        cls.x = np.arange(100)
+        cls.y = np.ones(cls.x.shape[0])
+        cls.model = models.Const1D()
+        cls.p = PosteriorClassDummy(cls.x, cls.y, cls.model)
+        p_alpha = lambda alpha: ((-1. <= alpha) & (alpha <= 5.))/6.0
+
+        priors = {"amplitude":p_alpha}
+        cls.p.logprior = set_logprior(cls.p, priors)
+
+    def test_inputs(self):
+        assert np.allclose(self.p.x, self.x)
+        assert np.allclose(self.p.y, self.y)
+        assert isinstance(self.p.model, models.Const1D)
+
+    def test_call_method_positive(self):
+        t0 = [1]
+        post = self.p(t0, neg=False)
+        assert post == 1.0 + np.log(1./6.0)
+
+    def test_call_method_negative(self):
+        t0 = [1]
+        post = self.p(t0, neg=True)
+        assert post == -(1.0 + np.log(1./6.0))
+
+
+class TestPSDPosterior(object):
+
+    @classmethod
+    def setup_class(cls):
+        m = 1
+        nfreq = 1000000
+        freq = np.arange(nfreq)
+        noise = np.random.exponential(size=nfreq)
+        power = noise*2.0
+
+        ps = Powerspectrum()
+        ps.freq = freq
+        ps.power = power
+        ps.m = m
+        ps.df = freq[1]-freq[0]
+        ps.norm = "leahy"
+
+        cls.ps = ps
+        cls.a_mean, cls.a_var = 2.0, 1.0
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(amplitude)
+
+        cls.priors = {"amplitude":p_amplitude}
+
+    def test_logprior_fails_without_prior(self):
+        lpost = PSDPosterior(self.ps, self.model)
+
+        with pytest.raises(AttributeError):
+            lpost.logprior([1])
+
+    def test_making_posterior(self):
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert lpost.x.all() == self.ps.freq.all()
+        assert lpost.y.all() == self.ps.power.all()
+
+    def test_correct_number_of_parameters(self):
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        with pytest.raises(IncorrectParameterError):
+            lpost([2,3])
+
+    def test_logprior(self):
+        t0 = [2.0]
+
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        lp_test = lpost.logprior(t0)
+        lp = np.log(scipy.stats.norm(2.0, 1.0).pdf(t0))
+        assert lp == lp_test
+
+    def test_loglikelihood(self):
+        t0 = [2.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.ps.freq)
+
+        loglike = -np.sum(np.log(mean_model)) - np.sum(self.ps.power/mean_model)
+
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=False)
+
+        assert np.isclose(loglike, loglike_test)
+
+
+    def test_negative_loglikelihood(self):
+        t0 = [2.0]
+        m = self.model(self.ps.freq[1:], t0)
+        loglike = np.sum(self.ps.power[1:]/m + np.log(m))
+
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=True)
+
+        assert np.isclose(loglike, loglike_test)
+
+    def test_posterior(self):
+        t0 = [2.0]
+        m = self.model(self.ps.freq[1:], t0)
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=False)
+
+        loglike = -np.sum(self.ps.power[1:]/m + np.log(m))
+        logprior = np.log(scipy.stats.norm(2.0, 1.0).pdf(t0))
+        post = loglike + logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_negative_posterior(self):
+        t0 = [2.0]
+        m = self.model(self.ps.freq[1:], t0)
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=True)
+
+        loglike = -np.sum(self.ps.power[1:]/m + np.log(m))
+        logprior = np.log(scipy.stats.norm(2.0, 1.0).pdf(t0))
+        post = -loglike - logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+
+class TestPoissonPosterior(object):
+
+    @classmethod
+    def setup_class(cls):
+
+        nx = 1000000
+        cls.x = np.arange(nx)
+        cls.countrate = 10.0
+        cls.y = np.random.poisson(cls.countrate, size=cls.x.shape[0])
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.countrate, scale=cls.countrate).pdf(amplitude)
+
+        cls.priors = {"amplitude":p_amplitude}
+
+    def test_logprior_fails_without_prior(self):
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+
+        with pytest.raises(AttributeError):
+            lpost.logprior([10])
+
+    def test_making_posterior(self):
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert lpost.x.all() == self.x.all()
+        assert lpost.y.all() == self.y.all()
+
+    def test_correct_number_of_parameters(self):
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        with pytest.raises(IncorrectParameterError):
+            lpost([2,3])
+
+    def test_logprior(self):
+        t0 = [10.0]
+
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        lp_test = lpost.logprior(t0)
+        lp = np.log(scipy.stats.norm(self.countrate, self.countrate).pdf(t0))
+        assert lp == lp_test
+
+    def test_loglikelihood(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        loglike = np.sum(-mean_model + self.y*np.log(mean_model) - scipy_gammaln(self.y+1))
+
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=False)
+
+        assert np.isclose(loglike, loglike_test)
+
+
+    def test_negative_loglikelihood(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        loglike = -np.sum(-mean_model + self.y*np.log(mean_model) - scipy_gammaln(self.y+1))
+
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=True)
+
+        assert np.isclose(loglike, loglike_test)
+
+    def test_posterior(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=False)
+
+        loglike = np.sum(-mean_model + self.y*np.log(mean_model) - scipy_gammaln(self.y+1))
+        logprior = np.log(scipy.stats.norm(self.countrate, self.countrate).pdf(t0))
+
+        post = loglike + logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_negative_posterior(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = PoissonPosterior(self.x, self.y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=True)
+
+        loglike = np.sum(-mean_model + self.y*np.log(mean_model) - scipy_gammaln(self.y+1))
+        logprior = np.log(scipy.stats.norm(self.countrate, self.countrate).pdf(t0))
+
+        post = -loglike - logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_counts_are_nan(self):
+        y = np.nan * np.ones(self.x.shape[0])
+
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = PoissonPosterior(self.x, y, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert np.isclose(lpost(t0), logmin, 1e-5)
+
+class TestGaussianPosterior(object):
+
+    @classmethod
+    def setup_class(cls):
+
+        nx = 1000000
+        cls.x = np.arange(nx)
+        cls.countrate = 10.0
+        cls.cerr = 2.0
+        cls.y = np.random.normal(cls.countrate, cls.cerr, size=cls.x.shape[0])
+        cls.yerr = np.ones_like(cls.y)*cls.cerr
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.countrate, scale=cls.cerr).pdf(amplitude)
+
+        cls.priors = {"amplitude":p_amplitude}
+
+    def test_logprior_fails_without_prior(self):
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+
+        with pytest.raises(AttributeError):
+            lpost.logprior([10])
+
+    def test_making_posterior(self):
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert lpost.x.all() == self.x.all()
+        assert lpost.y.all() == self.y.all()
+
+    def test_correct_number_of_parameters(self):
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        with pytest.raises(IncorrectParameterError):
+            lpost([2,3])
+
+    def test_logprior(self):
+        t0 = [10.0]
+
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        lp_test = lpost.logprior(t0)
+        lp = np.log(scipy.stats.norm(self.countrate, self.cerr).pdf(t0))
+        assert lp == lp_test
+
+    def test_loglikelihood(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        loglike = np.sum(-0.5*np.log(2.*np.pi) - np.log(self.yerr) - \
+                         0.5*((self.y - mean_model)/self.yerr)**2.0)
+
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=False)
+
+        assert np.isclose(loglike, loglike_test)
+
+
+    def test_negative_loglikelihood(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        loglike = -np.sum(-0.5*np.log(2.*np.pi) - np.log(self.yerr) - \
+                         0.5*((self.y - mean_model)/self.yerr)**2.0)
+
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=True)
+
+        assert np.isclose(loglike, loglike_test)
+
+    def test_posterior(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=False)
+
+        loglike = np.sum(-0.5*np.log(2.*np.pi) - np.log(self.yerr) - \
+                         0.5*((self.y - mean_model)/self.yerr)**2.0)
+        logprior = np.log(scipy.stats.norm(self.countrate, self.cerr).pdf(t0))
+
+        post = loglike + logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_negative_posterior(self):
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.x)
+
+        lpost = GaussianPosterior(self.x, self.y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=True)
+
+        loglike = np.sum(-0.5*np.log(2.*np.pi) - np.log(self.yerr) - \
+                         0.5*((self.y - mean_model)/self.yerr)**2.0)
+        logprior = np.log(scipy.stats.norm(self.countrate, self.cerr).pdf(t0))
+
+        post = -loglike - logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_counts_are_nan(self):
+        y = np.nan * np.ones(self.x.shape[0])
+
+        t0 = [10.0]
+        self.model.amplitude = t0[0]
+
+        lpost = GaussianPosterior(self.x, y, self.yerr, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert np.isclose(lpost(t0), logmin, 1e-5)
+
+
+
+class TestPerPosteriorAveragedPeriodogram(object):
+
+    @classmethod
+    def setup_class(cls):
+
+        cls.m = 10
+        nfreq = 1000000
+        freq = np.arange(nfreq)
+        noise = scipy.stats.chi2(2.*cls.m).rvs(size=nfreq)/np.float(cls.m)
+        power = noise
+
+        ps = Powerspectrum()
+        ps.freq = freq
+        ps.power = power
+        ps.m = cls.m
+        ps.df = freq[1]-freq[0]
+        ps.norm = "leahy"
+
+
+        cls.ps = ps
+        cls.a_mean, cls.a_var = 2.0, 1.0
+
+        cls.model = models.Const1D()
+
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=cls.a_mean, scale=cls.a_var).pdf(amplitude)
+
+        cls.priors = {"amplitude":p_amplitude}
+
+    def test_logprior_fails_without_prior(self):
+        lpost = PSDPosterior(self.ps, self.model)
+
+        with pytest.raises(AttributeError):
+            lpost.logprior([1])
+
+    def test_making_posterior(self):
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert lpost.x.all() == self.ps.freq.all()
+        assert lpost.y.all() == self.ps.power.all()
+
+    def test_correct_number_of_parameters(self):
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        with pytest.raises(IncorrectParameterError):
+            lpost([2,3])
+
+    def test_logprior(self):
+        t0 = [2.0]
+
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        lp_test = lpost.logprior(t0)
+        lp = np.log(scipy.stats.norm(2.0, 1.0).pdf(t0))
+        assert lp == lp_test
+
+    def test_loglikelihood(self):
+        t0 = [2.0]
+        self.model.amplitude = t0[0]
+        mean_model = self.model(self.ps.freq)
+
+        loglike = -2.0*self.m*(np.sum(np.log(mean_model)) +
+                               np.sum(self.ps.power/mean_model) +
+                               np.sum((2.0 / (2. * self.m) - 1.0) *
+                                      np.log(self.ps.power)))
+
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=False)
+
+        assert np.isclose(loglike, loglike_test)
+
+
+    def test_negative_loglikelihood(self):
+        t0 = [2.0]
+        self.model.amplitude = t0[0]
+
+        mean_model = self.model(self.ps.freq)
+
+        loglike = 2.0*self.m*(np.sum(np.log(mean_model)) +
+                               np.sum(self.ps.power/mean_model) +
+                               np.sum((2.0 / (2. * self.m) - 1.0) *
+                                      np.log(self.ps.power)))
+
+
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        loglike_test = lpost.loglikelihood(t0, neg=True)
+
+        assert np.isclose(loglike, loglike_test)
+
+    def test_posterior(self):
+        t0 = [2.0]
+        self.model.amplitude = t0[0]
+
+        mean_model = self.model(self.ps.freq)
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=False)
+
+        loglike = -2.0*self.m*(np.sum(np.log(mean_model)) +
+                               np.sum(self.ps.power/mean_model) +
+                               np.sum((2.0 / (2. * self.m) - 1.0) *
+                                      np.log(self.ps.power)))
+
+        logprior = np.log(scipy.stats.norm(2.0, 1.0).pdf(t0))
+        post = loglike + logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_negative_posterior(self):
+        t0 = [2.0]
+        self.model.amplitude = t0[0]
+
+        mean_model = self.model(self.ps.freq)
+        lpost = PSDPosterior(self.ps, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        post_test = lpost(t0, neg=True)
+
+        loglike = -2.0*self.m*(np.sum(np.log(mean_model)) +
+                               np.sum(self.ps.power/mean_model) +
+                               np.sum((2.0 / (2. * self.m) - 1.0) *
+                                      np.log(self.ps.power)))
+
+        logprior = np.log(scipy.stats.norm(2.0, 1.0).pdf(t0))
+        post = -loglike - logprior
+
+        assert np.isclose(post_test, post, atol=1.e-10)
+
+    def test_counts_are_nan(self):
+        y = np.nan * np.ones_like(self.ps.freq)
+
+        ps_nan = copy.copy(self.ps)
+        ps_nan.power = np.nan*np.ones_like(self.ps.freq)
+
+        t0 = [2.0]
+        m = self.model(self.ps.freq[1:], t0)
+        lpost = PSDPosterior(ps_nan, self.model)
+        lpost.logprior = set_logprior(lpost, self.priors)
+
+        assert np.isclose(lpost(t0), logmin, 1e-5)
diff --git a/stingray/modeling/tests/test_scripts.py b/stingray/modeling/tests/test_scripts.py
new file mode 100644
index 00000000..66e389d6
--- /dev/null
+++ b/stingray/modeling/tests/test_scripts.py
@@ -0,0 +1,79 @@
+import numpy as np
+import scipy.stats
+import copy
+
+from astropy.tests.helper import pytest
+from astropy.modeling import models
+from scipy.special import gammaln as scipy_gammaln
+
+from stingray import Lightcurve, Powerspectrum
+from stingray.modeling import Posterior, PSDPosterior, PoissonPosterior, GaussianPosterior
+from stingray.modeling import set_logprior
+from stingray.modeling.posterior import logmin
+from stingray.modeling.posterior import IncorrectParameterError
+
+from stingray.modeling.scripts import fit_lorentzians
+
+
+class TestFitLorentzians(object):
+
+    @classmethod
+    def setup_class(cls):
+
+        np.random.seed(150)
+        cls.nlor = 3
+
+        cls.x_0_0 = 0.5
+        cls.x_0_1 = 2.0
+        cls.x_0_2 = 7.5
+
+        cls.amplitude_0 = 200.0
+        cls.amplitude_1 = 100.0
+        cls.amplitude_2 = 50.0
+
+        cls.fwhm_0 = 0.1
+        cls.fwhm_1 = 1.0
+        cls.fwhm_2 = 0.5
+
+        cls.whitenoise = 2.0
+
+        cls.model = models.Lorentz1D(cls.amplitude_0, cls.x_0_0, cls.fwhm_0) + \
+                    models.Lorentz1D(cls.amplitude_1, cls.x_0_1, cls.fwhm_1) + \
+                    models.Lorentz1D(cls.amplitude_2, cls.x_0_2, cls.fwhm_2) + \
+                    models.Const1D(cls.whitenoise)
+
+        freq = np.linspace(0.01, 10.0, 10.0 / 0.01)
+        p = cls.model(freq)
+        noise = np.random.exponential(size=len(freq))
+
+        power = p * noise
+        cls.ps = Powerspectrum()
+        cls.ps.freq = freq
+        cls.ps.power = power
+        cls.ps.df = cls.ps.freq[1] - cls.ps.freq[0]
+        cls.ps.m = 1
+
+        cls.t0 = [200.0, 0.5, 0.1, 100.0, 2.0, 1.0, 50.0, 7.5, 0.5, 2.0]
+
+        cls.parest, cls.res = fit_lorentzians(cls.ps, cls.nlor, cls.t0)
+
+    def test_function_creates_right_number_of_lorentians(self):
+        assert (len(self.parest.lpost.model.parameters)-1)/3 == self.nlor
+
+    def test_correct_parameters_without_whitenoise(self):
+        parest, res = fit_lorentzians(self.ps, self.nlor, self.t0[:-1],
+                                      fit_whitenoise=False)
+
+        assert len(parest.lpost.model.parameters)/3 == self.nlor
+
+    def test_parameters_in_right_ballpark(self):
+
+        true_pars = [self.amplitude_0, self.x_0_0, self.fwhm_0,
+                     self.amplitude_1, self.x_0_1, self.fwhm_1,
+                     self.amplitude_2, self.x_0_2, self.fwhm_2,
+                     self.whitenoise]
+
+        for t, p in zip(true_pars, self.res.p_opt):
+            print(str(t) + "\t" + str(p))
+
+        assert np.all(np.isclose(true_pars, self.res.p_opt, rtol=0.5))
\ No newline at end of file
diff --git a/stingray/pulse/tests/test_pulse.py b/stingray/pulse/tests/test_pulse.py
index 99308763..34a29e13 100644
--- a/stingray/pulse/tests/test_pulse.py
+++ b/stingray/pulse/tests/test_pulse.py
@@ -81,7 +81,7 @@ class TestAll(object):
         nbin = 16
         expo = phase_exposure(start_time, stop_time, period, nbin)
         expected = np.ones(nbin)
-        expected[nbin/2:] = 0
+        expected[nbin//2:] = 0
         np.testing.assert_array_almost_equal(expo, expected)
 
     def test_phase_exposure3(self):
@@ -92,7 +92,7 @@ class TestAll(object):
         nbin = 16
         expo = phase_exposure(start_time, stop_time, period, nbin, gtis=gtis)
         expected = np.ones(nbin)
-        expected[nbin/2:] = 0
+        expected[nbin//2:] = 0
         np.testing.assert_array_almost_equal(expo, expected)
 
     def test_phase_exposure4(self):