docs: 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