码迷,mamicode.com
首页 > 其他好文 > 详细

MCMC例子

时间:2020-06-29 00:47:14      阅读:124      评论:0      收藏:0      [点我收藏+]

标签:ipy   figure   mic   随机   append   sci   矩阵   alpha   sub   

技术图片
‘‘‘‘

假设目标分布是: 一维的正太分布, i.e., N(0, 1)
建议的转移矩阵Q(i,j)也是正太分布, j 服从 N(i, 2^2)

pi(i)Q(i,j)*alpha(i,j) = pi(j)Q(i,j)*alpha(i,j)
where alpha(i,j) = pi(j)Q(j,i) 接受率
alpha(i,j) 可能很小,导致接受率比较小,采样效率低下
修正: alpha(i,j) =  min(1, alpha(i,j ) / alpha(j, i))
                 =  min(1, pi(j)Q(j,i) / pi(i)Q(i,j))
                 =  min(1, pi(j) / pi(i))   # 通常Q(i,j)是对称的, 并且此时分布pi都无需归一化


Gibbs采样:
    考虑二维的情况,  A = (x1, y1), B = (x1, y2), C = (x2, y2)
    那么有: pi(A)pi(B|A) = p(x1)p(y1|x1) * p(y2|x1)
            pi(B)pi(A|B) = p(x1)p(y2|x1) * p(y1|x1)
    上面两个式子相等:
        ==>  pi(A)*p(y2|x1) = pi(B)*p(y1|x1),  A和B具有相同的x坐标x1
    同理,我们有
            pi(B)*p(x2|y2) = pi(C)*p(x1|y2),  B和C具有相同的y坐标y2
    因此,我们可以构造出接受率为1的状态转移矩阵,如上定义,沿着某个坐标轴进行采样,其余坐标保持不变
    高维的情况类似
‘‘‘

from scipy import stats
import numpy as np
import matplotlib.pyplot as plt

def MCMC_01():
    def target_pdf(x):
        return stats.norm.pdf(x, loc=0, scale=1)  # 随机变量的概率密度函数

    X_0 = 2.0
    Sample_X = [X_0]
    for t in range(100000):
        X_i = Sample_X[-1]
        X_j = np.random.randn() * 2 + X_i  # 服从分布 N(i, 2)
        alpha = min(1, target_pdf(x=X_j) / target_pdf(x=X_i))
        p = np.random.rand()
        if p < alpha: #接受
            Sample_X.append(X_j)
        else:
            Sample_X.append(X_i)
    show_cnt = 10000
    Sample_X = Sample_X[-show_cnt:]
    Standard_X = np.random.randn(show_cnt)

    plt.figure(figsize=(15, 9))
    plt.subplot(1, 2, 1)
    plt.hist(Sample_X, bins=50, normed=1)
    plt.xlabel(x)
    plt.xlabel(p)
    plt.title(MCMC)
    plt.subplot(1, 2, 2)
    plt.hist(Standard_X, bins=50, normed=1)
    plt.xlabel(x)
    plt.title(Normal(0, 1))
    plt.show()

MCMC_01()



def MCMC_02():
    ‘‘‘
        目标分布是二维的正太分布:N (u, sigma^2),
                u= (u1, u2) = (3, 5),
                sigma = [[sigma_1^2,          rho*sigma_1*sigma_2],
                        [rho*sigma_1*sigma_2, sigma_2^2]]
                      = [[1, 2],
                        [2, 4]]
        那么对应的条件分布为:
            P(x1|x2) = N(u1 + (x2 - u2) * rho * sigma_1 / sigma_2, (1-rho^2) * sigma_1^2)
            P(x2|x1) = N(u2 + (x1 - u1) * rho * sigma_2 / sigma_1, (1-rho^2) * sigma_2^2)
    ‘‘‘
    u1, u2 = 3, 5
    sigma_1, sigma_2, rho = 1, 2, 0.5
    def gibbs_sampling_x2_given_x1(x1):
        # return: p(x2|x1)
        now_u = u2 + (x1 - u1) * rho * sigma_2 / sigma_1
        now_sigma = np.sqrt(1 - rho ** 2) * sigma_2 ** 2
        now_sigma = np.sqrt(now_sigma)
        return np.random.randn() * now_sigma + now_u
    def gibbs_sampling_x1_given_x2(x2):
        # return: p(x1|x2)
        now_u = u1 + (x2 - u2) * rho * sigma_1 / sigma_2
        now_sigma = (1-rho ** 2) * sigma_1 ** 2
        now_sigma = np.sqrt(now_sigma)
        return np.random.randn() * now_sigma + now_u


    X_0 = [-2.0, -3.0]
    Sample_X = [X_0]
    for t in range(100000):
        x_2_pre = Sample_X[-1][-1]
        x_1 = gibbs_sampling_x1_given_x2(x2=x_2_pre)  #根据上一个采样点的x2, 采样现在的x1
        x_2 = gibbs_sampling_x2_given_x1(x1=x_1) #根据刚刚采样点的x1, 采样现在的x2
        Sample_X.append([x_1, x_2])
        if t<10:
            print(x_1:{}  x_2:{}.format(x_1, x_2))

    show_cnt = 10000
    Sample_X = np.array(Sample_X[-show_cnt:])
    # Standard_X = np.random.randn(show_cnt)
    print(Sample_X:, Sample_X.shape)
    plt.figure(figsize=(10, 9))
    plt.subplot(1, 1, 1)
    plt.plot(Sample_X[:, 0], Sample_X[:, 1], r.)
    plt.xlabel($x_1$)
    plt.ylabel($x_2$)
    plt.title(Gibbs)
    plt.show()

MCMC_02()
View Code

技术图片

 

技术图片

 

MCMC例子

标签:ipy   figure   mic   随机   append   sci   矩阵   alpha   sub   

原文地址:https://www.cnblogs.com/skykill/p/13205692.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!