Python关于mcmc

原创

mob649e8168b406 2024-09-13 06:19:05 ©著作权

文章标签 Python 马尔科夫链 sed 文章分类 Python 后端开发

©著作权归作者所有：来自51CTO博客作者mob649e8168b406的原创作品，请联系作者获取转载授权，否则将追究法律责任

MCMC（马尔科夫链蒙特卡洛方法）概述与Python代码示例

引言

马尔科夫链蒙特卡洛（MCMC）是一种用于从概率分布中抽样的统计方法，广泛应用于科学研究、机器学习及数据分析等多个领域。MCMC方法通过构建马尔科夫链，并利用其稳态分布与目标分布相同的特性，实现有效的采样。

在本文中，我们将探讨MCMC的基本原理及其Python实现，并通过实例来展示代码的使用。

MCMC的基本原理

MCMC的核心思想是利用马尔科夫链的性质，生成一系列随机样本，使其分布接近目标分布。简单来说，MCMC算法可以分为以下几个步骤：

初始化样本。
通过某种机制（如接受-拒绝法）进行样本更新。
重复步骤2直至收敛。

通过这种方式，经过多次迭代后，生成的样本将近似于目标分布。

MCMC的流程图

下面是MCMC流程的示意图：

flowchart TD
    A[初始化样本] --> B{样本更新}
    B --> C[接受样本]
    B --> D[拒绝样本]
    C -->|迭代| B
    D -->|迭代| B
    C --> E[收敛]

Python实现

接下来，我们将通过Python库numpy和matplotlib来演示MCMC的基本实现。我们将以一个简单的一维正态分布为例，进行样本的生成。

import numpy as np
import matplotlib.pyplot as plt

def target_distribution(x):
    """目标分布：标准正态分布"""
    return (1/np.sqrt(2 * np.pi)) * np.exp(-0.5 * x ** 2)

def mcmc(num_samples, initial_sample, step_size):
    """MCMC采样"""
    samples = [initial_sample]
    current_sample = initial_sample
    
    for _ in range(num_samples):
        proposed_sample = current_sample + np.random.normal(0, step_size)
        # 计算接受概率
        acceptance_ratio = target_distribution(proposed_sample) / target_distribution(current_sample)
        
        if np.random.rand() < acceptance_ratio:
            current_sample = proposed_sample
            
        samples.append(current_sample)
    
    return samples

# 主程序
initial_sample = 0
num_samples = 10000
step_size = 1.0

# 采样
samples = mcmc(num_samples, initial_sample, step_size)

# 绘制结果
plt.figure(figsize=(12, 6))
plt.hist(samples, bins=30, density=True, alpha=0.6, color='g')
x = np.linspace(-4, 4, 100)
plt.plot(x, target_distribution(x), 'k', linewidth=2)
plt.title('MCMC Sampling from Target Distribution')
plt.xlabel('Sample Value')
plt.ylabel('Density')
plt.show()

在上面的代码中，我们定义了一个标准正态分布作为目标分布，并通过mcmc函数进行MCMC采样。结果通过直方图和目标分布的密度函数图示。