Markov Chain Monte Carlo方法的一个解释 贝叶斯推理

main.py

#
# 代码来自 http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/
#
#

import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from scipy.stats import norm


def sampler(data, samples=4, mu_init=.5, proposal_width=.5, plot=False, mu_prior_mu=0, mu_prior_sd=1.):
    '''

    mcmc通过在后验分布中进行采样，得到一个大样本，然后使用直方图的方法来得到后验分布的形式。

    在此处，后验分布是关于mu的一个分布。
    最开始随机选择一个值，认为这个值是从后验分布采样得来。

    以这个值为标准，使用一个方法得到一个新的值，如果这个新的值的后验概率比旧的值的后验概率大于一个随机数，
    那么将这个新的值认为是一个新的采样，同时以这个新的值作为标准，重复上面的过程。

    可以有很多方法来从当前值获得一个新的值（markov chain阶段），metropolis方法是其中一种方法， 它在一个以当前值为mu，以
    proposal_width为std的正态分布中进行采样得到新的值。其中proposal_width是一个超参，任意指定proposal_width都会收敛到
    同样的结果，但是收敛的速度不一样。

    如果这个新的值的后验概率比旧的值的后验概率大于一个随机数，接受新值作为一个采样，否则旧值是一个采样（monte carlo阶段）。
    计算后验概率的比值，可以去除p(x)的计算，只需要计算似然函数和先验。
    我们会接受一些新后验比当前后验低的样本， 目的是为了得到后验的概率分布。

    Parameters:
    ------------------
    data: list
        样本点
    samples: int
        markov chain monte carlo 采样次数
    mu_init: float
        第一个采样点
    proposal_width: float
        metropolis sampler要求的参数， metropolis在一个以当前位置为mu，proposal_width为std的正态分布中进行采样得到下一个采样点
    plot: boolean
        是否绘图
    mu_prior_mu: float
        关于参数的先验分布， 这里假设参数的先验是正态分布
    mu_prior_std: float
        关于参数的先验分布， 这里假设参数的先验是正态分布

    Reference:
    ---------------------
    http://twiecki.github.io/blog/2015/11/10/mcmc-sampling/

    '''
    mu_current = mu_init
    posterior = [mu_current]
    for i in range(samples):
        # suggest new position
        mu_proposal = norm(mu_current, proposal_width).rvs()

        # Compute likelihood by multiplying probabilities of each data point
        likelihood_current = norm(mu_current, 1).pdf(data).prod()
        likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()
        
        # Compute prior probability of current and proposed mu        
        prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)
        prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)
        
        p_current = likelihood_current * prior_current
        p_proposal = likelihood_proposal * prior_proposal
        
        # Accept proposal?
        p_accept = p_proposal / p_current
        
        # Usually would include prior probability, which we neglect here for simplicity
        accept = np.random.rand() < p_accept
        
        if plot:
            plot_proposal(mu_current, mu_proposal, mu_prior_mu, mu_prior_sd, data, accept, posterior, i)
            pass
        
        if accept:
            # Update position
            mu_current = mu_proposal
        
        posterior.append(mu_current)
        
    return posterior

def calc_posterior_analytical(data, x, mu_0, sigma_0):
    '''
    如果先验是正态分布，
    似然函数也是正态分布（数据服从正态分布），
    那么后验也服从正态分布。
    这种后验和先验的分布形式一致的现象，叫做共轭。
    这样一种后验，可以使用数学分析的方法得到它的封闭解。
    

    '''
    sigma = 1.
    n = len(data)
    mu_post = (mu_0 / sigma_0**2 + data.sum() / sigma**2) / (1. / sigma_0**2 + n / sigma**2)
    sigma_post = (1. / sigma_0**2 + n / sigma**2)**-1
    return norm(mu_post, np.sqrt(sigma_post)).pdf(x)

# Function to display
def plot_proposal(mu_current, mu_proposal, mu_prior_mu, mu_prior_sd, data, accepted, trace, i):
    from copy import copy
    trace = copy(trace)
    fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, figsize=(16, 4))
    fig.suptitle('Iteration %i' % (i + 1))
    x = np.linspace(-3, 3, 5000)
    color = 'g' if accepted else 'r'
        
    # Plot prior
    prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)
    prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)
    prior = norm(mu_prior_mu, mu_prior_sd).pdf(x)
    ax1.plot(x, prior)
    ax1.plot([mu_current] * 2, [0, prior_current], marker='o', color='b')
    ax1.plot([mu_proposal] * 2, [0, prior_proposal], marker='o', color=color)
    ax1.annotate("", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),
                 arrowprops=dict(arrowstyle="->", lw=2.))
    ax1.set(ylabel='Probability Density', title='current: prior(mu=%.2f) = %.2f\nproposal: prior(mu=%.2f) = %.2f' % (mu_current, prior_current, mu_proposal, prior_proposal))
    
    # Likelihood
    likelihood_current = norm(mu_current, 1).pdf(data).prod()
    likelihood_proposal = norm(mu_proposal, 1).pdf(data).prod()
    y = norm(loc=mu_proposal, scale=1).pdf(x)
    sns.distplot(data, kde=False, norm_hist=True, ax=ax2)
    ax2.plot(x, y, color=color)
    ax2.axvline(mu_current, color='b', linestyle='--', label='mu_current')
    ax2.axvline(mu_proposal, color=color, linestyle='--', label='mu_proposal')
    #ax2.title('Proposal {}'.format('accepted' if accepted else 'rejected'))
    ax2.annotate("", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),
                 arrowprops=dict(arrowstyle="->", lw=2.))
    ax2.set(title='likelihood(mu=%.2f) = %.2f\nlikelihood(mu=%.2f) = %.2f' % (mu_current, 1e14*likelihood_current, mu_proposal, 1e14*likelihood_proposal))
    
    # Posterior
    posterior_analytical = calc_posterior_analytical(data, x, mu_prior_mu, mu_prior_sd)
    ax3.plot(x, posterior_analytical)
    posterior_current = calc_posterior_analytical(data, mu_current, mu_prior_mu, mu_prior_sd)
    posterior_proposal = calc_posterior_analytical(data, mu_proposal, mu_prior_mu, mu_prior_sd)
    ax3.plot([mu_current] * 2, [0, posterior_current], marker='o', color='b')
    ax3.plot([mu_proposal] * 2, [0, posterior_proposal], marker='o', color=color)
    ax3.annotate("", xy=(mu_proposal, 0.2), xytext=(mu_current, 0.2),
                 arrowprops=dict(arrowstyle="->", lw=2.))
    #x3.set(title=r'prior x likelihood $\propto$ posterior')
    ax3.set(title='posterior(mu=%.2f) = %.5f\nposterior(mu=%.2f) = %.5f' % (mu_current, posterior_current, mu_proposal, posterior_proposal))
    
    if accepted:
        trace.append(mu_proposal)
    else:
        trace.append(mu_current)
    ax4.plot(trace)
    ax4.set(xlabel='iteration', ylabel='mu', title='trace')
    plt.tight_layout()
    #plt.legend()
    plt.show()


if __name__ == '__main__':
    data = np.random.randn(20)
    sampler(data, plot=True)