一种基于AEPD的贝叶斯参数估计方法

A Technique of Bayesian Parameter Estimation Based on AEPD

摘要: 本文基于非对称指数幂分布(AEPD)，采用MCMC中的M-H算法进行贝叶斯参数估计，针对后验密度的复杂性，采取随机游走链的建议分布为中心对称分布。通过数值模拟，在M个大样本下进行N次迭代，模拟结果显示，该算法对具有非对称指数幂分布的贝叶斯参数估计具有实用性。

Abstract: In this paper, based on the asymmetric exponential power distribution (AEPD), the M-H algorithm in MCMC is used to estimate the Bayesian parameters. For the complexity of the posterior density, the recommended distribution of the random walk chain is the center symmetric distribution. Through numerical simulation, N iterations are performed under M large samples. The simulation results show that the algorithm is practical for Bayesian parameter estimation with asymmetric exponential power distribution.

1. 引言

随着社会的发展，数据的复杂性越来越大。在经济领域，已经很难找到服从正态分布的数据集，基于正态分布所建立的一些模型方法也不再适用于对复杂数据的分析。2009年由Zhu Dongming [1] 等提出的非对称指数幂分布(Asymmetric Exponential Power Distribution, AEPD)对数据的适应性更强，该分布是一个新的分布族，比如正态分布、均匀分布都是AEPD的特例。

基于此，针对AEPD的贝叶斯参数估计采用MCMC中的M-H算法作进一步探讨。首先，1995年Chib [2] 详细介绍了Metropolis-Hastings算法理论，并对该算法在应用上做了讨论，即使在特殊情况下，该算法也能体现出它的优良性；其次，2002年Aans Altaleb [3] 使用MCMC中的M-H算法在Logit回归模型进行参数估计，以及2003年John Geweke [4] 采用独立链、随机游走链和泰勒链对贝叶斯状态空间模型进行估计；最后，2015年David D.L. Minh [5] 对M-H算法进行两种互补和直观的方式来进行推导和理解该算法，并得出该算法具有较高的接受概率。

本文以非对称指数幂分布为基础，利用MCMC中的M-H算法思想，对此进行贝叶斯参数估计。

2. Metropolis-Hastings算法

在统计领域中主要有两大计算方法，一种是极大似然估计计算，另一种是Bayes计算。在Bayes计算中，最重要的算法就是Markov chain Monte Carlo (MCMC)中的Metropolis-Hastings和Gibbs算法。在本文中主要研究Metropolis-Hastings在贝叶斯参数估计的应用。

在统计学中，对连续变量求数学期望，公式如下：

(1)

其中，X是连续随机变量， $p\left(x\right)$ 为随机变量X的密度函数。

假设我们要计算某分布函数 $F\left(x\right)$ 关于概率密度函数 $\pi \left(x\right)$ 的数学期望，根据公式(1)可得：

$E\left(F\right)={\displaystyle {\int }_{X}F\left(x\right)\pi \left(x\right)\text{d}x}$ (2)

其中 $F\left(x\right)$ 和 $\pi \left(x\right)$ 较为简单，可以由上式求出。但是，在实践中，由分布函数或概率密度函数推出的后验分布或后验密度函数都是复杂、高维且非标准形式的，对此就要进行MCMC计算。

在实际中，假设 $\pi \left(x\right)$ 是未知的，我们可以计算归一化常数范围内的密度来替换 $\pi \left(x\right)$，公式如下：

$\begin{array}{l}\pi \left(x\right)=Cg\left(x\right)\\ C=1/{\displaystyle {\int }_{X}g\left(x\right)\text{d}x}\end{array}$ (3)

其中， $g\left(x\right)$ 是先验分布乘以似然函数，C是归一化常数。

在贝叶斯统计中，Metropolis-Hastings是最常用的算法。假设X的分布依赖我们想要估计的参数 $\vartheta$，对于 $F\left(x\right)\ge 0$，该分布的似然函数为 $L\left(x|\vartheta \right)$，根据(1)、(2)可知，则数学期望为：

$E\left(F\left(X\right)\right)={\displaystyle \int L\left(x|\vartheta \right)F\left(X\right)\text{d}x}$ (4)

