def lnprob_fun(pars, x, y, sigy):
    """
    The model is a straight line.
    This function returns the natural logarithm of

        P(model | data) ~ P(data | model) P(model)

    In this case we assume:

        1) P(model) = cost (i.e. uninformative prior)
        2) P(data | model) = likelihood for Gaussian errors

    This implies:

        P(model | data) ~ -0.5*chi^2

    """
    ymod = pars[0] + pars[1] * x
    resid = (y - ymod) / sigy
    chi2 = resid @ resid
    return -0.5 * chi2
