def bces(x1, x2, x1err=[], x2err=[], cerr=[], logify=True, model='yx', \
         bootstrap=5000, verbose='normal', full_output=True):
    """
    Bivariate, Correlated Errors and intrinsic Scatter (BCES)
    translated from the FORTRAN code by Christina Bird and Matthew Bershady
    (Akritas & Bershady, 1996)

    Linear regression in the presence of heteroscedastic errors on both
    variables 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
      cerr      : array of floats (optional)
                  Covariances of the uncertainties in the dependent and
                  independent variables
      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
      model     : {'yx', 'xy', 'bi', 'orth'}
                  BCES model with which to calculate regression. See Notes
                  below for details.
      bootstrap : False or int (default 5000)
                  get the errors from bootstrap resampling instead of the
                  analytical prescription? if bootstrap is an int, it is the
                  number of bootstrap resamplings
      verbose   : str (default 'normal')
                  Verbose level. Options are {'quiet', 'normal', 'debug'}
      full_output : bool (default True)
                  If True, return also the covariance between the
                  normalization and slope of the regression.

    Returns
    -------
      a         : tuple of length 2
                  Best-fit normalization and its uncertainty (a, da)
      b         : tuple of length 2
                  Best-fit slope and its uncertainty (b, db)

    Optional outputs
    ----------------
      cov_ab    : 2x2 array of floats
                  covariance between a and b. Returned if full_output is set to
                  True.

    Notes
    -----
      If verbose is normal or debug, the results from all the BCES models will
      be printed (still, only the one selected in *model* will be returned).

      the *model* parameter:
        -'yx' stands for BCES(Y|X)
        -'xy' stands for BCES(X|Y)
        -'bi' stands for BCES Bisector
        -'orth' stands for BCES Orthogonal

    """

    def _bess_bootstrap(npts, x1, x2, x1err, x2err, cerr,nsim):
        ##added by Gerrit, July 2014
        ##Unfortunately I needed a copy of the _bess function for bootstrapping.
        #Would be nicer if those two could be combined
        """
        Do the entire regression calculation for 4 slopes:
        OLS(Y|X), OLS(X|Y), bisector, orthogonal
        """
        #calculate sigma's for datapoints using length of confidence intervals
        sig11var = numpy.sum(x1err ** 2,axis=1,keepdims=True) / npts
        sig22var = numpy.sum(x2err ** 2,axis=1,keepdims=True) / npts
        sig12var = numpy.sum(cerr,axis=1,keepdims=True) / npts
        
        # calculate means and variances
        x1av = numpy.mean(x1,axis=1,keepdims=True)
        x1var = x1.var(axis=1,keepdims=True)
        x2av = numpy.mean(x2,axis=1,keepdims=True)
        x2var = x2.var(axis=1,keepdims=True)
        covar_x1x2 = numpy.mean((x1-numpy.mean(x1,axis=1,keepdims=True)) * \
                             (x2-numpy.mean(x2,axis=1,keepdims=True)),
                             axis=1,keepdims=True)
        
        # compute the regression slopes for OLS(X2|X1), OLS(X1|X2), 
        # bisector and orthogonal
        if model == 'yx':
            modelint = 1
        else:
            modelint = 4
        b = numpy.zeros((modelint,nsim))
        b[0] = ((covar_x1x2 - sig12var) / (x1var - sig11var)).flatten()
        if model != 'yx':
            b[1] = ((x2var - sig22var) / (covar_x1x2 - sig12var)).flatten()
            b[2] = ((b[0] * b[1] - 1 + numpy.sqrt((1 + b[0] ** 2) * \
                   (1 + b[1] ** 2))) / (b[0] + b[1])).flatten()
            b[3] = 0.5 * ((b[1] - 1 / b[0]) + numpy.sign(covar_x1x2).flatten()* \
                   numpy.sqrt(4 + (b[1] - 1 / b[0]) ** 2))
        
        # compute intercepts for above 4 cases:
        a = x2av.flatten() - b * x1av.flatten()
        
        # set up variables to calculate standard deviations of slope and 
        # intercept
        xi = []
        xi.append(((x1 - x1av) * (x2 - b[0].reshape(nsim,1) * x1 - \
                                  a[0].reshape(nsim,1)) + \
                   b[0].reshape(nsim,1) * x1err ** 2) / \
                  (x1var - sig11var))
        if model != 'yx':
            xi.append(((x2 - x2av) * (x2 - b[1].reshape(nsim,1) * x1 - \
                                      a[1].reshape(nsim,1)) + x2err ** 2) / \
                      covar_x1x2)
            xi.append((xi[0] * (1 + b[1].reshape(nsim,1) ** 2) + \
                       xi[1] * (1 + b[0].reshape(nsim,1) ** 2)) / \
                      ((b[0].reshape(nsim,1) + \
                       b[1].reshape(nsim,1)) * \
                       numpy.sqrt((1 + b[0].reshape(nsim,1) ** 2) * \
                               (1 + b[1].reshape(nsim,1) ** 2))))
            xi.append((xi[0] / b[0].reshape(nsim,1) ** 2 + xi[1]) * \
                      b[3].reshape(nsim,1) / \
                      numpy.sqrt(4 + (b[1].reshape(nsim,1) - \
                              1 / b[0].reshape(nsim,1)) ** 2))
        zeta = []
        for i in xrange(modelint):
            zeta.append(x2 - b[i].reshape(nsim,1) * x1 - x1av * xi[i])
        
        # calculate  variance for all a and b
        bvar = numpy.zeros((4,nsim))
        avar = numpy.zeros((4,nsim))
        for i in xrange(modelint):
            bvar[i] = xi[i].var(axis=1,keepdims=False)/ npts
            avar[i] = zeta[i].var(axis=1,keepdims=False) / npts
        return a, b, avar, bvar, xi, zeta

    def _bess(npts, x1, x2, x1err, x2err, cerr):
        """
        Do the entire regression calculation for 4 slopes:
          OLS(Y|X), OLS(X|Y), bisector, orthogonal
        """
        # calculate sigma's for datapoints using length of confidence
        # intervals
        sig11var = sum(x1err ** 2) / npts
        sig22var = sum(x2err ** 2) / npts
        sig12var = sum(cerr) / npts
        # calculate means and variances
        x1av = numpy.average(x1)
        x1var = numpy.std(x1) ** 2
        x2av = numpy.average(x2)
        x2var = numpy.std(x2) ** 2
        covar_x1x2 = sum((x1 - x1av) * (x2 - x2av)) / npts
        # compute the regression slopes for OLS(X2|X1), OLS(X1|X2), 
        # bisector and orthogonal
        b = numpy.zeros(4)
        b[0] = (covar_x1x2 - sig12var) / (x1var - sig11var)
        b[1] = (x2var - sig22var) / (covar_x1x2 - sig12var)
        b[2] = (b[0] * b[1] - 1 + numpy.sqrt((1 + b[0] ** 2) * \
               (1 + b[1] ** 2))) / (b[0] + b[1])
        b[3] = 0.5 * ((b[1] - 1 / b[0]) + numpy.sign(covar_x1x2) * \
               numpy.sqrt(4 + (b[1] - 1 / b[0]) ** 2))
        # compute intercepts for above 4 cases:
        a = x2av - b * x1av
        # set up variables to calculate standard deviations of slope
        # and intercept
        xi = []
        xi.append(((x1 - x1av) * \
                   (x2 - b[0] * x1 - a[0]) + b[0] * x1err ** 2) / \
                  (x1var - sig11var))
        xi.append(((x2 - x2av) * (x2 - b[1] * x1 - a[1]) + x2err ** 2) / \
                  covar_x1x2)
        xi.append((xi[0] * (1 + b[1] ** 2) + xi[1] * (1 + b[0] ** 2)) / \
                  ((b[0] + b[1]) * \
                   numpy.sqrt((1 + b[0] ** 2) * (1 + b[1] ** 2))))
        xi.append((xi[0] / b[0] ** 2 + xi[1]) * b[3] / \
                  numpy.sqrt(4 + (b[1] - 1 / b[0]) ** 2))
        zeta = []
        for i in xrange(4):
            zeta.append(x2 - b[i]*x1 - x1av*xi[i])
        # calculate  variance for all a and b
        bvar = numpy.zeros(4)
        avar = numpy.zeros(4)
        for i in xrange(4):
            bvar[i] = numpy.std(xi[i]) ** 2 / npts
            avar[i] = numpy.std(zeta[i]) ** 2 / npts
        return a, b, avar, bvar, xi, zeta

    def _bootspbec(npts, x, y, xerr, yerr, cerr):
        """
        Bootstrap samples
        """
        j = numpy.random.randint(npts, size = npts)
        xboot = x[j]
        xerrboot = xerr[j]
        yboot = y[j]
        yerrboot = yerr[j]
        cerrboot = cerr[j]
        return xboot, yboot, xerrboot, yerrboot, cerrboot

    # ----  Main routine starts here  ---- #
    # convert to numpy arrays just in case
    x1 = numpy.array(x1)
    x2 = numpy.array(x2)
    x1err = numpy.array(x1err)
    x2err = numpy.array(x2err)
    if logify:
        x1, x2, x1err, x2err = to_log(x1, x2, x1err, x2err)
    cerr = numpy.array(cerr)
    models = [['yx', 'xy', 'bi', 'orth'],
              ['BCES(Y|X)', 'BCES(X|Y)', 'BCES Bisector', 'BCES Orthogonal']]
    # which to return?
    j = models[0].index(model)
    npts = len(x1)
    # are the errors defined?
    if len(x1err) == 0:
        x1err = numpy.zeros(npts)
    if len(x2err) == 0:
        x2err = numpy.zeros(npts)
    if len(cerr) == 0:
        cerr = numpy.zeros(npts)
    if verbose == 'debug':
        print 'x1 =', x1
        print 'x1err =', x1err
        print 'x2 =', x2
        print 'x2err =', x2err
        print 'cerr =', cerr
        print '\n ** Returning values for', models[1][j], '**'
        if bootstrap is not False:
            print '    with errors from %d bootstrap resamplings' %bootstrap
        print ''

    # calculate nominal fits
    bessresults = _bess(npts, x1, x2, x1err, x2err, cerr)
    (a, b, avar, bvar, xi, zeta) = bessresults
    # covariance between normalization and slope
    if full_output:
        covar_ab = numpy.cov(xi[j], zeta[j])

    if bootstrap is not False:
        # make bootstrap simulated datasets, and compute averages and
        # standard deviations of regression coefficients
        asum = numpy.zeros(4)
        assum = numpy.zeros(4)
        bsum = numpy.zeros(4)
        bssum = numpy.zeros(4)
        sda = numpy.zeros(4)
        sdb = numpy.zeros(4)
        for i in xrange(bootstrap):
            samples = _bootspbec(npts, x1, x2, x1err, x2err, cerr)
            (x1sim, x2sim, x1errsim, x2errsim, cerrsim) = samples
            besssim = _bess(npts, x1sim, x2sim, x1errsim, x2errsim, cerrsim)
            (asim, bsim, avarsim, bvarsim, xi, zeta) = besssim
            asum += asim
            assum += asim ** 2
            bsum += bsim
            bssum += bsim ** 2

        aavg = asum / bootstrap
        bavg = bsum / bootstrap
        for i in range(4):
            sdtest = assum[i] - bootstrap * aavg[i] ** 2
            if sdtest > 0:
                sda[i] = numpy.sqrt(sdtest / (bootstrap - 1))
            sdtest = bssum[i] - bootstrap * bavg[i] ** 2
            if sdtest > 0:
                sdb[i] = numpy.sqrt(sdtest / (bootstrap - 1))

    if verbose in ('normal', 'debug'):
        print '%s   B          err(B)' %('Fit'.ljust(19)),
        print '         A          err(A)'
        for i in range(4):
            print '%s  %9.2e +/- %8.2e    %10.3e +/- %9.3e' \
                  %(models[1][i].ljust(16), b[i], 
                    numpy.sqrt(bvar[i]), a[i], numpy.sqrt(avar[i]))
            if bootstrap is not False:
                print '%s  %9.2e +/- %8.2e    %10.3e +/- %9.3e' \
                      %('bootstrap'.ljust(16), bavg[i],
                        sdb[i], aavg[i], sda[i])
            print ''
        if verbose == 'debug':
            print 'cov[%s] =' %models[model]
            print covar_ab

    if bootstrap is not False:
        if full_output:
            return (a[j], sda[j]), (b[j], sdb[j]), covar_ab
        else:
            return (a[j], sda[j]), (b[j], sdb[j])
    if full_output:
        out = ((a[j], numpy.sqrt(avar[j])),
               (b[j], numpy.sqrt(bvar[j])),
               covar_ab)
    else:
        out = ((a[j], numpy.sqrt(avar[j])),
               (b[j], numpy.sqrt(bvar[j])))
    return out
