Repository: StingraySoftware/stingray
Pull request #228 (2017-07-26 21:00:29)

Contains all changes including additions of convenience functions to …
…search for periods and compare models, and associated tests.

Only changes this time, no rebase. This is a continuation from PR #179. 


Patch:
diff --git a/.travis.yml b/.travis.yml
index 0209d94c..9c45139e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -30,7 +30,7 @@ env:
         # to repeat them for all configurations.
         - ASTROPY_VERSION=stable
         - SETUP_CMD='test'
-        - PIP_DEPENDENCIES='emcee statsmodels corner'
+        - PIP_DEPENDENCIES='emcee statsmodels corner pytest-catchlog'
         # 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=''
@@ -45,10 +45,10 @@ matrix:
     fast_finish: true
     include:
         - python: 3.5
-          env: PIP_DEPENDENCIES='git+https://git@github.com/nanograv/pint.git jplephem emcee statsmodels corner' SETUP_CMD='test --remote-data --coverage'
+          env: PIP_DEPENDENCIES='git+https://git@github.com/nanograv/pint.git jplephem emcee statsmodels corner pytest-catchlog' SETUP_CMD='test --remote-data --coverage'
 
         - python: 2.7
-          env: PIP_DEPENDENCIES='git+https://git@github.com/nanograv/pint.git jplephem emcee statsmodels corner' SETUP_CMD='test --remote-data'
+          env: PIP_DEPENDENCIES='git+https://git@github.com/nanograv/pint.git jplephem emcee statsmodels corner pytest-catchlog' SETUP_CMD='test --remote-data'
 
         #Try without importing h5py but importing numba
         - python: 2.7
@@ -59,7 +59,15 @@ matrix:
 
         # Try without using corner
         - python: 3.5
-          env: PIP_DEPENDENCIES='emcee statsmodels' SETUP_CMD='test --coverage'
+          env: PIP_DEPENDENCIES='emcee statsmodels pytest-catchlog' SETUP_CMD='test --coverage'
+
+        # Try without using emcee
+        - python: 3.5
+          env: PIP_DEPENDENCIES='statsmodels pytest-catchlog' SETUP_CMD='test --coverage'
+
+        # Try without using statsmodels
+        - python: 3.5
+          env: PIP_DEPENDENCIES='emcee pytest-catchlog' SETUP_CMD='test --coverage'
 
         - python: 3.4
           env: NUMPY_VERSION=1.11
diff --git a/appveyor.yml b/appveyor.yml
index 352b9aef..e6f20dfe 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -11,7 +11,7 @@ environment:
                           # 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: "emcee statsmodels"
+        PIP_DEPENDENCIES: "emcee statsmodels pytest-catchlog"
 
     matrix:
 
diff --git a/stingray/modeling/parameterestimation.py b/stingray/modeling/parameterestimation.py
index f48bc4fc..9f01d807 100644
--- a/stingray/modeling/parameterestimation.py
+++ b/stingray/modeling/parameterestimation.py
@@ -26,6 +26,8 @@ try:
 except ImportError:
     use_corner = False
 
+import logging
+
 import numpy as np
 import scipy
 import scipy.optimize
@@ -39,10 +41,12 @@ try:
 except ImportError:
     comp_hessian = False
 
-from stingray.modeling.posterior import Posterior, LogLikelihood
 from astropy.modeling.fitting import _fitter_to_model_params, \
     _model_to_fit_params, _validate_model, _convert_input
 
+from stingray.modeling.posterior import Posterior, PSDPosterior, \
+    LogLikelihood, PSDLogLikelihood
+
 
 class OptimizationResults(object):
 
@@ -84,7 +88,7 @@ class OptimizationResults(object):
         else:
             if comp_hessian:
                 # calculate Hessian approximating with finite differences
-                print("Approximating Hessian with finite differences ...")
+                logging.info("Approximating Hessian with finite differences ...")
 
                 phess = approx_hess(self.p_opt, lpost)
 
@@ -131,7 +135,7 @@ class OptimizationResults(object):
 
     def print_summary(self, lpost):
 
-        print("The best-fit model parameters plus errors are:")
+        logging.info("The best-fit model parameters plus errors are:")
 
         fixed = [lpost.model.fixed[n] for n in lpost.model.param_names]
         tied = [lpost.model.tied[n] for n in lpost.model.param_names]
@@ -139,38 +143,39 @@ class OptimizationResults(object):
 
         parnames = [n for n, f in zip(lpost.model.param_names,
                                       np.logical_or(fixed, tied)) \
-                    if f is False]
+                    if not f]
 
         all_parnames = [n for n in lpost.model.param_names]
         for i, par in enumerate(all_parnames):
-            print("{:3}) Parameter {:<20}: ".format(i, par), end="")
+            logging.info("{:3}) Parameter {:<20}: ".format(i, par))
 
             if par in parnames:
                 idx = parnames.index(par)