在贝叶斯推断中，统计模型的观测值和参数变量都是随机的。给定样本数据 $X=\left(x_1,x_2,\cdots ,x_n\right)$ 、统计模型参数 $\vartheta$ 和先验分布，根据贝叶斯定理，模型的参数联合后验密度为：

$p\left(\vartheta |x\right)=\frac{L\left(x|\vartheta \right)p\left(\vartheta \right)}{{\displaystyle \int L\left(x|s\right)p\left(s\right)\text{d}s}}=C_{X}L\left(x|\vartheta \right)p\left(\vartheta \right)$ (5)

其中， $L\left(x|\vartheta \right)$ 为统计模型的似然函数， $C_X$ 依赖于X的归一化常数，在M-H算法中，可以不用计算 $C_X$，则该参数联合后验密度可表示为：

$p\left(\vartheta |x\right)=L\left(x|\vartheta \right)p\left(\vartheta \right)$ (6)

在M-H算法中，选择一个建议分布函数对(6)进行参数更新，建议分布如下：

$q\left(\vartheta ,x\right)=\pi \left(x\right)$ (7)

$q\left(\vartheta ,y\right)$ 是与目标分布 $p\left(\vartheta |x\right)$ 较接近的已知分布，该建议分布称为备选生成密度，应满足：

${\displaystyle \int q\left(\vartheta ,x\right)\text{d}x}=1$ (8)

在进行M = H算法时，应选择一个接受概率 $\alpha \left(\vartheta ,{\vartheta }^{\prime }\right)$ 使相应的 $q\left(\vartheta ,y\right)$ 以 $p\left(\vartheta |x\right)$ 为其平稳分布，最常用的选择是：

$\alpha \left(\vartheta ,{\vartheta }^{\prime }\right)=\min \left\{\frac{p\left({\vartheta }^{\prime }|x\right)q\left({\vartheta }^{\prime },x\right)}{p\left(\vartheta |x\right)q\left(\vartheta ,x\right)},1\right\}$ (9)

3. 非对称指数幂分布贝叶斯参数估计

在经济领域，很难找到服从正态分布的数据，往往都是尾重或者带有偏度的数据。如果采取假设正态分布的方法进行研究，那么对数据的分析就会稍差一些，很难完全解读数据带给我们的信息。2009年Zhu Dongming [1] 提出了非对称指数幂分布(Asymmetric Exponential Power Distribution, AEPD)，该分布对于数据的适应性很强，其概率密度函数为：

