写一个mcmc拟合程序
时间: 2023-05-27 13:02:54 浏览: 65
由于MCMC算法在不同的问题中应用非常广泛,因此编写一个通用的MCMC程序是比较困难的。但是,以下是一个简单的Python示例程序,用于在一维线性模型中使用MCMC拟合数据。
```python
import numpy as np
import matplotlib.pyplot as plt
# True model parameters
theta_true = [2, 1]
# Generate some synthetic data
N = 50
x = np.linspace(0, 1, N)
y = theta_true[0] + theta_true[1]*x + np.random.normal(0, 0.2, N)
# Define prior distribution
def prior(theta):
return 1
# Define likelihood function
def likelihood(theta, x, y):
return np.prod(1/np.sqrt(2*np.pi*0.2**2)*np.exp(-(y - (theta[0] + theta[1]*x))**2/(2*0.2**2)))
# Define posterior distribution
def posterior(theta, x, y):
return likelihood(theta, x, y)*prior(theta)
# Define MCMC algorithm
def mcmc(posterior, x, y, theta_0, num_samples):
theta_current = theta_0
samples = []
for i in range(num_samples):
# Sample a new theta from proposal distribution
theta_proposal = theta_current + np.random.normal(0, 0.1, 2)
# Calculate acceptance probability
accept_prob = min(1, posterior(theta_proposal, x, y)/posterior(theta_current, x, y))
# Accept or reject the proposal
if np.random.uniform(0, 1) < accept_prob:
theta_current = theta_proposal
# Add current theta to samples
samples.append(theta_current)
return samples
# Run MCMC algorithm
theta_0 = [0, 0]
num_samples = 10000
samples = mcmc(posterior, x, y, theta_0, num_samples)
# Plot results
plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.plot(x, y, 'o', label='Data')
for i in range(0, num_samples, 500):
plt.plot(x, samples[i][0] + samples[i][1]*x,
'--', color='gray', alpha=0.2)
plt.plot(x, theta_true[0] + theta_true[1]*x, '-', color='k', label='True model')
plt.legend()
plt.subplot(122)
plt.hist2d([s[0] for s in samples], [s[1] for s in samples],
bins=50, cmap='Blues')
plt.colorbar()
plt.plot(theta_true[0], theta_true[1], 'ro', markersize=10)
plt.xlabel(r'$\theta_0$')
plt.ylabel(r'$\theta_1$')
plt.tight_layout()
plt.show()
```
在此示例中,我们生成了一些具有噪声的线性数据,并使用MCMC算法拟合了数据。我们使用了一个简单的高斯先验分布,并且没有使用任何高级技术(例如,自适应步长)。然而,此示例展示了MCMC的基本思想和实现方式。