emcee——Quickstart

emcee是一个强大的、易上手的python第三方的库，使用基于蒙特卡洛的算法实现对贝叶斯参数的估计（Bayesian Parameter Estimation）。

首先需要bear in mind的是：采样算法仅仅完成数据的采样（比如使用蒙特卡洛算法实现圆的面积的求解，进而近似计算 $\pi$ 的值），通过采样数据进一步得到相关信息。

如何对一个多维高斯密度函数进行采样

$$p(\vec{x})\propto \exp\left [ -\frac12(\vec x-\vec \mu)^T\Sigma^{-1}(\vec x-\vec \mu)\right]$$

首先导入必要的库：

import numpy as np
import emcee

编码实现概率密度（关于 $\vec \mu,\Sigma$）的计算，需首先转换为对数形式：

$$\ln p(\vec x)=-\frac12(\vec x-\vec \mu)^T\Sigma^{-1}(\vec x-\vec \mu)$$

def lnpdf(x, mu, icov):
    diff = x - mu
    return -np.dot(diff, np.dot(icov, diff))/2
        # numpy中的二维数组与一维数组进行矩阵乘法时的灵活处理

构造协方差矩阵（$\Sigma$），以50维为例：

ndim = 50
mu = np.random.rand(ndim)

# 首先构造对称矩阵
cov = np.random.rand(ndim**2).reshape((ndim, ndim))-.5
cov = np.triu(cov)
cov += cov.T - np.diag(cov.diagonal())
cov = np.dot(cov, cov)

# 协方差矩阵的逆
icov = np.linalg.inv(cov)

定义emcee库最为核心的类EnsembleSampler构造函数所需的参数：

ndim = 50
mu = np.random.rand(ndim)
nwalkers = 250
nsteps = 1000
p0 = np.random.rand(nwalkers*ndim).reshape((nwalkers, ndim))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnpdf, args=(mu, icov))
sampler.run_mcmc(p0, nsteps)

注意：

samples = sampler.chain  
        # samples表示采样数据
        # shape = (nwakers, nsteps, ndim)
samples = sampler.chain.reshape((-1, ndim))
        # 等价于 sampler.flatchain

显示每一维的数据：

plt.figure(facecolor='w', edgecolor='k')
for i in range(ndim):
    plt.subplot(nim//5, 5, i+1)
    plt.xticks([])
    plt.yticks([])
    plt.hist(sampler.flatchain[:, i], 100, histtype='step')
plt.show()

[1] emcee Quickstart