$f_{\text{AEP}}\left(y|\mu ,\sigma ,\alpha ,p_1,p_2\right)=\{\begin{array}{l}\frac{1}{\sigma }\exp \left(-{\left|\frac{y-\mu }{\alpha \sigma /\Gamma \left(1+p_1\right)}\right|}^{p_1}\right),y\le \mu \\ \frac{1}{\sigma }\exp \left(-{\left|\frac{y-\mu }{\left(1-\alpha \right)\sigma /\Gamma \left(1+p_2\right)}\right|}^{p_2}\right),y>\mu \end{array}$ (10)

其中， $\text{Γ}\left(x\right)$ 为gamma函数， $0<\alpha <1$ 为偏度参数， $\sigma >0$ 为尺度参数， $p_1>0$ 和 $p_2>0$ 分别控制左尾和

右尾参数， $-\infty <\mu <+\infty$ 为位置参数。当 $\alpha =\frac{1}{2}$ 和 $p_1=p_2$ 时，AEPD分布为对称的指数幂分布；当 $\alpha =\frac{1}{2}$

和 $p_1=p_2=2$ 时，AEPD分布为标准的正态分布；当 $p_1=p_2=1$ 时便是ALD分布； $p_1\left(p_2\right)$ 值越小，则左尾(右尾)越厚， $\alpha$ 值越大，则分布呈现为左偏。

图1给出了非对称指数幂分布(AEPD)在参数, $\sigma =1$, $p_1=2$, $p_2=2$ 下不同偏度参数 $\alpha$ 的概率密度函数。从图中可以看出，随着 $\alpha$ 的不断增大，该分布呈现左偏现象。当 $\alpha =0.5$ 时，该分布表现为正态分布，可以看到，非对称指数幂分布随着不同参数的变化，可以呈现很多常见的分布，比如正态分布、非对称拉普拉斯分布等。

Figure 1. AEPD probability density function graph as a function of $\alpha$

图1. 随 $\alpha$ 变化的AEPD概率密度函数图

图2给出了非对称指数幂分布在参数 $\mu =0$, $\sigma =1$, $\alpha =0.5$, $p_2=2$ 下随着 $p_1$ 变化的图像。可以看到当 $p_1$ 不断减小时，该分布呈现左尾厚重现象。当 $p_1=2$ 时，该分布呈现正态分布现象。

Figure 2. AEPD probability density function as a function of $p_1$

图2. 随 $p_1$ 变化的AEPD概率密度函数图

图3给出了非对称指数幂分布在参数 $\mu =0$, $\sigma =1$, $\alpha =0.5$, $p_1=2$ 下随着 $p_2$ 变化的图像。可以看到当 $p_2$ 不断减小时，该分布呈现右尾厚重现象。当 $p_2=2$ 时，该分布呈现正态分布现象。

假设总体的密度函数为 $f_{\text{AEP}}\left(y,\vartheta \right)$，其中 $\vartheta =\left(\mu ,\sigma ,\alpha ,p_1,p_2\right)$ 为未知参数向量， $y_i\left(i=1,2,\cdots ,n\right)$ 是来自该总体的样本，由(10)可得，该密度函数的似然函数为：

$L\left(y|\vartheta \right)=\underset{i}{\overset{n}{{\displaystyle \prod }}}\frac{1}{\sigma }\left[\exp \left(-{\left(\frac{\left(\mu -y_i\right)\Gamma \left(1+1/p_1\right)}{\alpha \sigma }\right)}^{p_1}\right)+\exp \left(-{\left(\frac{\left(y_i-\mu \right)\Gamma \left(1+1/p_2\right)}{\left(1-\alpha \right)\sigma }\right)}^{p_2}\right)\right]$ (11)

Figure 3. AEPD probability density function as a function of $p_2$

图3. 随 $p_2$ 变化的AEPD概率密度函数图

对数似然函数为：

$\begin{array}{c}\ln L\left(y|\vartheta \right)=-n\ln \sigma -\underset{i}{\overset{n}{{\displaystyle \sum }}}{\left(\frac{\Gamma \left(1+\frac{1}{p_1}\right)\left(\mu -y_i\right)}{\alpha \sigma }\right)}^{p_1}I\left(y_i\le \mu \right)\\ -\underset{i}{\overset{n}{{\displaystyle \sum }}}{\left(\frac{\Gamma \left(1+\frac{1}{p_2}\right)\left(y_i-\mu \right)}{\left(1-\alpha \right)\sigma }\right)}^{p_2}I\left(y_i>\mu \right)\end{array}$ (12)

在贝叶斯推断中，选择无信息先验分布 $p\left(\vartheta \right)\propto 1$。根据(6)可得联合后验分布为：

$\begin{array}{l}p\left(\vartheta |y\right)\propto L\left(y|\vartheta \right)p\left(\vartheta \right)\\ \propto \frac{1}{{\sigma }^n}\left[\exp \left(-{\left(\frac{{{\displaystyle \sum }}_i^n\left(\mu -y_i\right)\Gamma \left(1+1/p_1\right)}{\alpha \sigma }\right)}^{p_1}\right)+\exp \left(-{\left(\frac{{{\displaystyle \sum }}_i^n\left(y_i-\mu \right)\Gamma \left(1+1/p_2\right)}{\left(1-\alpha \right)\sigma }\right)}^{p_2}\right)\right]\end{array}$ (13)

其中， $y=y_1,\cdots ,y_n$ 为观测值。由于联合后验分布的复杂性，要想对参数进行贝叶斯估计，采用MCMC方法进行模拟。MCMC方法的思想就是构造一个平稳分布对目标分布进行后验抽样。根据(7)、(8)、(9)进行M-H算法，对参数进行贝叶斯估计。

4. 数值模拟

非对称指数幂分布是一个新的分布族，随着参数的变化，呈现出不同的分布，比如正态分布、均匀分布都是该分布的一个特例。对此，利用随机模拟生成1000个服从正态分布的数据，即 $y_i~N\left(0,10\right)$。利用这些模拟数据来进行非对称指数幂分布的贝叶斯参数估计。在进行贝叶斯参数估计的过程中，用到MCMC中的M-H算法，具体步骤如下：

对于M-H算法中的建议分布，针对非对称指数幂分布，采用正态分布作为建议分布，在R语言中将从建议分布中产生的第一个值设为初始化链条的。

a) 选择建议分布 $q\left(\vartheta ,y_i\right)$ 为正态分布的密度，设 $i=\text{2}$，迭代b、c、d步；

