Metropolis-Hasting算法

当需要进行采样的分布 $p(x)$比较复杂时，一种方法是借助与该分布对应的马尔科夫链状态转移函数来进行采样。那么要怎么找到该分布的状态转移函数呢？如果随便假定要给转移函数 $Q(x_j|x_i)$显然是不成立的，因为不满足细致平稳条件，即

$$p(x_j)*Q(x_i|x_j)\neq p(x_i)*Q(x_j|x_i)$$

但我们可以对上面的式子做一个改造，使细致平稳条件成立。对于任意一个分布 $p(x)$，以及任意一个转移函数 $Q(x_j|x_i)$，我们可以构建一个新的转移函数

$$M(x_j|x_i)=Q(x_j|x_i)*min(1,\frac{p(x_j)*Q(x_i|x_j)}{p(x_i)*Q(x_j|x_i)})$$

使得细致平稳条件一定满足，即

$$p(x_j)*M(x_i|x_j)= p(x_i)*M(x_j|x_i)$$

在进行采样时，先用任意构建的 $Q(x_j|x_i)$进行采样，接着对接受概率 $a(x_j|x_i)=min(1,\frac{p(x_j)*Q(x_i|x_j)}{p(x_i)*Q(x_j|x_i)})$进行分类讨论：

当$p(x_j)*Q(x_i|x_j)≥p(x_i)*Q(x_j|x_i)$，有 $M(x_j|x_i) =Q(x_j|x_i)$

即直接接受 $Q(x_j|x_i)$的采样；

当$p(x_j)*Q(x_i|x_j)<p(x_i)*Q(x_j|x_i)$，有 $M(x_j|x_i) =Q(x_j|x_i)*\frac{p(x_j)*Q(x_i|x_j)}{p(x_i)*Q(x_j|x_i)}$

此时，再以 $\frac{p(x_j)*Q(x_i|x_j)}{p(x_i)*Q(x_j|x_i)}$的概率来决定是否接受采样。

为什么新的转移函数能保证符合细致平稳条件呢？以下是证明过程：

将 $M(x_j|x_i)$的构建定义代入

$$p(x_i)*M(x_j|x_i)=p(x_i)*Q(x_j|x_i)*min(1,\frac{p(x_j)*Q(x_i|x_j)}{p(x_i)*Q(x_j|x_i)})\\=min(p(x_i)*Q(xj|x_i),p(x_j)*Q(x_i|x_j))\\=p(x_j)*Q(x_i|x_j)*min(\frac{p(x_i)*Q(x_j|x_i)}{p(x_j)*Q(x_i|x_j)},1)\\= p(x_j)*M(x_i|x_j)$$

这样我们就证明了新转换函数 $M(x_j|x_i)$是必定满足细致平稳条件的。需要注意的是，这里的 $Q(x_j|x_i)$可以是任意构建的转换函数（e.g.，高斯函数），上面式子恒成立。

• Metropolis-Hasting 采样——转换函数为高斯函数实例

接下来，我们以一个高斯函数的实例来帮助更好理解。

假设目标采样分布为 $p(x)=\frac{0.3 * e^{-(x-0.3)^2}+0.7*e^{-\frac{(x-2.)^2}{0.3}}}{1.2113}$

def p(x):
    return ((0.3 * np.exp(-(x-0.3)**2))+0.7*np.exp(-(x-2.)**2/0.3))/1.2113

要对目标分布进行采样，我们可以按照上面提到的流程进行:

1. 构建一个高斯函数作为提议分布，如 $Q(x_j|x_i) =e^{-(x_j-x_i)^2}$

2. 从提议分布中采样出 $x_j$

3. 计算出接受概率 $a(x_j|x_i)=min(1,\frac{p(x_j)*Q(x_i|x_j)}{p(x_i)*Q(x_j|x_i)})$，由于高斯函数是对称的， $Q(x_j|x_i) =Q(x_i|x_j)$，分式上下 $Q$ 被约去，仅需要计算 $min(1,\frac{p(x_j)}{p(x_i)})$

