https://qiita.com/niikura/items/79dc6837f017c05afaa7 を参考に lmfit + emcee でフィットした場合にどうなるか検証してみました

test_mcmc_emcee.py

#!/usr/bin/env python

import numpy as np
import lmfit
import matplotlib.pyplot as plt

# reference: https://qiita.com/niikura/items/79dc6837f017c05afaa7
data = np.loadtxt("data.txt")
xx = data.T[0]
yy = data.T[1]
ey = data.T[2]

p = lmfit.Parameters()
p.add_many(("m", 0.5), ("b", -0.5))


def log_likelihood(p):
    v = p.valuesdict()
    model = v["m"] * xx + v["b"]
    s2 = ey * ey
    return -0.5 * np.sum((yy - model) ** 2 / s2 + np.log(2 * np.pi * s2))


res = lmfit.minimize(
    log_likelihood,
    method="emcee",
    burn=300,
    steps=2500,
    thin=20,
    params=p,
    is_weighted=True,
)
lmfit.printfuncs.report_fit(res.params)

# plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax.errorbar(xx, yy, ey, fmt="o", label="data")
ax.plot(xx, xx * res.params["m"] + res.params["b"], label="best fit")
ax.legend(loc="upper left")
plt.show()