emcee_mcmc_linear_regression_parameter_estimation.py

1import numpy as np
2import emcee
3
4def log_prob(theta, x, y, yerr):
5    m, b, log_f = theta
6    model = m * x + b
7    sigma2 = yerr**2 + model**2 * np.exp(2 * log_f)
8    return -0.5 * np.sum((y - model)**2 / sigma2 + np.log(sigma2))
9
10# Generate some synthetic data from a model.
11np.random.seed(123)
12m_true = -0.9594
13b_true = 4.294
14f_true = 0.534
15x = np.sort(10 * np.random.rand(50))
16yerr = 0.1 + 0.5 * np.random.rand(50)
17y = m_true * x + b_true
18y += np.abs(f_true * y) * np.random.randn(50)
19y += yerr * np.random.randn(50)
20
21# Initialize the walkers
22nwalkers = 32
23ndim = 3
24pos = np.array([m_true, b_true, np.log(f_true)]) + 1e-4 * np.random.randn(nwalkers, ndim)
25
26# Run the sampler
27sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=(x, y, yerr))
28sampler.run_mcmc(pos, 5000, progress=True)
29
30# Access the results
31samples = sampler.get_chain()