How to use the mc3.utils.burn function in mc3

def test_burn_dict():
    Zdict = {'posterior':Z, 'zchain':zchain, 'burnin':1}
    posterior, chain, mask = mu.burn(Zdict)
    np.testing.assert_equal(posterior, np.array([[11.,12.,21.,22.,31.,32.]]).T)
    np.testing.assert_equal(chain, np.array([0, 0, 1, 1, 2, 2]))
    np.testing.assert_equal(mask,  np.array([5, 8, 7, 9, 6, 10]))

def test_burn_Z():
    burnin = 1
    posterior, chain, mask = mu.burn(Z=Z, zchain=zchain, burnin=burnin)
    np.testing.assert_equal(posterior, np.array([[11.,12.,21.,22.,31.,32.]]).T)
    np.testing.assert_equal(chain, np.array([0, 0, 1, 1, 2, 2]))
    np.testing.assert_equal(mask,  np.array([5, 8, 7, 9, 6, 10]))

def test_burn_Z_unburn():
    # Only remove pre-MCMC samples (zchain==-1):
    burnin = 0
    posterior, chain, mask = mu.burn(Z=Z, zchain=zchain, burnin=burnin)
    np.testing.assert_equal(posterior,
        np.array([[10., 11., 12., 20., 21., 22., 30., 31., 32.]]).T)
    np.testing.assert_equal(chain, np.array([0, 0, 0, 1, 1, 1, 2, 2, 2]))
    np.testing.assert_equal(mask,  np.array([2, 5, 8, 3, 7, 9, 4, 6, 10]))

def test_burn_unsort():
    Zdict = {'posterior':Z, 'zchain':zchain, 'burnin':1}
    posterior, chain, mask = mu.burn(Zdict, sort=False)
    np.testing.assert_equal(posterior, np.array([[11.,31.,21.,12.,22.,32.]]).T)
    np.testing.assert_equal(chain, np.array([0, 2, 1, 0, 1, 2]))
    np.testing.assert_equal(mask,  np.array([5, 6, 7, 8, 9, 10]))

def test_burn_override_burnin():
    Zdict = {'posterior':Z, 'zchain':zchain, 'burnin':1}
    posterior, chain, mask = mu.burn(Zdict, burnin=0)
    np.testing.assert_equal(posterior,
        np.array([[10., 11., 12., 20., 21., 22., 30., 31., 32.]]).T)
    np.testing.assert_equal(chain, np.array([0, 0, 0, 1, 1, 1, 2, 2, 2]))
    np.testing.assert_equal(mask,  np.array([2, 5, 8, 3, 7, 9, 4, 6, 10]))

fitpars[ifree] = np.copy(bestp[ifree])
  for s in ishare:
      fitpars[s] = fitpars[-int(pstep[s])-1]
  best_model = chains[0].eval_model(fitpars)

  # Remove pre-MCMC and post-MCMC alocated samples:
  zvalid = zchain>=0
  Z = Z[zvalid]
  zchain = zchain[zvalid]
  log_post = log_post[zvalid]
  log_prior = ms.log_prior(Z, prior, priorlow, priorup, pstep)
  chisq = -2*(log_post - log_prior)
  best_log_prior = ms.log_prior(bestp[ifree], prior, priorlow, priorup, pstep)
  best_chisq = -2*(best_log_post.value - best_log_prior)
  # And remove burn-in samples:
  posterior, _, zmask = mu.burn(Z=Z, zchain=zchain, burnin=zburn)

  # Number of evaluated and kept samples:
  nsample  = len(Z)*thinning
  nzsample = len(posterior)

  # Print out Summary:
  log.msg('\nMCMC Summary:'
          '\n-------------')
  fmt = len(str(nsample))
  log.msg("Number of evaluated samples:        {:{}d}".
          format(nsample,  fmt), indent=2)
  log.msg("Number of parallel chains:          {:{}d}".
          format(nchains,  fmt), indent=2)
  log.msg("Average iterations per chain:       {:{}d}".
          format(nsample//nchains, fmt), indent=2)
  log.msg("Burned-in iterations per chain:     {:{}d}".

if leastsq is not None:
      if (output['best_log_post'] - fit_output['best_log_post'] > 5.0e-8
          and np.any((output['bestp'] - fit_output['bestp'])!=0.0)):
          log.warning("MCMC found a better fit than the minimizer:\n"
              "MCMC best-fitting parameters:        (chisq={:.8g})\n{}\n"
              "Minimizer best-fitting parameters:   (chisq={:.8g})\n{}".
              format(-2*output['best_log_post'], output['bestp'],
                  -2*fit_output['best_log_post'], fit_output['bestp']))
      else:
          output['best_log_post'] = fit_output['best_log_post']
          output['best_chisq'] = fit_output['best_chisq']
          output['best_model'] = fit_output['best_model']
          output['bestp'] = fit_output['bestp']

  # And remove burn-in samples:
  posterior, zchain, zmask = mu.burn(Z=output['posterior'],
      zchain=output['zchain'], burnin=output['burnin'])

  # Get some stats:
  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  = []