b) 从建议分布 $N\left({\vartheta }_{i-1},2\right)$ 为中生成 ${{\vartheta }^{\prime }}_i$ ；

c) 令 ${{\vartheta }^{\prime }}_i$ 的接受概率为公式(9)，其中，选择 $q\left(\vartheta ,y\right)=q\left({\vartheta }^{\prime },y\right)$ 为随机游走的M算法；

d) 产生一个服从均分布 $U\left(0,1\right)$ 的随机数u，令 ${\vartheta }_{i+1}=\{\begin{array}{l}{\vartheta }_i,u\le \alpha \left(\vartheta ,{\vartheta }^{\prime }\right)\\ {\vartheta }_1,u>\alpha \left(\vartheta ,{\vartheta }^{\prime }\right)\end{array}$ ；

e) 令 $i=i+\text{1}$，返回a步。

表1给出了采用MCMC中的随机游走M-H算法，将1000个数据样本进行10,000次迭代，剔除前1000个样本后的后验均值及95%的置信区间，这表明M-H算法在贝叶斯估计中的良好表现，解决了贝叶斯估计中高维数据处理困难的问题。从表中可以看到在不同的分位点每个待估参数不同的模拟估计值，随着分位点的增大，模拟的估计值也在逼近真值，例如在分位点 $\alpha =0.95$ 时，待估参数的真值为： $\mu =0$， $\sigma =1$， $p_1=0.5$， $p_2=2$，相对应的贝叶斯参数估计模拟值为： $\mu =0.002295$， $\sigma =1.002321$， $\alpha =0.948325$， $p_1=0.502165$， $p_2=2.001092$。在表中也可以看到每个待估参数在各个分位点的95%的置信区间。

图4、图5、图6、图7和图8中给出了每个待估参数的模拟轨迹图、参数模拟值直方图、不同分位点轨迹图和自相关图。例如，图4中给出了 $\mu$ 相关的轨迹图和检验收敛图，从左往右依次为：剔除前1000个样本估计值的轨迹图，图中可以看到，均值 $\mu$ 的马尔科夫链收敛；第二个图为参数模拟值的直方图，可以看到抽到的估计值基本都在真值周围；第三个图为在0.025、0.5和0.975分位数下的轨迹图，从图中可以看到在样本1500以后都趋于收敛，其中上下两条虚线分别代表0.025和0.975分位数线的变化，实线代表0.5分位数的变化；第四个图为链的等价自相关函数图，从图中可以看到，最终的自相关度都趋于收敛。

Table 1. Simulated estimates of Bayesian parameters of AEP distributed at different quantiles

表1. AEP分布在不同分位点的贝叶斯参数模拟估计值

Figure 4. $\mu$ -parameter Bayesian estimate MCMC related image

图4. $\mu$ 参数贝叶斯估计值MCMC相关图像

Figure 5. $\sigma$ -parameter Bayesian estimate MCMC related image

图5. $\sigma$ 参数贝叶斯估计值MCMC相关图像

Figure 6. $\alpha$ -parameter Bayesian estimate MCMC related image

图6. $\alpha$ 参数贝叶斯估计值MCMC相关图像

Figure 7. $p_1$ -parameter Bayesian estimate MCMC related image

图7. $p_1$ 参数贝叶斯估计值MCMC相关图像

Figure 8. $p_2$ -parameter Bayesian estimate MCMC related image

图8. $p_2$ 参数贝叶斯估计值MCMC相关图像

5. 结论

本文针对非对称指数幂分布后验密度的复杂性，采用MCMC中的M-H算法对参数后验密度进行抽取。对AEP分布进行贝叶斯参数估计，假定先验分布为无信息先验，利用随机游动的M-H算法，对待估参数进行后验密度抽取，从MCMC模拟结果可以看到，各个分位点待估参数估计值的MCMC迭代轨迹图均收敛，这表明随机游动M-H算法很好的模拟了非对称指数幂分布的后验密度；得到的待估参数估计值与真值非常接近，说明该算法对具有非对称指数幂分布的数据集来说非常适用。

NOTES

^*通讯作者。