Bivariate, Correlated Errors and intrinsic Scatter (BCES)
    translated from the FORTRAN code by Christina Bird and Matthew Bershady
    (Akritas & Bershady, 1996)

    Linear regression in the presence of heteroscedastic errors on both
    variables 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
      cerr      : array of floats (optional)
                  Covariances of the uncertainties in the dependent and
                  independent variables
      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
      model     : {'yx', 'xy', 'bi', 'orth'}
                  BCES model with which to calculate regression. See Notes
                  below for details.
      bootstrap : False or int (default 5000)
                  get the errors from bootstrap resampling instead of the
                  analytical prescription? if bootstrap is an int, it is the
                  number of bootstrap resamplings
      verbose   : str (default 'normal')
                  Verbose level. Options are {'quiet', 'normal', 'debug'}
      full_output : bool (default True)
                  If True, return also the covariance between the
                  normalization and slope of the regression.

    Returns
    -------
      a         : tuple of length 2
                  Best-fit normalization and its uncertainty (a, da)
      b         : tuple of length 2
                  Best-fit slope and its uncertainty (b, db)

    Optional outputs
    ----------------
      cov_ab    : 2x2 array of floats
                  covariance between a and b. Returned if full_output is set to
                  True.

    Notes
    -----
      If verbose is normal or debug, the results from all the BCES models will
      be printed (still, only the one selected in *model* will be returned).

      the *model* parameter:
        -'yx' stands for BCES(Y|X)
        -'xy' stands for BCES(X|Y)
        -'bi' stands for BCES Bisector
        -'orth' stands for BCES Orthogonal