-                print("{:<20.5f} +/- {:<20.5f} ".format(self.p_opt[idx],
-                                                  self.err[idx]), end="")
-                print("[{:>10} {:>10}]".format(str(bounds[i][0]),
+
+                logging.info("{:<20.5f} +/- {:<20.5f} ".format(self.p_opt[idx],
+                                                  self.err[idx]))
+                logging.info("[{:>10} {:>10}]".format(str(bounds[i][0]),
                                                str(bounds[i][1])))
             elif fixed[i]:
-                print("{:<20.5f} (Fixed) ".format(lpost.model.parameters[i]))
+                logging.info("{:<20.5f} (Fixed) ".format(lpost.model.parameters[i]))
             elif tied[i]:
-                print("{:<20.5f} (Tied) ".format(lpost.model.parameters[i]))
+                logging.info("{:<20.5f} (Tied) ".format(lpost.model.parameters[i]))
 
-        print("\n")
+        logging.info("\n")
 
-        print("Fitting statistics: ")
-        print(" -- number of data points: %i"%(len(lpost.x)))
+        logging.info("Fitting statistics: ")
+        logging.info(" -- 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: " +
+        logging.info(" -- Deviance [-2 log L] D = %f.3"%self.deviance)
+        logging.info(" -- The Akaike Information Criterion of the model is: " +
               str(self.aic) + ".")
 
-        print(" -- The Bayesian Information Criterion of the model is: " +
+        logging.info(" -- The Bayesian Information Criterion of the model is: " +
               str(self.bic) + ".")
 
         try:
@@ -178,14 +183,14 @@ class OptimizationResults(object):
         except AttributeError:
             self._compute_statistics(lpost)
 
-        print(" -- The figure-of-merit function for this model " +
+        logging.info(" -- 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)
+        logging.info(" -- Summed Residuals S = %f.5f"%self.sobs)
+        logging.info(" -- Expected S ~ %f.5 +/- %f.5"%(self.sexp, self.ssd))
+        logging.info(" -- merit function (SSE) M = %f.5f \n\n"%self.merit)
 
         return
 
@@ -252,8 +257,9 @@ class ParameterEstimation(object):
                             "Posterior or LogLikelihoood.")
 
         newmod = lpost.model.copy()
-        newmod.parameters = t0
-        p0, _ = _model_to_fit_params(newmod)
+
+        p0=t0
+
         # p0 will be shorter than t0, if there are any frozen/tied parameters
         # this has to match with the npar attribute.
         if not len(p0) == lpost.npar:
@@ -281,38 +287,76 @@ class ParameterEstimation(object):
             # perturb parameters slightly
             t0_p = np.random.multivariate_normal(p0, np.diag(np.abs(p0)/100.))
 
-            # print(lpost.model, dir(lpost.model), lpost.model.parameter_constraints, lpost.model.param_names)
             params = [getattr(newmod,name) for name in newmod.param_names]
             bounds = [p.bounds for p in params if not np.any([p.tied, p.fixed])]
-            # print(params, bounds)
+
+            if len(bounds) > 0 and self.fitmethod not in ["L-BFGS-B",
+                                                          "TNC",
+                                                          "SLSQP"]:
+                logging.warning("Fitting method %s "%self.fitmethod +
+                                "cannot incorporate the bounds you set!")
+
+            if len(bounds) == 0 or self.fitmethod not in ["L-BFGS-B",
+                                                          "TNC",
+                                                          "SLSQP"]:
+                use_bounds = False
+            else:
+                use_bounds = True
+
+
             # 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,
-                                              bounds=bounds,
-                                              **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,
+                if use_bounds:
+                    opt = scipy.optimize.minimize(lpost, t0_p,
                                                   method=self.fitmethod,
                                                   args=args, tol=1.e-10,
                                                   bounds=bounds,
                                                   **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,
+                else:
+                    opt = scipy.optimize.minimize(lpost, t0_p,
                                                   method=self.fitmethod,
                                                   args=args, tol=1.e-10,
-                                                  bounds=bounds,
                                                   **scipy_optimize_options)
 
+
+            # if max_post is False, then do a Maximum Likelihood Fit
+            else:
+                if isinstance(lpost, Posterior):
+                    if use_bounds:
+                        # This could be a `Posterior` object
+                        opt = scipy.optimize.minimize(lpost.loglikelihood, t0_p,
+                                                      method=self.fitmethod,
+                                                      args=args, tol=1.e-10,
+                                                      bounds=bounds,
+                                                      **scipy_optimize_options)
+                    else:
+                        opt = scipy.optimize.minimize(lpost.loglikelihood, t0_p,
+                                                      method=self.fitmethod,
+                                                      args=args, tol=1.e-10,
+                                                      **scipy_optimize_options)
+
+
+                elif isinstance(lpost, LogLikelihood):
+
+                    if use_bounds:
+                        # 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,
+                                                      #bounds=bounds,
+                                                      **scipy_optimize_options)
+
+
+                    else:
+                        opt = scipy.optimize.minimize(lpost.evaluate, t0_p,
+                                                      method=self.fitmethod,
+                                                      args=args, tol=1.e-10,
+                                                      **scipy_optimize_options)
+
             funcval = opt.fun
             i += 1
 
@@ -347,6 +391,17 @@ class ParameterEstimation(object):
             If True, set the internal state to do the optimization with the
             log-likelihood rather than the log-posterior.
 
+        Returns
+        -------
+        lrt : float
+            The likelihood ratio for model 2 and model 1
+
+        res1 : OptimizationResults object
+            Contains the result of fitting `lpost1`
+
+        res2 : OptimizationResults object
+            Contains the results of fitting `lpost2`
+
         """
 
         self.max_post = max_post
@@ -358,7 +413,7 @@ class ParameterEstimation(object):
         # compute log likelihood ratio as difference between the deviances
         lrt = res1.deviance - res2.deviance
 
-        return lrt
+        return lrt, res1, res2
 
     def sample(self, lpost, t0, cov=None,
                nwalkers=500, niter=100, burnin=100, threads=1,
@@ -370,11 +425,15 @@ class ParameterEstimation(object):
 
         Parameters
         ----------
+        lpost : instance of a Posterior subclass
+            and instance of class Posterior or one of its subclasses
+            that defines the function to be minized (either in loglikelihood
+            or logposterior)
+
         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
@@ -446,10 +505,213 @@ class ParameterEstimation(object):
             res.print_results()
 
         if plot:
-            res.plot_results(namestr + "_corner.pdf")
+            fig = res.plot_results(fig=None, save_plot=True,
+                                   filename=namestr + "_corner.pdf")
 
         return res
 
+    def _generate_model(self, lpost, pars):
+        """
+        Helper function that generates a fake PSD similar to the
+        one in the data, but with different parameters.
+
+        Parameters
+        ----------
+        lpost : instance of a Posterior or LogLikelihood subclass
+            The object containing the relevant information about the
+            data and the model
+
+        pars : iterable
+            A list of parameters to be passed to lpost.model in oder
+            to generate a model data set.
+
+        Returns:
+        --------
+        model_data : numpy.ndarray
+            An array of model values for each bin in lpost.x
+
+        """
+
+        assert isinstance(lpost, LogLikelihood) or isinstance(lpost, Posterior), \
+            "lpost must be of type LogLikelihood or Posterior or one of its " \
+            "subclasses!"
+
+        # assert pars is of correct length
+        assert len(pars) == lpost.npar, "pars must be a list " \
+                                        "of %i parameters"%lpost.npar
+        # get the model
+        m = lpost.model
+
+        # reset the parameters
+        _fitter_to_model_params(m, pars)
+
+        # make a model spectrum
+        model_data = lpost.model(lpost.x)
+
+        return model_data
+
+    @staticmethod
+    def _compute_pvalue(obs_val, sim):
+        """
+        Compute the p-value given an observed value of a test statistic
+        and some simulations of that same test statistic.
+
+        Parameters
+        ----------
+        obs_value : float
+            The observed value of the test statistic in question
+
+        sim: iterable
+            A list or array of simulated values for the test statistic
+
+        Returns
+        -------
+        pval : float [0, 1]
+            The p-value for the test statistic given the simulations.
+
+        """
+
+        # cast the simulations as a numpy array
+        sim = np.array(sim)
+
+        # find all simulations that are larger than
+        # the observed value
+        ntail = sim[sim > obs_val].shape[0]
+
+        # divide by the total number of simulations
+        pval = np.float(ntail) / np.float(sim.shape[0])
+
+        return pval
+
+
+    def simulate_lrts(self, s_all, lpost1, t1, lpost2, t2, max_post=True,
+                      seed=None):
+        """
+        Simulate likelihood ratios.
+        For details, see definitions in the subclasses that implement this
+        task.
+        """
+        raise NotImplementedError("The behaviour of `simulate_lrts` should be defined "
+                        "in the subclass appropriate for your problem, not in "
+                        "this super class!")
+
+    def calibrate_lrt(self, lpost1, t1, lpost2, t2, sample=None, neg=True,
+                      max_post=False,
+                      nsim=1000, niter=200, nwalkers=500, burnin=200,
+                      namestr="test", seed=None):
+
+        """
+        Calibrate the outcome of a Likelihood Ratio Test via MCMC.
+
+        In order to compare models via likelihood ratio test, one generally
+        aims to compute a p-value for the null hypothesis (generally the
+        simpler model). There are two special cases where the theoretical
+        distribution used to compute that p-value analytically given the
+        observed likelihood ratio (a chi-square distribution) is not
+        applicable:
+        * the models are not nested (i.e. Model 1 is not a special, simpler
+        case of Model 2),
+        * the parameter values fixed in Model 2 to retrieve Model 1 are at the
+        edges of parameter space (e.g. if one must set, say, an amplitude to
+        zero in order to remove a component in the more complex model, and
+        negative amplitudes are excluded a priori)
+
+        In these cases, the observed likelihood ratio must be calibrated via
+         simulations of the simpler model (Model 1), using MCMC to take into
+        account the uncertainty in the parameters. This function does
+        exactly that: it computes the likelihood ratio for the observed data,
+        and produces simulations to calibrate the likelihood ratio and
+        compute a p-value for observing the data under the assumption that
+        Model 1 istrue.
+
+        If `max_post=True`, the code will use MCMC to sample the posterior
+        of the parameters and simulate fake data from there.
+
+        If `max_post=False`, the code will use the covariance matrix derived
+        from the fit to simulate data sets for comparison.
+
+        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.
+
+
+        Returns
+        -------
+        pvalue : float [0,1]
+            p-value 'n stuff
+
+        """
+
+        # compute the observed likelihood ratio
+        lrt_obs, res1, res2 = self.compute_lrt(lpost1, t1,
+                                               lpost2, t2,
+                                               neg=neg,
+                                               max_post=max_post)
+
+        if sample is None:
+            # simulate parameter sets from the simpler model
+            if not max_post:
+                # using Maximum Likelihood, so I'm going to simulate parameters
+                # from a multivariate Gaussian
+
+                # set up the distribution
+                mvn = scipy.stats.multivariate_normal(mean=res1.p_opt,
+                                                      cov=res1.cov)
+
+
+                # sample parameters
+                s_all = mvn.rvs(size=nsim)
+                if lpost1.npar == 1:
+                    s_all = np.atleast_2d(s_all).T
+
+            else:
+                # sample the posterior using MCMC
+                s_mcmc = self.sample(lpost1, res1.p_opt, cov=res1.cov,
+                                       nwalkers=nwalkers, niter=niter,
+                                       burnin=burnin, namestr=namestr)
+
+
+                # pick nsim samples out of the posterior sample
+                s_all = s_mcmc.samples[
+                    np.random.choice(s_mcmc.samples.shape[0], nsim,
+                                     replace=False)]
+
+                #if lpost1.npar == 1:
+                #    s_all = np.atleast_2d(s_all).T
+
+
+        else:
+            s_all = sample
+
+
+        # simulate LRTs
+        # this method is defined in the subclasses!
+        lrt_sim = self.simulate_lrts(s_all, lpost1, t1, lpost2, t2,
+                                      max_post=max_post, seed=seed)
+
+        # now I can compute the p-value:
+        pval = ParameterEstimation._compute_pvalue(lrt_obs, lrt_sim)
+
+        return pval
+
 
 class SamplingResults(object):
 
@@ -496,7 +758,7 @@ class SamplingResults(object):
         try:
             self.acor = sampler.acor
         except emcee.autocorr.AutocorrError:
-            print("Chains too short to compute autocorrelation lengths.")
+            logging.info("Chains too short to compute autocorrelation lengths.")
 
         self.rhat = self._compute_rhat(sampler)
 
@@ -535,22 +797,24 @@ class SamplingResults(object):
 
         """
 
-        print("-- The acceptance fraction is: %f.5"%self.acceptance)
+        logging.info("-- The acceptance fraction is: %f.5"%self.acceptance)
         try:
-            print("-- The autocorrelation time is: %f.5"%self.acor)
+            logging.info("-- The autocorrelation time is: %f.5"%self.acor)
         except AttributeError:
             pass
-        print("R_hat for the parameters is: " + str(self.rhat))
+        logging.info("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")
+        logging.info("-- Posterior Summary of Parameters: \n")
+        logging.info("parameter \t mean \t\t sd \t\t 5% \t\t 95% \n")
+        logging.info("---------------------------------------------\n")
         for i in range(self.ndim):
-            print("theta[" + str(i) + "] \t " +
+            logging.info("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):
+    def plot_results(self, nsamples=1000, fig=None, save_plot=False,
+                     filename="test.pdf"):
+
         """
         Plot some results in a triangle plot.
         If installed, will use `corner` for the plotting
@@ -558,23 +822,37 @@ class SamplingResults(object):
         through pip), if not, uses its own code to make a triangle
         plot.
 
+        By default, this method returns a matplotlib.Figure object, but
+        if `save_plot=True`, the plot can be saved to file automatically,
+
         Parameters
         ----------
 
+        nsamples: int, default 1000
+            The maximum number of samples used for plotting.
+
+        fig: matplotlib.Figure instance, default None
+            If created externally, you can pass a Figure instance to this method.
+            If none is passed, the method will create one internally.
+
+        save_plot: bool, default False
+            If True, save the plot to file with a file name specified by the
+            keyword `filename`. If False, just return the `Figure` object
+
         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,
+            corner.corner(self.samples, labels=None, fig=fig, bins=int(20),
                           quantiles=[0.16, 0.5, 0.84],
                           show_titles=True, title_args={"fontsize": 12})
 
         else:
-            fig = plt.figure(figsize=(15, 15))
+            if fig is None:
+                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)
@@ -617,11 +895,12 @@ class SamplingResults(object):
 
                             ax.contour(xx, yy, zz, 7)
                         except ValueError:
-                            print("Not making contours.")
+                            logging.info("Not making contours.")
 
-        plt.savefig(filename, format='pdf')
-        plt.close()
-        return
+        if save_plot:
+            plt.savefig(filename, format='pdf')
+
+        return fig
 
 
 class PSDParEst(ParameterEstimation):
@@ -661,6 +940,214 @@ class PSDParEst(ParameterEstimation):
 
         return res
 
+    def _generate_data(self, lpost, pars, rng=None):
+        """
+        Generate a fake power spectrum from a model.
+
+        Parameters:
+        ----------
+        lpost : instance of a Posterior or LogLikelihood subclass
+            The object containing the relevant information about the
+            data and the model
+
+        pars : iterable
+            A list of parameters to be passed to lpost.model in oder
+            to generate a model data set.
+
+        Returns:
+        --------
+        sim_ps : stingray.Powerspectrum object
+            The simulated Powerspectrum object
+
+        """
+        # create own random state object
+        if rng is None:
+            rng = np.random.RandomState(None)
+
+        model_spectrum = self._generate_model(lpost, pars)
+
+        # use chi-square distribution to get fake data
+        model_powers = model_spectrum * \
+                       rng.chisquare(2 * self.ps.m,
+                                           size=model_spectrum.shape[0]) \
+                                                / (2. * self.ps.m)
+
+        sim_ps = copy.copy(self.ps)
+
+        sim_ps.power = model_powers
+
+        return sim_ps
+
+    def simulate_lrts(self, s_all, lpost1, t1, lpost2, t2, max_post=True,
+                      seed=None):
+        """
+        Simulate likelihood ratios for two given models based on MCMC samples
+        for the simpler model (i.e. the null hypothesis).
+
+        Parameters
+        ----------
+        s_all : numpy.ndarray of shape (nsamples, lpost1.npar)
+            An array with MCMC samples derived from the null hypothesis model in
+            `lpost1`. Its second dimension must match the number of free
+             parameters in `lpost1.model`.
+
+        lpost1 : LogLikelihood or Posterior subclass object
+            Object containing the null hypothesis model
+
+        t1 : iterable of length `lpost1.npar`
+            A starting guess for fitting the model in `lpost1`
+
+        lpost2 : LogLikelihood or Posterior subclass object
+            Object containing the alternative hypothesis model
+
+        t2 : iterable of length `lpost2.npar`
+            A starting guess for fitting the model in `lpost2`
+
+        max_post : bool, optional, default True
+            If True, then `lpost1` and `lpost2` should be Posterior subclass
+            objects; if False, then `lpost1` and `lpost2` should be
+            LogLikelihood subclass objects
+
+        seed : int, optional default None
+            A seed to initialize the numpy.random.RandomState object to be
+            passed on to `_generate_data`. Useful for producing exactly
+            reproducible results
+
+        Returns
+        -------
+        lrt_sim : numpy.ndarray
+            An array with the simulated likelihood ratios for the simulated
+            data
+
+        """
+
+        assert lpost1.__class__ == lpost2.__class__, "Both LogLikelihood or " \
+                                                     "Posterior objects must be " \
+                                                     "of the same class!"
+
+        nsim = s_all.shape[0]
+        lrt_sim = np.zeros(nsim)
+
+        rng = np.random.RandomState(seed)
+
+        # now I can loop over all simulated parameter sets to generate a PSD
+        for i, s in enumerate(s_all):
+
+            # generate fake PSD
+            sim_ps = self._generate_data(lpost1, s, rng)
+
+            # make LogLikelihood objects for both:
+            if isinstance(lpost1, LogLikelihood):
+                sim_lpost1 = PSDLogLikelihood(sim_ps.freq, sim_ps.power,
+                                              model=lpost1.model)
+                sim_lpost2 = PSDLogLikelihood(sim_ps.freq, sim_ps.power,
+                                              model=lpost2.model, m=sim_ps.m)
+                neg = True
+            else:
+                # make a Posterior object
+                sim_lpost1 = PSDPosterior(sim_ps.freq, sim_ps.power,
+                                          lpost1.model, m=sim_ps.m)
+                sim_lpost1.logprior = lpost1.logprior
+
+                sim_lpost2 = PSDPosterior(sim_ps.freq, sim_ps.power,
+                                          lpost2.model, m=sim_ps.m)
+
+                sim_lpost2.logprior = lpost2.logprior
+                neg=False
+
+            parest_sim = PSDParEst(sim_ps, max_post=max_post)
+
+            lrt_sim[i], _, _ = parest_sim.compute_lrt(sim_lpost1, t1,
+                                                      sim_lpost2, t2,
+                                                      neg=neg,
+                                                      max_post=max_post)
+        return lrt_sim
+
+
+    def calibrate_highest_outlier(self, lpost, t0, sample=None,
+                                  max_post=False,
+                                  nsim=1000, niter=200, nwalkers=500,
+                                  burnin=200, namestr="test", seed=None):
+
+        """
+
+        """
+        # fit the model to the data
+        res = self.fit(lpost, t0, neg=True)
+
+        # find the highest data/model outlier:
+        out_high, _, _ = self._compute_highest_outlier(lpost, res)
+        # simulate parameter sets from the simpler model
+        if not max_post:
+            # using Maximum Likelihood, so I'm going to simulate parameters
+            # from a multivariate Gaussian
+
+            # set up the distribution
+            mvn = scipy.stats.multivariate_normal(mean=res.p_opt,
+                                                  cov=res.cov)
+
+            if lpost.npar == 1:
+                # sample parameters
+                s_all = np.atleast_2d(mvn.rvs(size=nsim)).T
+
+            else:
+                s_all = mvn.rvs(size=nsim)
+
+        else:
+            if sample is None:
+                # sample the posterior using MCMC
+                sample = self.sample(lpost, res.p_opt, cov=res.cov,
+                                       nwalkers=nwalkers, niter=niter,
+                                       burnin=burnin, namestr=namestr)
+
+            # pick nsim samples out of the posterior sample
+            s_all = sample.samples[
+                np.random.choice(sample.samples.shape[0], nsim, replace=False)]
+
+        # simulate LRTs
+        # this method is defined in the subclasses!
+        out_high_sim = self.simulate_highest_outlier(s_all, lpost, t0,
+                                                max_post=max_post, seed=seed)
+        # now I can compute the p-value:
+        pval = ParameterEstimation._compute_pvalue(out_high, out_high_sim)
+
+        return pval
+
+    def simulate_highest_outlier(self, s_all, lpost, t0, max_post=True,
+                                 seed=None):
+
+        # the number of simulations
+        nsim = s_all.shape[0]
+
+        # empty array for the simulation results
+        max_y_all = np.zeros(nsim)
+
+        rng = np.random.RandomState(seed)
+
+        # now I can loop over all simulated parameter sets to generate a PSD
+        for i, s in enumerate(s_all):
+
+            # generate fake PSD
+            sim_ps = self._generate_data(lpost, s, rng=rng)
+
+            # make LogLikelihood objects for both:
+            if not max_post:
+                sim_lpost = PSDLogLikelihood(sim_ps.freq, sim_ps.power,
+                                              model=lpost.model, m=sim_ps.m)
+            else:
+                # make a Posterior object
+                sim_lpost = PSDPosterior(sim_ps.freq, sim_ps.power,
+                                         lpost.model, m=sim_ps.m)
+                sim_lpost.logprior = lpost.logprior
+
+            parest_sim = PSDParEst(sim_ps, max_post=max_post)
+
+            res = parest_sim.fit(sim_lpost, t0, neg=True)
+            max_y_all[i], maxfreq, maxind = self._compute_highest_outlier(sim_lpost,
+                                                               res,
+                                                               nmax=1)
+        return np.hstack(max_y_all)
+
     def _compute_highest_outlier(self, lpost, res, nmax=1):
 
         residuals = 2.0 * lpost.y/ res.mfit
@@ -679,7 +1166,7 @@ class PSDParEst(ParameterEstimation):
 
     @staticmethod
     def _find_outlier(xdata, ratio, max_y):
-        max_ind = np.where(ratio == max_y)[0][0]+1
+        max_ind = np.where(ratio == max_y)[0][0]
         max_x = xdata[max_ind]
 
         return max_x, max_ind
@@ -688,7 +1175,7 @@ class PSDParEst(ParameterEstimation):
                  namestr='test', log=False):
 
         if not can_plot:
-            print("No matplotlib imported. Can't plot!")
+            logging.info("No matplotlib imported. Can't plot!")
         else:
             # make a figure
             f = plt.figure(figsize=(12, 10))
@@ -767,7 +1254,7 @@ class PSDParEst(ParameterEstimation):
             else:
                 s2.plot(self.ps.freq, pldif, color='black',
                         linestyle='steps-mid')
-                s2.plot(self.ps.power, np.ones(self.ps.power.shape[0]),
+                s2.plot(self.ps.freq, np.ones_like(self.ps.power),
                         color='blue', lw=2)
 
                 s2.set_xscale("log")
diff --git a/stingray/modeling/posterior.py b/stingray/modeling/posterior.py
index 81553bf7..bc979fd0 100644
--- a/stingray/modeling/posterior.py
+++ b/stingray/modeling/posterior.py
@@ -1,9 +1,13 @@
 from __future__ import division
 
 import abc
+import warnings
 
 import numpy as np
 import six
+
+np.seterr('warn')
+
 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
@@ -11,6 +15,7 @@ from astropy.modeling import models
 
 from stingray import Lightcurve, Powerspectrum
 
+
 # TODO: Add checks and balances to code
 
 #from stingray.modeling.parametricmodels import logmin
@@ -75,7 +80,7 @@ def set_logprior(lpost, priors):
 
     Instantiate the posterior:
 
-    >>> lpost = PSDPosterior(ps, pl)
+    >>> lpost = PSDPosterior(ps.freq, ps.power, pl, m=ps.m)
 
     Define the priors:
 
@@ -126,7 +131,12 @@ def set_logprior(lpost, priors):
         # correctly!
         for pname in lpost.model.param_names:
             if not lpost.model.fixed[pname]:
-                logp += np.log(priors[pname](t0[ii]))
+                with warnings.catch_warnings(record=True) as out:
+                    logp += np.log(priors[pname](t0[ii]))
+                    if len(out) > 0:
+                        if isinstance(out[0].message, RuntimeWarning):
+                            logp = np.nan
+
                 ii += 1
 
         if not np.isfinite(logp):
@@ -143,7 +153,7 @@ def set_logprior(lpost, priors):
 @six.add_metaclass(abc.ABCMeta)
 class LogLikelihood(object):
 
-    def __init__(self, x, y, model):
+    def __init__(self, x, y, model, **kwargs):
         """
         x : iterable
             x-coordinate of the data. Could be multi-dimensional.
@@ -151,8 +161,13 @@ class LogLikelihood(object):
         y : iterable
             y-coordinate of the data. Could be multi-dimensional.
 
-        model: probably astropy.modeling.FittableModel instance
+        model : probably astropy.modeling.FittableModel instance
             Your model
+
+        kwargs :
+            keyword arguments specific to the individual sub-classes. For
+            details, see the respective docstrings for each subclass
+
         """
         self.x = x
         self.y = y
@@ -276,17 +291,18 @@ class PoissonLogLikelihood(LogLikelihood):
 
 class PSDLogLikelihood(LogLikelihood):
 
-    def __init__(self, x, y, model, m=1):
+    def __init__(self, freq, power, model, m=1):
         """
         A Gaussian likelihood.
 
         Parameters
         ----------
-        x : iterable
-            x-coordinate of the data
+        freq: iterable
+            Array with frequencies
 
-        y : iterable
-            y-coordinte of the data
+        power: iterable
+            Array with (averaged/singular) powers corresponding to the
+            frequencies in `freq`
 
         model: an Astropy Model instance
             The model to use in the likelihood.
@@ -296,9 +312,8 @@ class PSDLogLikelihood(LogLikelihood):
 
         """
 
-        self.x = x
-        self.y = y
-        self.model = model
+        LogLikelihood.__init__(self, freq, power, model)
+
         self.m = m
         self.npar = 0
         for pname in self.model.param_names:
@@ -315,16 +330,18 @@ class PSDLogLikelihood(LogLikelihood):
 
         mean_model = self.model(self.x)
 
-        if self.m == 1:
-            loglike = -np.sum(np.log(mean_model)) - \
-                      np.sum(self.y/mean_model)
+        with warnings.catch_warnings(record=True) as out:
 
-        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 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
@@ -336,7 +353,7 @@ class PSDLogLikelihood(LogLikelihood):
 
 class Posterior(object):
 
-    def __init__(self, x, y, model):
+    def __init__(self, x, y, model, **kwargs):
         """
         Define a posterior object.
 
@@ -369,6 +386,9 @@ class Posterior(object):
             The parametric model supposed to represent the data. For details
             see the astropy.modeling documentation
 
+        kwargs :
+            keyword arguments related to the subclases of `Posterior`. For
+            details, see the documentation of the individual subclasses
 
         References
         ----------
@@ -416,7 +436,7 @@ class Posterior(object):
 
 class PSDPosterior(Posterior):
 
-    def __init__(self, ps, model, priors=None):
+    def __init__(self, freq, power, model, priors=None, m=1):
         """
         Posterior distribution for power spectra.
         Uses an exponential distribution for the errors in the likelihood,
@@ -447,16 +467,16 @@ class PSDPosterior(Posterior):
             this module. Note that it is impossible to call the posterior object
             itself or the `self.logposterior` method without defining a prior.
 
+        m: int, default 1
+            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!
+
         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
 
@@ -471,12 +491,11 @@ class PSDPosterior(Posterior):
                for a maximum likelihood-style analysis, no prior is required.
 
         """
-        self.loglikelihood = PSDLogLikelihood(ps.freq,
-                                              ps.power,
-                                              model, m=ps.m)
+        self.loglikelihood = PSDLogLikelihood(freq, power,
+                                              model, m=m)
 
-        self.m = ps.m
-        Posterior.__init__(self, ps.freq, ps.power, model)
+        self.m = m
+        Posterior.__init__(self, freq, power, model)
 
         if not priors is None:
             self.logprior = set_logprior(self, priors)
diff --git a/stingray/modeling/scripts.py b/stingray/modeling/scripts.py
index 8aa5d058..81a7a4a3 100644
--- a/stingray/modeling/scripts.py
+++ b/stingray/modeling/scripts.py
@@ -99,7 +99,7 @@ def fit_powerspectrum(ps, model, starting_pars, max_post=False, priors=None,
     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., 1.5, 2.5]
+    >>> t0 = [80, 1.5, 2.5]
 
     Let's also make a model to test:
     >>> model_to_test = models.PowerLaw1D() + models.Const1D()
@@ -114,10 +114,12 @@ def fit_powerspectrum(ps, model, starting_pars, max_post=False, priors=None,
 
     """
     if priors:
-        lpost = PSDPosterior(ps, model, priors=priors)
+        lpost = PSDPosterior(ps.freq, ps.power, model, priors=priors,
+                             m=ps.m)
     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)
 
@@ -242,4 +244,3 @@ def fit_lorentzians(ps, nlor, starting_pars, fit_whitenoise=True,
 
     return fit_powerspectrum(ps, model, starting_pars, max_post=max_post,
                              priors=priors, fitmethod=fitmethod)
-


Test patch:
diff --git a/stingray/modeling/tests/test_parameterestimation.py b/stingray/modeling/tests/test_parameterestimation.py
index 468666bb..677ff3ce 100644
--- a/stingray/modeling/tests/test_parameterestimation.py
+++ b/stingray/modeling/tests/test_parameterestimation.py
@@ -2,7 +2,7 @@ from __future__ import division, print_function
 import numpy as np
 import scipy.stats
 import os
-import re
+import logging
 
 from astropy.tests.helper import pytest
 from astropy.modeling import models
@@ -12,7 +12,8 @@ from stingray import Powerspectrum
 
 from stingray.modeling import ParameterEstimation, PSDParEst, \
     OptimizationResults, SamplingResults
-from stingray.modeling import PSDPosterior, set_logprior, PSDLogLikelihood
+from stingray.modeling import PSDPosterior, set_logprior, PSDLogLikelihood, \
+    LogLikelihood
 
 try:
     from statsmodels.tools.numdiff import approx_hess
@@ -27,6 +28,13 @@ except ImportError:
     can_sample = False
 
 
+class LogLikelihoodDummy(LogLikelihood):
+    def __init__(self, x, y, model):
+        LogLikelihood.__init__(self, x, y, model)
+
+    def evaluate(self, parse, neg=False):
+        return np.nan
+
 class TestParameterEstimation(object):
 
     @classmethod
@@ -53,7 +61,8 @@ class TestParameterEstimation(object):
             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 = PSDPosterior(cls.ps.freq, cls.ps.power, cls.model,
+                                 m=cls.ps.m)
         cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
 
     def test_par_est_initializes(self):
@@ -69,7 +78,8 @@ class TestParameterEstimation(object):
         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)
+        llike = PSDLogLikelihood(self.ps.freq, self.ps.power,
+                                 self.model, m=self.ps.m)
         pe = ParameterEstimation()
         res = pe.fit(llike, [2.0])
 
@@ -110,6 +120,14 @@ class TestParameterEstimation(object):
         assert isinstance(res,
                           OptimizationResults), "res must be of type OptimizationResults"
 
+    def test_fit_method_fails_with_too_many_tries(self):
+        lpost = LogLikelihoodDummy(self.ps.freq, self.ps.power, self.model)
+        pe = ParameterEstimation()
+        t0 = [2.0]
+
+        with pytest.raises(Exception):
+            res = pe.fit(lpost, t0, neg=True)
+
     def test_compute_lrt_fails_when_garbage_goes_in(self):
         pe = ParameterEstimation()
         t0 = [2.0]
@@ -133,10 +151,12 @@ class TestParameterEstimation(object):
         t0 = [2.0]
         pe = ParameterEstimation()
 
-        delta_deviance = pe.compute_lrt(self.lpost, t0, self.lpost, t0)
+        delta_deviance, opt1, opt2 = pe.compute_lrt(self.lpost, t0,
+                                                    self.lpost, t0)
 
         assert delta_deviance < 1e-7
 
+    @pytest.mark.skipif("not can_sample")
     def test_sampler_runs(self):
         pe = ParameterEstimation()
         if os.path.exists("test_corner.pdf"):
@@ -148,6 +168,19 @@ class TestParameterEstimation(object):
         assert sample_res.acceptance > 0.25
         assert isinstance(sample_res, SamplingResults)
 
+    @pytest.mark.skipif("can_sample")
+    def test_sample_raises_error_without_emcee(self):
+        pe = ParameterEstimation()
+
+        with pytest.raises(ImportError):
+            sample_res = pe.sample(self.lpost, [2.0])
+
+
+    def test_simulate_lrt_fails_in_superclass(self):
+
+        pe = ParameterEstimation()
+        with pytest.raises(NotImplementedError):
+            pe.simulate_lrts(None, None, None, None, None)
 
 class TestOptimizationResults(object):
 
@@ -175,7 +208,8 @@ class TestOptimizationResults(object):
             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 = PSDPosterior(cls.ps.freq, cls.ps.power,
+                                 cls.model, m=cls.ps.m)
         cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
 
         cls.fitmethod = "BFGS"
@@ -244,13 +278,16 @@ class TestOptimizationResults(object):
 
         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
+        if res is not None:
+            self.result = res.fun
+            self.p_opt = res.x
+        else:
+            self.result = None
+            self.p_opt = None
         self.model = lpost.model
 
 
@@ -304,7 +341,8 @@ class TestOptimizationResultInternalFunctions(object):
                       "amplitude_0": p_amplitude_0,
                       "alpha_0": p_alpha_0}
 
-        cls.lpost = PSDPosterior(cls.ps, cls.model)
+        cls.lpost = PSDPosterior(cls.ps.freq, cls.ps.power,
+                                 cls.model, m=cls.ps.m)
         cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
 
         cls.fitmethod = "BFGS"
@@ -318,16 +356,6 @@ class TestOptimizationResultInternalFunctions(object):
         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)
@@ -343,6 +371,15 @@ class TestOptimizationResultInternalFunctions(object):
 
         assert np.all(optres.mfit == mfit_test)
 
+    @pytest.mark.skipif("comp_hessian")
+    def test_compute_covariance_without_comp_hessian(self):
+        optres = OptimizationResultsSubclassDummy(self.lpost, None,
+                                                  neg=True)
+
+        optres._compute_covariance(self.lpost, None)
+        assert optres.cov is None
+        assert optres.err is None
+
     def test_compute_statistics_computes_mfit(self):
         optres = OptimizationResultsSubclassDummy(self.lpost, self.opt,
                                                   neg=True)
@@ -432,6 +469,40 @@ class TestOptimizationResultInternalFunctions(object):
 
             assert np.all(optres.cov == hess_inv)
             assert np.all(optres.err == np.sqrt(np.diag(np.abs(hess_inv))))
+        else:
+            assert optres.cov is None
+            assert optres.err is None
+
+    def test_print_summary_works(self, logger, caplog):
+
+        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, self.opt)
+
+        optres.print_summary(self.lpost)
+
+        assert 'Parameter amplitude_0' in caplog.text
+        assert "Parameter x_0_0" in caplog.text
+        assert "Parameter alpha_0" in caplog.text
+        assert "Parameter amplitude_1" in caplog.text
+        assert "(Fixed)" in caplog.text
+        assert "Fitting statistics" in caplog.text
+        assert "number of data points" in caplog.text
+        assert "Deviance [-2 log L] D =" in caplog.text
+        assert "The Akaike Information Criterion of " \
+               "the model is" in caplog.text
+        assert "The Bayesian Information Criterion of " \
+               "the model is" in caplog.text
+        assert "The figure-of-merit function for this model" in caplog.text
+        assert "Summed Residuals S =" in caplog.text
+        assert "Expected S" in caplog.text
+        assert "merit function" in caplog.text
 
 
 if can_sample:
@@ -477,7 +548,8 @@ if can_sample:
                     amplitude)
 
             cls.priors = {"amplitude": p_amplitude}
-            cls.lpost = PSDPosterior(cls.ps, cls.model)
+            cls.lpost = PSDPosterior(cls.ps.freq, cls.ps.power,
+                                     cls.model, m=cls.ps.m)
             cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
 
             cls.fitmethod = "BFGS"
@@ -502,6 +574,10 @@ if can_sample:
 
             _, _, _ = cls.sampler.run_mcmc(p0, cls.niter)
 
+
+        def test_can_sample_is_true(self):
+            assert can_sample
+
         def test_sample_results_object_initializes(self):
             SamplingResults(self.sampler)
 
@@ -550,6 +626,14 @@ if can_sample:
             assert np.all(np.isclose(test_ci, s.ci, atol=0.01, rtol=0.01))
 
 
+@pytest.fixture()
+def logger():
+    logger = logging.getLogger('Some.Logger')
+    logger.setLevel(logging.INFO)
+
+    return logger
+
+
 class TestPSDParEst(object):
 
     @classmethod
@@ -559,8 +643,8 @@ class TestPSDParEst(object):
         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)
+        rng = np.random.RandomState(100)  # set the seed for the random number generator
+        noise = rng.exponential(size=nfreq)
 
         cls.model = models.Lorentz1D() + models.Const1D()
 
@@ -574,7 +658,6 @@ class TestPSDParEst(object):
         cls.model.fwhm_0 = cls.fwhm_0
         cls.model.amplitude_0 = cls.amplitude_0
         cls.model.amplitude_1 = cls.amplitude_1
-        cls.model.x_0_0.fixed = True
 
         p = cls.model(freq)
 
@@ -609,7 +692,8 @@ class TestPSDParEst(object):
                       "x_0_0": p_x_0_0,
                       "fwhm_0": p_fwhm_0}
 
-        cls.lpost = PSDPosterior(cls.ps, cls.model)
+        cls.lpost = PSDPosterior(cls.ps.freq, cls.ps.power,
+                                 cls.model, m=cls.ps.m)
         cls.lpost.logprior = set_logprior(cls.lpost, cls.priors)
 
         cls.fitmethod = "BFGS"
@@ -620,7 +704,7 @@ class TestPSDParEst(object):
     def test_fitting_with_ties_and_bounds(self, capsys):
         double_f = lambda model : model.x_0_0 * 2
         model = self.model.copy()
-        model =  self.model + models.Lorentz1D(amplitude=model.amplitude_0,
+        model += models.Lorentz1D(amplitude=model.amplitude_0,
                                    x_0 = model.x_0_0 * 2,
                                    fwhm = model.fwhm_0)
         model.x_0_0 = self.model.x_0_0
@@ -643,28 +727,22 @@ class TestPSDParEst(object):
         ps.df = self.ps.df
         ps.norm = "leahy"
 
-        pe = PSDParEst(ps)
+        pe = PSDParEst(ps, fitmethod="TNC")
         llike = PSDLogLikelihood(ps.freq, ps.power, model)
 
-        true_pars = [self.amplitude_0, self.x_0_0, self.fwhm_0,
+        true_pars = [self.x_0_0, self.fwhm_0,
                      self.amplitude_1,
-                     model.amplitude_2.value, model.x_0_2.value,
+                     model.amplitude_2.value,
                      model.fwhm_2.value]
-        res = pe.fit(llike, true_pars)
 
-        res.print_summary(llike)
-        out, err = capsys.readouterr()
-        assert "100.00000            (Fixed)" in out
-        pattern = \
-            re.compile(r"5\) Parameter x_0_2\s+: [0-9]\.[0-9]{5}\s+\(Tied\)")
-        assert pattern.search(out)
+        res = pe.fit(llike, true_pars, neg=True)
 
         compare_pars = [self.x_0_0, self.fwhm_0,
                         self.amplitude_1,
                         model.amplitude_2.value,
                         model.fwhm_2.value]
 
-        assert np.all(np.isclose(compare_pars, res.p_opt, rtol=0.5))
+        assert np.allclose(compare_pars, res.p_opt, rtol=0.5)
 
     def test_par_est_initializes(self):
         pe = PSDParEst(self.ps)
@@ -703,13 +781,14 @@ class TestPSDParEst(object):
 
     def test_fit_method_works_with_correct_parameter(self):
         pe = PSDParEst(self.ps)
-        lpost = PSDPosterior(self.ps, self.model, self.priors)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, self.priors, m=self.ps.m)
         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 " \
@@ -718,7 +797,8 @@ class TestPSDParEst(object):
     def test_plotfits_leahy(self):
         pe = PSDParEst(self.ps)
         t0 = [2.0, 1, 1, 1]
-        lpost = PSDPosterior(self.ps, self.model, self.priors)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, self.priors, m=self.ps.m)
 
         res = pe.fit(lpost, t0)
 
@@ -730,7 +810,6 @@ class TestPSDParEst(object):
     def test_plotfits_log_leahy(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)
 
@@ -748,7 +827,6 @@ class TestPSDParEst(object):
         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)
 
@@ -765,7 +843,6 @@ class TestPSDParEst(object):
         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)
@@ -784,7 +861,6 @@ class TestPSDParEst(object):
         ps.df = self.ps.df
         ps.norm = "none"
         pe = PSDParEst(ps)
-        lpost = PSDPosterior(self.ps, self.model, self.priors)
 
         res = pe.fit(self.lpost, t0)
 
@@ -801,7 +877,6 @@ class TestPSDParEst(object):
         ps.df = self.ps.df
         ps.norm = "none"
         pe = PSDParEst(ps)
-        lpost = PSDPosterior(self.ps, self.model, self.priors)
 
         t0 = [2.0, 1, 1, 1]
         res = pe.fit(self.lpost, t0)
@@ -826,7 +901,7 @@ class TestPSDParEst(object):
         pe = PSDParEst(self.ps, max_post=True)
 
         assert pe.max_post is True
-        delta_deviance = pe.compute_lrt(self.lpost, t0, self.lpost, t0)
+        delta_deviance, _, _ = pe.compute_lrt(self.lpost, t0, self.lpost, t0)
 
         assert pe.max_post is False
 
@@ -836,17 +911,17 @@ class TestPSDParEst(object):
         pe = PSDParEst(self.ps, max_post=True)
 
         # MB: This is a little too random
-        delta_deviance = pe.compute_lrt(self.lpost, t0,
+        delta_deviance, _, _ = pe.compute_lrt(self.lpost, t0,
                                         self.lpost, t0)
 
         assert np.absolute(delta_deviance) < 1.5e-4
 
+    @pytest.mark.skipif("not can_sample")
     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,
+        sample_res = pe.sample(self.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")
@@ -854,3 +929,459 @@ class TestPSDParEst(object):
         assert sample_res.acceptance > 0.25
         assert isinstance(sample_res, SamplingResults)
 
+    def test_generate_model_data(self):
+        pe = PSDParEst(self.ps)
+
+        m = self.model
+        _fitter_to_model_params(m, self.t0)
+
+        model = m(self.ps.freq)
+
+        pe_model = pe._generate_model(self.lpost, [2.0, 0.1, 100, 2.0])
+
+        assert np.allclose(model, pe_model)
+
+
+    def generate_data_rng_object_works(self):
+        pe = PSDParEst(self.ps)
+
+        sim_data1 = pe._generate_data(self.lpost,
+                                      [2.0, 0.1, 100.0, 2.0],
+                                      seed=1)
+        sim_data2 = pe._generate_data(self.lpost,
+                                      [2.0, 0.1, 100.0, 2.0],
+                                      seed=1)
+
+        assert np.allclose(sim_data1.power, sim_data2.power)
+
+    def test_generate_data_produces_correct_distribution(self):
+        model = models.Const1D()
+
+        model.amplitude = 2.0
+
+        p = model(self.ps.freq)
+
+        seed = 100
+        rng = np.random.RandomState(seed)
+
+        noise = rng.exponential(size=len(p))
+        power = noise*p
+
+        ps = Powerspectrum()
+        ps.freq = self.ps.freq
+        ps.power = power
+        ps.m = 1
+        ps.df = self.ps.freq[1]-self.ps.freq[0]
+        ps.norm = "leahy"
+
+        lpost = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        pe = PSDParEst(ps)
+
+        rng2 = np.random.RandomState(seed)
+        sim_data = pe._generate_data(lpost, [2.0], rng2)
+
+        assert np.allclose(ps.power, sim_data.power)
+
+    def test_generate_model_breaks_with_wrong_input(self):
+
+        pe = PSDParEst(self.ps)
+
+        with pytest.raises(AssertionError):
+            pe_model = pe._generate_model([1, 2, 3, 4], [1, 2, 3, 4])
+
+    def  test_generate_model_breaks_for_wrong_number_of_parameters(self):
+        pe = PSDParEst(self.ps)
+
+        with pytest.raises(AssertionError):
+            pe_model = pe._generate_model(self.lpost, [1, 2, 3])
+
+
+    def test_compute_lrt_works(self):
+
+        m = 1
+        nfreq = 100000
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(100)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        s_all = np.atleast_2d(np.ones(10) * 2.0).T
+
+        model2 = models.PowerLaw1D() + models.Const1D()
+        model2.x_0_0.fixed = True
+        loglike2 = PSDLogLikelihood(ps.freq, ps.power, model2, 1)
+
+        pe = PSDParEst(ps)
+
+        lrt_obs, res1, res2 = pe.compute_lrt(loglike, [2.0], loglike2,
+                                             [2.0, 1.0, 2.0], neg=True)
+        lrt_sim = pe.simulate_lrts(s_all, loglike, [2.0], loglike2,
+                                           [2.0, 1.0, 2.0], max_post=False,
+                                           seed=100)
+
+        assert (lrt_obs > 0.4) and (lrt_obs < 0.6)
+        assert np.all(lrt_sim < 10.0) and np.all(lrt_sim > 0.01)
+
+
+    def test_compute_lrt_fails_with_wrong_input(self):
+        pe = PSDParEst(self.ps)
+        with pytest.raises(AssertionError):
+            lrt_sim = pe.simulate_lrts(np.arange(10), self.lpost, [1, 2, 3, 4],
+                                       [1, 2, 3, 4], [1, 2, 3, 4])
+
+
+    def test_pvalue_calculated_correctly(self):
+        a = [1, 1, 1, 2]
+        obs_val = 1.5
+
+        pe = PSDParEst(self.ps)
+        pval = pe._compute_pvalue(obs_val, a)
+
+        assert np.isclose(pval, 1./len(a))
+
+    def test_calibrate_lrt_fails_without_lpost_objects(self):
+        pe = PSDParEst(self.ps)
+
+        with pytest.raises(TypeError):
+            pval = pe.calibrate_lrt(self.lpost, [1, 2, 3, 4],
+                                    np.arange(10), np.arange(4))
+
+    def test_calibrate_lrt_fails_with_wrong_parameters(self):
+        pe = PSDParEst(self.ps)
+
+        with pytest.raises(ValueError):
+            pval = pe.calibrate_lrt(self.lpost, [1, 2, 3, 4],
+                                    self.lpost, [1, 2, 3])
+
+    def test_calibrate_lrt_works_as_expected(self):
+
+        m = 1
+        nfreq = 100000
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(100)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        s_all = np.atleast_2d(np.ones(10) * 2.0).T
+
+        model2 = models.PowerLaw1D() + models.Const1D()
+        model2.x_0_0.fixed = True
+        loglike2 = PSDLogLikelihood(ps.freq, ps.power, model2, 1)
+
+        pe = PSDParEst(ps)
+
+        pval = pe.calibrate_lrt(loglike, [2.0], loglike2,
+                                [2.0, 1.0, 2.0], sample=s_all,
+                                max_post=False, nsim=10,
+                                seed=100)
+
+        assert pval > 0.001
+
+    def test_calibrate_lrt_works_with_mvn(self):
+
+        m = 1
+        nfreq = 100000
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(100)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        model2 = models.PowerLaw1D() + models.Const1D()
+        model2.x_0_0.fixed = True
+        loglike2 = PSDLogLikelihood(ps.freq, ps.power, model2, 1)
+
+        pe = PSDParEst(ps)
+
+        pval = pe.calibrate_lrt(loglike, [2.0], loglike2,
+                                [2.0, 1.0, 2.0], sample=None,
+                                max_post=False, nsim=10,
+                                seed=100)
+
+        assert pval > 0.001
+
+    @pytest.mark.skipif("not can_sample")
+    def test_calibrate_lrt_works_with_sampling(self):
+        m = 1
+        nfreq = 100000
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(100)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        lpost = PSDPosterior(ps.freq, ps.power, model, m=1)
+
+        p_amplitude_1 = lambda amplitude: \
+            scipy.stats.norm(loc=2.0, scale=1.0).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=self.a2_mean, scale=self.a2_var).pdf(
+                amplitude)
+
+
+        priors = {"amplitude": p_amplitude_1}
+
+        priors2 = {"amplitude_1": p_amplitude_1,
+                      "amplitude_0": p_amplitude_0,
+                      "alpha_0": p_alpha_0}
+
+
+        lpost.logprior = set_logprior(lpost, priors)
+
+        model2 = models.PowerLaw1D() + models.Const1D()
+        model2.x_0_0.fixed = True
+        lpost2 = PSDPosterior(ps.freq, ps.power, model2, 1)
+        lpost2.logprior = set_logprior(lpost2, priors2)
+
+        pe = PSDParEst(ps)
+
+        pval = pe.calibrate_lrt(lpost, [2.0], lpost2,
+                                [2.0, 1.0, 2.0], sample=None,
+                                max_post=True, nsim=10, nwalkers=100,
+                                burnin=100, niter=20,
+                                seed=100)
+
+        assert pval > 0.001
+
+    def test_find_highest_outlier_works_as_expected(self):
+
+        mp_ind = 5
+        max_power = 1000.0
+
+        ps = Powerspectrum()
+        ps.freq = np.arange(10)
+        ps.power = np.ones_like(ps.freq)
+        ps.power[mp_ind] = max_power
+        ps.m = 1
+        ps.df = ps.freq[1]-ps.freq[0]
+        ps.norm = "leahy"
+
+        pe = PSDParEst(ps)
+
+        max_x, max_ind = pe._find_outlier(ps.freq, ps.power, max_power)
+
+        assert np.isclose(max_x, ps.freq[mp_ind])
+        assert max_ind == mp_ind
+
+    def test_compute_highest_outlier_works(self):
+
+        mp_ind = 5
+        max_power = 1000.0
+
+        ps = Powerspectrum()
+        ps.freq = np.arange(10)
+        ps.power = np.ones_like(ps.freq)
+        ps.power[mp_ind] = max_power
+        ps.m = 1
+        ps.df = ps.freq[1]-ps.freq[0]
+        ps.norm = "leahy"
+
+        model = models.Const1D()
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=1.0, scale=1.0).pdf(
+                amplitude)
+
+        priors = {"amplitude": p_amplitude}
+
+        lpost = PSDPosterior(ps.freq, ps.power, model, 1)
+        lpost.logprior = set_logprior(lpost, priors)
+
+        pe = PSDParEst(ps)
+
+        res = pe.fit(lpost, [1.0])
+
+        res.mfit = np.ones_like(ps.freq)
+
+        max_y, max_x, max_ind = pe._compute_highest_outlier(lpost, res)
+
+        assert np.isclose(max_y[0], 2*max_power)
+        assert np.isclose(max_x[0], ps.freq[mp_ind])
+        assert max_ind == mp_ind
+
+    def test_simulate_highest_outlier_works(self):
+        m = 1
+        nfreq = 100000
+        seed = 100
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(seed)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        nsim = 10
+
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        s_all = np.atleast_2d(np.ones(nsim) * 2.0).T
+
+        pe = PSDParEst(ps)
+
+        res = pe.fit(loglike, [2.0], neg=True)
+
+        maxpow_sim = pe.simulate_highest_outlier(s_all, loglike, [2.0],
+                                                 max_post=False, seed=seed)
+
+        assert maxpow_sim.shape[0] == nsim
+        assert np.all(maxpow_sim > 20.00) and np.all(maxpow_sim < 31.0)
+
+
+    def test_calibrate_highest_outlier_works(self):
+        m = 1
+        nfreq = 100000
+        seed = 100
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(seed)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        nsim = 10
+
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        s_all = np.atleast_2d(np.ones(nsim) * 2.0).T
+
+        pe = PSDParEst(ps)
+
+        pval = pe.calibrate_highest_outlier(loglike, [2.0], sample=s_all,
+                                            max_post=False, seed=seed)
+
+        assert pval > 0.001
+
+    def test_calibrate_highest_outlier_works_with_mvn(self):
+        m = 1
+        nfreq = 100000
+        seed = 100
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(seed)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        nsim = 10
+
+        loglike = PSDLogLikelihood(ps.freq, ps.power, model, m=1)
+
+        pe = PSDParEst(ps)
+
+        pval = pe.calibrate_highest_outlier(loglike, [2.0], sample=None,
+                                            max_post=False, seed=seed,
+                                            nsim=nsim)
+
+        assert pval > 0.001
+
+    @pytest.mark.skipif("not can_sample")
+    def test_calibrate_highest_outlier_works_with_sampling(self):
+        m = 1
+        nfreq = 100000
+        seed = 100
+        freq = np.linspace(1, 10, nfreq)
+        rng = np.random.RandomState(seed)
+        noise = rng.exponential(size=nfreq)
+        model = models.Const1D()
+        model.amplitude = 2.0
+        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"
+
+        nsim = 10
+
+        lpost = PSDPosterior(ps.freq, ps.power, model, m=1)
+        p_amplitude = lambda amplitude: \
+            scipy.stats.norm(loc=1.0, scale=1.0).pdf(
+                amplitude)
+
+        priors = {"amplitude": p_amplitude}
+        lpost.logprior = set_logprior(lpost, priors)
+
+        pe = PSDParEst(ps)
+
+        pval = pe.calibrate_highest_outlier(lpost, [2.0], sample=None,
+                                            max_post=True, seed=seed,
+                                            nsim=nsim, niter=20, nwalkers=100,
+                                            burnin=100)
+
+        assert pval > 0.001
diff --git a/stingray/modeling/tests/test_posterior.py b/stingray/modeling/tests/test_posterior.py
index 21572d9c..336262d7 100644
--- a/stingray/modeling/tests/test_posterior.py
+++ b/stingray/modeling/tests/test_posterior.py
@@ -49,7 +49,7 @@ class TestSetPrior(object):
         pl = models.PowerLaw1D()
         pl.x_0.fixed = True
 
-        cls.lpost = PSDPosterior(cls.ps, pl)
+        cls.lpost = PSDPosterior(cls.ps.freq, cls.ps.power, pl, m=cls.ps.m)
 
     def test_set_prior_runs(self):
         p_alpha = lambda alpha: ((-1. <= alpha) & (alpha <= 5.))/6.0
@@ -152,20 +152,23 @@ class TestPSDPosterior(object):
         cls.priors = {"amplitude":p_amplitude}
 
     def test_logprior_fails_without_prior(self):
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power, self.model,
+                             m=self.ps.m)
 
         with pytest.raises(AttributeError):
             lpost.logprior([1])
 
     def test_making_posterior(self):
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         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 = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         with pytest.raises(IncorrectParameterError):
@@ -174,7 +177,8 @@ class TestPSDPosterior(object):
     def test_logprior(self):
         t0 = [2.0]
 
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         lp_test = lpost.logprior(t0)
@@ -188,7 +192,8 @@ class TestPSDPosterior(object):
 
         loglike = -np.sum(np.log(mean_model)) - np.sum(self.ps.power/mean_model)
 
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         loglike_test = lpost.loglikelihood(t0, neg=False)
@@ -201,7 +206,8 @@ class TestPSDPosterior(object):
         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 = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         loglike_test = lpost.loglikelihood(t0, neg=True)
@@ -211,7 +217,8 @@ class TestPSDPosterior(object):
     def test_posterior(self):
         t0 = [2.0]
         m = self.model(self.ps.freq[1:], t0)
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         post_test = lpost(t0, neg=False)
@@ -225,7 +232,8 @@ class TestPSDPosterior(object):
     def test_negative_posterior(self):
         t0 = [2.0]
         m = self.model(self.ps.freq[1:], t0)
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         post_test = lpost(t0, neg=True)
@@ -518,20 +526,23 @@ class TestPerPosteriorAveragedPeriodogram(object):
         cls.priors = {"amplitude":p_amplitude}
 
     def test_logprior_fails_without_prior(self):
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
 
         with pytest.raises(AttributeError):
             lpost.logprior([1])
 
     def test_making_posterior(self):
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         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 = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         with pytest.raises(IncorrectParameterError):
@@ -540,7 +551,8 @@ class TestPerPosteriorAveragedPeriodogram(object):
     def test_logprior(self):
         t0 = [2.0]
 
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         lp_test = lpost.logprior(t0)
@@ -557,7 +569,8 @@ class TestPerPosteriorAveragedPeriodogram(object):
                                np.sum((2.0 / (2. * self.m) - 1.0) *
                                       np.log(self.ps.power)))
 
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         loglike_test = lpost.loglikelihood(t0, neg=False)
@@ -577,7 +590,8 @@ class TestPerPosteriorAveragedPeriodogram(object):
                                       np.log(self.ps.power)))
 
 
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         loglike_test = lpost.loglikelihood(t0, neg=True)
@@ -589,7 +603,8 @@ class TestPerPosteriorAveragedPeriodogram(object):
         self.model.amplitude = t0[0]
 
         mean_model = self.model(self.ps.freq)
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         post_test = lpost(t0, neg=False)
@@ -609,7 +624,8 @@ class TestPerPosteriorAveragedPeriodogram(object):
         self.model.amplitude = t0[0]
 
         mean_model = self.model(self.ps.freq)
-        lpost = PSDPosterior(self.ps, self.model)
+        lpost = PSDPosterior(self.ps.freq, self.ps.power,
+                             self.model, m=self.ps.m)
         lpost.logprior = set_logprior(lpost, self.priors)
 
         post_test = lpost(t0, neg=True)
@@ -632,7 +648,7 @@ class TestPerPosteriorAveragedPeriodogram(object):
 
         t0 = [2.0]
         m = self.model(self.ps.freq[1:], t0)
-        lpost = PSDPosterior(ps_nan, self.model)
+        lpost = PSDPosterior(ps_nan.freq, ps_nan.power, 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
index 9904066e..cb13dcec 100644
--- a/stingray/modeling/tests/test_scripts.py
+++ b/stingray/modeling/tests/test_scripts.py
@@ -78,9 +78,14 @@ class TestFitLorentzians(object):
         model.amplitude_0 = self.t0[0]
         # model.bounds = {}
 
+        t0 = np.array([self.amplitude_0, self.x_0_0, self.fwhm_0,
+              self.amplitude_1, self.fwhm_1,
+              self.amplitude_2, self.fwhm_2,
+              self.whitenoise])
+
         parest, res = fit_powerspectrum(self.ps, model,
-                                np.random.normal(self.t0,
-                                                 self.t0 / 10))
+                                np.random.normal(t0,
+                                                 t0 / 10))
 
         true_pars = [self.amplitude_0,
                      self.x_0_0, self.fwhm_0,
@@ -97,9 +102,15 @@ class TestFitLorentzians(object):
         model.amplitude_0.fixed = True
         # model.bounds = {}
 
+        t0 = np.array([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])
+
+
         parest, res = fit_powerspectrum(self.ps, model,
-                                  np.random.normal(self.t0,
-                                                   self.t0 / 10))
+                                  np.random.normal(t0,
+                                                   t0 / 10))
 
         true_pars = [self.x_0_0, self.fwhm_0,
                      self.amplitude_1, self.x_0_1, self.fwhm_1,