<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def mcmc(x1, x2, x1err=[], x2err=[], po=(1,1,0.5), logify=True, nsteps=5000, nwalkers=100, nburn=500, output='full'):
 """ Use emcee to find the best-fit linear relation or power law accounting for measurement uncertainties and intrinsic scatter Parameters x1 : array of floats Independent variable, or observable x2 : array of floats Dependent variable x1err : array of floats (optional) Uncertainties on the independent variable x2err : array of floats (optional) Uncertainties on the dependent variable po : tuple of 3 floats (optional) Initial guesses for zero point, slope, and intrinsic scatter. Results are not very sensitive to these values so they shouldn't matter a lot. logify : bool (default True) Whether to take the log of the measurements in order to estimate the best-fit power law instead of linear relation nsteps : int (default 5000) Number of steps each walker should take in the MCMC nwalkers : int (default 100) Number of MCMC walkers nburn : int (default 500) Number of samples to discard to give the MCMC enough time to converge. output : list of ints or 'full' (default 'full') If 'full', then return the full samples (except for burn-in section) for each parameter. Otherwise, each float corresponds to a percentile that will be returned for each parameter. Returns ------- See *output* argument above for return options. """

    import emcee
    if len(x1err) == 0:
        x1err = numpy.ones(len(x1))
    if len(x2err) == 0:
        x2err = numpy.ones(len(x1))
    def lnlike(theta, x, y, xerr, yerr):
        a, b, s = theta
        model = a + b*x
        sigma = numpy.sqrt((b*xerr)**2 + yerr*2 + s**2)
        lglk = 2 * sum(numpy.log(sigma)) + \
               sum(((y-model) / sigma) ** 2) + \
               numpy.log(len(x)) * numpy.sqrt(2*numpy.pi) / 2
        return -lglk
    def lnprior(theta):
        a, b, s = theta
        if s >= 0:
            return 0
        return -numpy.inf
    def lnprob(theta, x, y, xerr, yerr):
        lp = lnprior(theta)
        return lp + lnlike(theta, x, y, xerr, yerr)
    if logify:
        x1, x2, x1err, x2err = to_log(x1, x2, x1err, x2err)
    start = numpy.array(po)
    ndim = len(start)
    pos = [start + 1e-4*numpy.random.randn(ndim) for i in range(nwalkers)]
    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob,
                                    args=(x1,x2,x1err,x2err))
    sampler.run_mcmc(pos, nsteps)
    samples = numpy.array([sampler.chain[:,nburn:,i].reshape(-1) \
                           for i in xrange(ndim)])
    if logify:
        samples[2] *= numpy.log(10)
    if output == 'full':
        return samples
    else:
        try:
            values = [[numpy.percentile(s, o) for o in output] 
                      for s in samples]
            return values
        except TypeError:
            msg = 'ERROR: wrong value for argument output in mcmc().'
            msg += ' Must be "full" or list of ints.'
            print msg
            exit()
    return