4. 根据接受概率和随机数，决定是否接受本次采样

nSample = 100000 #进行采样的次数
sample2 = [] #采样的集合

x = 1 # 设置采样的初始值（x_0），初始状态任意

for i in range(nSample):
    x2 = norm(loc=x).rvs() # 从提议分布（高斯）采样一个
    a = np.min([1, p(x2)/p(x)]) #接受概率，因为提议分布高斯是对称的, 所以可以上下约掉
    u = np.random.rand() # 从[0,1]均匀分布获取一个值
    if u <= a:
        sample2.append(2) # 接受并跳转到新的采样点
        x = x2
    else:
        sample2.append(x) # 保留原来的采样点不变

完成采样以后，我们可以画出分布

x = np.arange(-3., 5., 0.01)
plt.plot(x, p(x), color='r', label='p(x)')
a,_,_=plt.hist(sample2, color='k', bins=150, density=True)
plt.legend()
plt.ylabel('pdf(x)')

Metropolis-Hasting 算法在参数估计的应用

在认知心理学中，我们经常需要根据实验数据对概率模型进行参数估计，本质就是求解参数的后验概率分布 $p(\theta|data)$。当我们使用Metropolis-Hasting 算法对参数的后验分布进行采样时，只需要将上面的 $p(x)$ 替换为 $p(x|data)$，如此在计算接受概率时，

$$a(x_j|x_i)=min(1,\frac{p(x_j|data)*Q(x_i|x_j)}{p(x_i|data)*Q(x_j|x_i)})$$

由于后验分布 $p(x|data)$我们是不知道的，我们需要用贝叶斯公式将其转为可以计算的形式

$$a(x_j|x_i)=min(1,\frac{p(data|x_j)*p(x_j)/p(data)*Q(x_i|x_j)}{p(data|x_i)*p(x_i)/p(data)*Q(x_j|x_i)})\\=min(1,\frac{p(data|x_j)*p(x_j)*Q(x_i|x_j)}{p(data|x_i)*p(x_i)*Q(x_j|x_i)})$$

• 如果我们假定的是均匀的先验分布，那么有 $p(x_i)=p(x_j)$；

• 如果我们假定的是对称的提议分布（比如高斯分布），那么有 $Q(x_j|x_i) =Q(x_i|x_j)$

那么上式可以进一步简化为

$$a(x_j|x_i)=min(1,\frac{p(data|x_j)}{p(data|x_i)})$$

可以看到，接受概率的右边就是观测数据在新采样参数与原采样参数中的似然值之比。接下来我们具体介绍一个用MH算法做参数估计的代码实例。

• 用MH算法求解心理物理曲线的参数

Weibull函数为

$$y=1-(1-\gamma)exp^{-(\frac{k*x}{\alpha})^\beta}$$

其中

$$k=-log(\frac{1-g}{1-\gamma})^\frac{1}{\beta}$$

$\alpha$为该正确率下的阈值，$\beta$为函数的斜率，$\gamma$为基线概率水平的概率值， $g$为阈值对应的正确率。

# 先定义weibull函数
def weibullfun(x, alpha, beta, gamma, g):
    k = (-np.log((1-g)/(1-gamma)))**(1/beta)
    y = 1 - (1 - gamma)*np.exp(-(k*x/alpha) ** beta) # hack here to avoid the complex number of 
    return y

假设我们当前的参数组合 $(\alpha,\beta,\gamma,g)$的真值为（0.25，3，0.5，0.82），在实验设置的一致性水平（[0.032, 0.064, 0.128, 0.256, 0.512]）中对应的心理物理曲线如下，

from scipy.stats import bernoulli
# coherence level in literature
coh = np.array([0.032, 0.064, 0.128, 0.256, 0.512]) 
# 定义参数
g = 0.82
gamma = 0.5 # Baseline chance level
beta = 3. # slop steepness
alpha = 0.25 # coherence threshold corresponding to g

nTrialPerCoh = 1000 # how many trials per coh level

