5.1.5 MCMC在参数估计中的应用

在第四章的心理物理曲线部分，我们利用最大似然估计的方法拟合了心理物理曲线。这种方法的问题在于，如果中间的隐变量无法通过积分或求和消除，就很难求解。

我们通过手写MH算法来估计心理物理曲线。

实验任务是一个随机点运动任务。实验设置的一致性水平为[0.032, 0.064, 0.128, 0.256, 0.512]。已知该被试的行为反应数据，包含$5个一致性水平 \times 1000个试次 = 5000 个试次$ 。我们的目标是用Weilbull函数拟合一位被试的心理物理曲线，并估计出82%阈值。

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$为阈值对应的正确率。

该被试的心理物理曲线参数组合 $(\alpha,\beta,\gamma,g)$的真值为（0.25，3，0.5，0.82），其产生的行为数据如下。具体的数据模拟代码见第四章的心理物理曲线部分

首先，与最大似然估计类似，我们先写出关于 $\alpha$的负对数似然函数：

# 定义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

# 定义负对数似然函数
eps = np.finfo(float).eps #大于0的极小值，防止后续计算对数时出现0的情况
def negloglikeli(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]，经过对数函数后会变为负数，取负的对数似然值能保证返回值是一个正数

使用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)])$$

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。

在采样数量较小时，分布可能很不稳定，可以通过提高采样数量来获得更高精度，但是运行速度会更加慢。

MCMC工具包选择指南

MCMC的算法本身和具体的模型无关，只是一种采样的方法，所以现在有许多工具包可以实现MCMC的算法。

目前的主流MCMC算法包括:

• Metropolis-Hasting sampler

• Gibbs sampler

• Hamiltonian sampler

• No-U-turn sampler

目前的主流Bayesian编程工具包包括：

• Stan (基于C++, 有pystan, rstan接口，推荐)

• PYMC5 (基于python，不方便写for循环)

• Numpyro (基于python，不方便写for循环)

• Jags, winbugs (已被淘汰，不推荐)