大数跨境
0
0

MCMC(Markov Chain Monte Carlo)的理解与实践(Python)

MCMC(Markov Chain Monte Carlo)的理解与实践(Python) 数据皮皮侠
2020-12-05
1
导读:微信公众号:数据皮皮侠如果你觉得该公众号对你有帮助,欢迎关注、推广和宣传内容目录:MCMC(Markov C

微信公众号:数据皮皮侠
如果你觉得该公众号对你有帮助,欢迎关注、推广和宣传

内容目录:MCMC(Markov Chain Monte Carlo)的理解与实践(Python)

Markov Chain Monte Carlo (MCMC) methods are a class of 
algorithms for sampling from a probability distribution
based on constructing a Markov chain that has the desired 
distribution as its stationary distribution. The state of 
the chain after a number of steps is then used as a sample 
of the desired distribution. The quality of the sample 
improves as a function of the number of steps.


With MCMC, we draw samples from a (simple) proposal 
distribution so that each draw depends only on the state of 
the previous draw (i.e. the samples form a Markov chain,
 θp=θ+Δθ,Δθ∼N(0,σ)θp=θ+Δθ,Δθ∼N(0,σ)).

MCMC(Markov Chain Monte Carlo)用MCMC采样算法实现对Beta 分布的采样

MCMC(Markov Chain Monte Carlo)

首先来看经典的MCMC采样算法:


用MCMC采样算法实现对Beta 分布的采样

关于Beta distribution更详尽的内容请参见 Beta函数与Gamma函数及其与Beta分布的关系。已知Beta distribution的概率密度函数(pdf)为:

import numpy as np
import scipy.special as ss
import matplotlib.pyplot as plt

def beta_s(x, a, b):
    return x**(a-1)*(1-x)**(b-1)
def beta(x, a, b):
    return beta_s(x, a, b)/ss.beta(a, b)

def plot_mcmc(a, b):
    cur = np.random.rand()
    states = [cur]
    for i in range(10**5):
        next, u = np.random.rand()
        if u < np.min((beta_s(next, a, b)/beta_s(cur, a, b), 1)):
            states.append(next)
            cur = next
    x = np.arange(01.01)
    plt.figure(figsize=(105))
    plt.plot(x, beta(x, a, b), lw=2, label='real dist: a={}, b={}'.format(a, b))
    plt.hist(states[-1000:], 25, normed=True, label='simu mcmc: a={}, b={}'.format(a, b))
    plt.show()

if __name__ == '__main__':
    plot_mcmc(0.10.1)
    plot_mcmc(11)
    plot_mcmc(23)



【声明】内容源于网络
0
0
数据皮皮侠
社科数据综合服务中心,立志服务百千万社科学者
内容 2137
粉丝 0
数据皮皮侠 社科数据综合服务中心,立志服务百千万社科学者
总阅读16
粉丝0
内容2.1k