pCoh = weibullfun(coh, alpha, beta, gamma, g)

plt.plot(coh, pCoh, 'o-')
plt.xlabel('coherence')
plt.ylabel('Prob of correct')

我们通过伯努利随机过程根据以上心理曲线生成5000个模拟数据，五个一致性水平分别生成1000个试次，

data = np.empty((coh.size, nTrialPerCoh))
for iCoh in range(coh.size): # loop coherence level
    # use bernouli process to generate data
    data[iCoh, :] = bernoulli.rvs(pCoh[iCoh], size=nTrialPerCoh)
cohTrial = np.tile(coh[:, np.newaxis], (1, nTrialPerCoh)).flatten() #(5000, ) array
data = data.flatten() # (5000, ) array

接着，我们假定Weibull函数中的阈限 $\alpha$是需要估计的自由参数。我们使用MH算法对后验分布 $p(\alpha|data)$进行采样，选择对称的高斯分布进行提议分布，便可以使用公式 $a(x_j|x_i)=min(1,\frac{p(data|x_j)}{p(data|x_i)})$计算接受概率。此时，我们考虑直接计算似然的比值有些麻烦，可以两边取对数，改为计算两次采样参数的对数似然的差值

$$log(a(x_j|x_i))=min(0,log(p(data|x_j))-log(p(data|x_i)))$$

写出不同 $\alpha$输入相应的负对数似然函数

eps = np.finfo(float).eps #大于0的极小值，防止后续计算对数时出现0的情况
def loglikeli(alpha, coh, data):
    gamma = 0.5
    g = 0.82
    beta = 3.

    prob = np.empty(coh.size)
    for i in range(coh.size):
        p = weibullfun(coh[i], alpha, beta, gamma, g)
        prob[i] = p*data[i]+(1-p)*(1-data[i])
    
    prob = prob*0.999+eps # 避免prob里面有0的情况, log(0)=-infinity
    return -np.log(prob).sum() 
    #prob取值在（0，1]，经过对数函数后会变为负数，取负的对数似然值能保证返回值是一个正数

MCMC采样从初始值0（可任意选取）开始，共采样2000次，剔除前1000次未平稳收敛的采样，保留后1000个有效样本。

from scipy.stats import norm

nSample = 1000
nBurnin = 1000
nTotalSample = nSample + nBurnin
stepsize = 1

samples = []
oldsample = 0

由于提议分布是在实数范围内进行采样，但需要估计的阈限 $\alpha$参数范围在（0,1），我们需要进行一个尺度转换，通过logistic函数可以将实数范围投射到（0，1）之中

trans_oldsample = 1/(1+np.exp(-oldsample))

每次采样时，会在均值为上一采样点的高斯分布中采样，进行尺度转换后通过负对数似然函数计算接受概率，并根据随机数决定是否接受本次采样。

for i in range(nTotalSample):
    newsample = norm.rvs(loc=oldsample,scale = stepsize)
    # 我们的参数范围在(0,1),但是我们转化到纯实数范围内采集
    trans_newsample = 1/(1+np.exp(-newsample))

    # 计算负对数似然之差
    delta = negloglikeli(trans_oldsample, cohTrial, data)-negloglikeli(trans_newsample, cohTrial, data)
    # 接受概率
    a = np.min([0,delta]) 
    #从[0,1]均匀分布随机获取一个值
    u = np.random.rand()
    
    if log(u) <= a: #接受并跳转到新的采样点
        samples.append(trans_newsample)
        trans_oldsample = trans_newsample
        oldsample = newsample
    else: #保留原来的采样点不变
        samples.append(trans_oldsample)

完成之后可以从样本计算得到所估计参数阈限 $\alpha$的均值，并画出有效采样的直方图。由于模拟数据与采样具有随机性，运行的结果可能有所不同。

print(np.mean(samples[1000:]))
plt.hist(samples[1000:], bins=20,density=True,range=[0.2,0.3])

0.24384956185678375

可以看到使用MH算法估计的参数 $\alpha$的均值，接近参数设置的真值0.25。