Fit transit signal to trapezoid model using MCMC

        .. note:: As currently implemented, this method creates a
            bunch of attributes relevant to the MCMC fit; I plan
            to refactor this to define those attributes as properties
            so as not to have their creation hidden away here.  I plan
            to refactor how this works.

def MCMC(self, niter=500, nburn=200, nwalkers=200, threads=1,
             fit_partial=False, width=3, savedir=None, refit=False,
             thin=10, conf=0.95, maxslope=MAXSLOPE, debug=False, p0=None):
        """
        Fit transit signal to trapezoid model using MCMC

        .. note:: As currently implemented, this method creates a
            bunch of attributes relevant to the MCMC fit; I plan
            to refactor this to define those attributes as properties
            so as not to have their creation hidden away here.  I plan
            to refactor how this works.
        """
        if fit_partial:
            wok = np.where((np.absolute(self.ts-self.center) < (width*self.dur)) &
                           ~np.isnan(self.fs))
        else:
            wok = np.where(~np.isnan(self.fs))

        if savedir is not None:
            if not os.path.exists(savedir):
                os.mkdir(savedir)

        alreadydone = True
        alreadydone &= savedir is not None
        alreadydone &= os.path.exists('%s/ts.npy' % savedir)
        alreadydone &= os.path.exists('%s/fs.npy' % savedir)

        if savedir is not None and alreadydone:
            ts_done = np.load('%s/ts.npy' % savedir)
            fs_done = np.load('%s/fs.npy' % savedir)
            alreadydone &= np.all(ts_done == self.ts[wok])
            alreadydone &= np.all(fs_done == self.fs[wok])

        if alreadydone and not refit:
            logging.info('MCMC fit already done for %s.  Loading chains.' % self.name)
            Ts = np.load('%s/duration_chain.npy' % savedir)
            ds = np.load('%s/depth_chain.npy' % savedir)
            slopes = np.load('%s/slope_chain.npy' % savedir)
            tcs = np.load('%s/tc_chain.npy' % savedir)
        else:
            logging.info('Fitting data to trapezoid shape with MCMC for %s....' % self.name)
            if p0 is None:
                p0 = self.trapfit.copy()
                p0[0] = np.absolute(p0[0])
                if p0[2] < 2:
                    p0[2] = 2.01
                if p0[1] < 0:
                    p0[1] = 1e-5
            logging.debug('p0 for MCMC = {}'.format(p0))
            sampler = traptransit_MCMC(self.ts[wok],self.fs[wok],self.dfs[wok],
                                        niter=niter,nburn=nburn,nwalkers=nwalkers,
                                        threads=threads,p0=p0,return_sampler=True,
                                        maxslope=maxslope)

            Ts,ds,slopes,tcs = (sampler.flatchain[:,0],sampler.flatchain[:,1],
                                sampler.flatchain[:,2],sampler.flatchain[:,3])

            self.sampler = sampler
            if savedir is not None:
                np.save('%s/duration_chain.npy' % savedir,Ts)
                np.save('%s/depth_chain.npy' % savedir,ds)
                np.save('%s/slope_chain.npy' % savedir,slopes)
                np.save('%s/tc_chain.npy' % savedir,tcs)
                np.save('%s/ts.npy' % savedir,self.ts[wok])
                np.save('%s/fs.npy' % savedir,self.fs[wok])

        if debug:
            print(Ts)
            print(ds)
            print(slopes)
            print(tcs)

        N = len(Ts)
        try:
            self.Ts_acor = integrated_time(Ts)
            self.ds_acor = integrated_time(ds)
            self.slopes_acor = integrated_time(slopes)
            self.tcs_acor = integrated_time(tcs)
            self.fit_converged = True
        except AutocorrError:
            self.fit_converged = False


        ok = (Ts > 0) & (ds > 0) & (slopes > 0) & (slopes < self.maxslope)
        logging.debug('trapezoidal fit has {} good sample points'.format(ok.sum()))
        if ok.sum()==0:
            if (Ts > 0).sum()==0:
                #logging.debug('{} points with Ts > 0'.format((Ts > 0).sum()))
                logging.debug('{}'.format(Ts))
                raise MCMCError('{}: 0 points with Ts > 0'.format(self.name))
            if (ds > 0).sum()==0:
                #logging.debug('{} points with ds > 0'.format((ds > 0).sum()))
                logging.debug('{}'.format(ds))
                raise MCMCError('{}: 0 points with ds > 0'.format(self.name))
            if (slopes > 0).sum()==0:
                #logging.debug('{} points with slopes > 0'.format((slopes > 0).sum()))
                logging.debug('{}'.format(slopes))
                raise MCMCError('{}: 0 points with slopes > 0'.format(self.name))
            if (slopes < self.maxslope).sum()==0:
                #logging.debug('{} points with slopes < maxslope ({})'.format((slopes < self.maxslope).sum(),self.maxslope))
                logging.debug('{}'.format(slopes))
                raise MCMCError('{} points with slopes < maxslope ({})'.format((slopes < self.maxslope).sum(),self.maxslope))


        durs,deps,logdeps,slopes = (Ts[ok],ds[ok],np.log10(ds[ok]),
                                              slopes[ok])


        inds = (np.arange(len(durs)/thin)*thin).astype(int)
        durs,deps,logdeps,slopes = (durs[inds],deps[inds],logdeps[inds],
                                              slopes[inds])

        self.durs,self.deps,self.logdeps,self.slopes = (durs,deps,logdeps,slopes)

        self._make_kde(conf=conf)

        self.hasMCMC = True