Python里面mcmc

MCMC（Markov Chain Monte Carlo）是一种基于马尔可夫链的随机模拟方法。在Python中，可以使用PyMC3和emcee等库来实现MCMC。

PyMC3是一个基于Theano的概率编程库，可以用于贝叶斯建模和MCMC模拟。它提供了灵活的模型定义和推断方法，并支持多种MCMC算法（如NUTS、Metropolis-Hastings等）。

以下是使用PyMC3进行MCMC的示例代码：

import pymc3 as pm
import numpy as np

# 构造数据
np.random.seed(123)
x = np.linspace(0, 5, 100)
y = 2 * x + np.random.normal(0, 1, size=100)

# 定义模型
with pm.Model() as model:
    # 定义先验分布
    beta = pm.Normal('beta', mu=0, sd=10)
    sigma = pm.HalfNormal('sigma', sd=1)
    
    # 定义似然函数
    y_obs = pm.Normal('y_obs', mu=beta * x, sd=sigma, observed=y)
    
    # 进行MCMC推断
    trace = pm.sample(1000, tune=1000)
    
# 分析结果
pm.traceplot(trace)
pm.summary(trace)

emcee是一个基于MCMC的贝叶斯推断库，使用Python实现。它提供了高效的并行计算和自适应步长控制等功能。

以下是使用emcee进行MCMC的示例代码：

import emcee
import numpy as np

# 构造数据
np.random.seed(123)
x = np.linspace(0, 5, 100)
y = 2 * x + np.random.normal(0, 1, size=100)

# 定义模型
def log_likelihood(theta, x, y):
    beta, ln_sigma = theta
    sigma = np.exp(ln_sigma)
    y_model = beta * x
    return -0.5 * np.sum((y - y_model)**2 / sigma**2 + np.log(2 * np.pi * sigma**2))

def log_prior(theta):
    beta, ln_sigma = theta
    if -10 < beta < 10 and -10 < ln_sigma < 10:
        return 0
    return -np.inf

def log_posterior(theta, x, y):
    return log_prior(theta) + log_likelihood(theta, x, y)

# 进行MCMC推断
ndim = 2
nwalkers = 100
pos = np.random.randn(nwalkers, ndim)
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x, y))
sampler.run_mcmc(pos, 1000, progress=True)

# 分析结果
import corner
samples = sampler.get_chain(discard=100, flat=True)
corner.corner(samples, labels=['beta', 'ln_sigma'])

在以上两个示例中，我们分别使用了PyMC3和emcee进行了MCMC推断，并得到了贝叶斯估计的后验分布。通过分析后验分布，我们可以得到模型参数的置信区间和概率密度估计等信息，从而对模型进行进一步的推断和分析。