How to use the mc3.stats.cred_region function in mc3

def test_cred_region():
    np.random.seed(2)
    posterior = np.random.normal(0, 1.0, 100000)
    pdf, xpdf, HPDmin = ms.cred_region(posterior)
    np.testing.assert_approx_equal(np.amin(xpdf[pdf>HPDmin]), -1.0,
        significant=3)
    np.testing.assert_approx_equal(np.amax(xpdf[pdf>HPDmin]), 1.0,
        significant=3)

def credregion(posterior=None, percentile=0.6827, pdf=None, xpdf=None):
    """
    Compute the highest-posterior-density credible region for a
    posterior distribution.

    This function has been deprecated.  Use mc3.stats.cred_region()
    instead.
    """
    with Log() as log:
        log.warning('Deprecation warning: mc3.utils.credregion() moved to '
                    'mc3.stats.cred_region().')
    from .. import stats as ms
    return ms.cred_region(posterior, percentile, pdf, xpdf)

output['chisq_factor'] = chisq_factor
  output['BIC'] = output['best_chisq'] + nfree*np.log(ndata)
  if ndata > nfree:
      output['red_chisq'] = output['best_chisq']/(ndata-nfree)
  else:
      output['red_chisq'] = np.nan
  output['stddev_residuals'] = np.std(output['best_model']-data)

  # Compute the credible region for each parameter:
  bestp = output['bestp']
  CRlo = np.zeros(nparams)
  CRhi = np.zeros(nparams)
  pdf  = []
  xpdf = []
  for post, idx in zip(posterior.T, ifree):
      PDF, Xpdf, HPDmin = ms.cred_region(post)
      pdf.append(PDF)
      xpdf.append(Xpdf)
      CRlo[idx] = np.amin(Xpdf[PDF>HPDmin])
      CRhi[idx] = np.amax(Xpdf[PDF>HPDmin])
  # CR relative to the best-fitting value:
  CRlo[ifree] -= bestp[ifree]
  CRhi[ifree] -= bestp[ifree]

  # Get the mean and standard deviation from the posterior:
  meanp = np.zeros(nparams, np.double) # Parameters mean
  stdp  = np.zeros(nparams, np.double) # Parameter standard deviation
  meanp[ifree] = np.mean(posterior, axis=0)
  stdp [ifree] = np.std(posterior,  axis=0)
  for s in ishare:
      bestp[s] = bestp[-int(pstep[s])-1]
      meanp[s] = meanp[-int(pstep[s])-1]

npages = 1  # Assume there's only one page
      figs = [axes[0].get_figure()]
      for ax in axes:
          ax.set_yticklabels([])

  maxylim = 0
  for ipar in range(npars):
      ax = axes[ipar]
      ax.tick_params(labelsize=fs-1)
      plt.setp(ax.xaxis.get_majorticklabels(), rotation=90)
      ax.set_xlabel(pnames[ipar], size=fs)
      vals, bins, h = ax.hist(posterior[0::thinning,ipar], bins=25,
          range=ranges[ipar], zorder=0, **hkw)
      # Plot HPD region:
      if quantile is not None:
          PDF, Xpdf, HPDmin = ms.cred_region(posterior[:,ipar], quantile,
                                             pdf[ipar], xpdf[ipar])
          vals = np.r_[0, vals, 0]
          bins = np.r_[bins[0] - (bins[1]-bins[0]), bins]
          # Interpolate xpdf into the histogram:
          f = si.interp1d(bins+0.5*(bins[1]-bins[0]), vals, kind='nearest')
          # Plot the HPD region as shaded areas:
          if ranges[ipar] is not None:
              xran = np.argwhere((Xpdf>ranges[ipar][0])
                               & (Xpdf=HPDmin,
              facecolor='0.75', edgecolor='none', interpolate=False,
              zorder=-2)
      if bestp is not None:
          ax.axvline(bestp[ipar], dashes=(7,4), lw=1.0, **bkw)