假设向量D_t=(y_t, x_t1, x_t2, …, x_tr)^T是模型在第t个时刻的观测向量, 其中：y_t是第t个时刻的因变量; x_t1, x_t2, …, x_tr是在第t个时刻的r个自变量; t为观测时间点, t=1, 2, …, T.记Z_t是第t个观测时间点的隐状态, 其可能的取值为{1, 2, …, K}, 称之为K个隐状态的隐马尔科夫模型.

隐状态转移过程满足的马尔科夫链条件如(1)式：

其中：u=1, 2, …, K, s=1, 2, …, K, t=2, 3, …, T; a_us表示从前一个观测时间点的隐状态u向后一个观测时间点的隐状态s转移的概率.所有可能的隐状态转移概率构成的矩阵称之为隐状态转移概率矩阵：

在隐状态给定的条件下, 自变量与因变量之间的多元线性回归模型表示如下：

其中：β为r维回归系数向量, ε_k为模型误差且[ε_k|Z_t=k]~N(0, σ_k²).

在贝叶斯理论体系中, 每个参数都是一个随机变量, 为每个参数选择先验分布, 是进行贝叶斯推断的前提^[10].文献[11]等对隐马尔科夫正态分布进行研究时, 运用对称的狄尼克莱分布作为状态转移概率矩阵中每一行的先验分布.文献[12]等对线性回归模型参数进行贝叶斯估计时, 运用多元正态分布和倒伽玛分布作为模型中参数的先验分布.本文参考已有文献中的研究经验, 选择(4)式为隐马尔科夫多元线性回归模型中相关参数的先验分布.

其中：A_k为转移概率矩阵的第k行, Diri为狄尼克莱分布, α表示狄尼克莱分布中的超参数; μ_ok与Λ_ok为多元正态分布中的超参数, Inv-Gamma为倒伽玛分布, υ_ok与s_ok²为倒伽玛分布中的超参数.

将模型中的隐状态转移概率矩阵记为A, 多元线性回归模型中的参数(β_k, σ_k²)记为θ_k, (θ₁, θ₂, …, θ_k)记为θ, 隐状态向量记为Z.在隐状态个数固定的情况下, 模型中的贝叶斯参数估计问题可以表示为P(Z, A, θ|D), 其中D=(Y, X).

由于P(Z, A, θ|D)的复杂性, 无法直接从中推导出每个参数的边际后验分布.因此, 本文选择使用MCMC算法来模拟该分布, 使用收敛之后的后验均值作为每个参数的贝叶斯估计. MCMC算法迭代的具体步骤如下：

① 抽样, 对隐状态Z进行更新;

② 抽样, 对隐状态转移概率矩阵A进行更新;

③ 抽样, 对线性回归模型中参数θ进行更新.

其中, 参数θ的抽样更新可进一步划分为：

① 抽样, 对误差正态分布的方差σ_k²进行更新;

② 抽样, 对多元线性回归系数β_k进行更新.

每一步参数的抽样更新, 都是借助于Gibbs抽样^[13]和MH算法^[14]从参数的全条件后验分布中取样^[15].

在混合模型的研究中, 模型识别是一个非常重要的问题, 隐马尔科夫模型的研究也同样存在这个问题.因为如果不处理好模型识别问题, 进行MCMC算法迭代处理时, 不同隐状态下的参数有可能发生混淆, 这样计算出的后验均值就没有意义了.目前最常见的解决模型识别问题的方法是交换抽样技术^[16], 但是也有很多研究通过对模型中的某个特定参数进行限制来解决这个问题^[17].本文的研究中, 采用比较第一个自变量的回归系数大小的方法, 即β₁₁ < β₂₁ < … < β_K1来解决标签转换问题.暂不考虑第一个自变量的回归系数有相同取值的情况.实际应用中如需考虑该情况, 可以参考文献[18]使用的方法.

贝叶斯后验推断的任务是运用样本的观测信息和参数的先验信息推断参数的后验信息.隐马尔科夫模型在隐状态个数确定的情况下, 后验推断为P(Z, A, θ|D), 也即根据观测变量集合和模型参数的先验分布去推断模型的隐状态、隐状态转移概率矩阵、多元线性回归模型的参数.其中, 转移概率矩阵A、多元线性回归模型的系数β是隐马尔科夫多元线性回归模型中的参数, 是需要进行估计的对象.然而由于P(Z, A, θ|D)后验分布的复杂性, 无法直接从中推导出每个参数的边际后验分布, 因此本文选择使用MCMC算法对该后验分布进行模拟, 运用Gibbs抽样、MH算法从每个参数的全条件后验分布中抽样.

隐状态向量Z中每一个隐状态Z_t=k的全条件后验分布如下：

即

其中f(D_t|θ_k)是观测变量在隐状态k下的似然函数.

由于隐状态转移概率矩阵只与隐状态的取值有关, 即P(A_k|Z, θ, D)=P(A_k|Z), 所以隐状态转移概率矩阵第k行的全条件后验分布如下：

其中n_kj(j=1, 2, …K)是前一个时间点隐状态为k且后一个时间点隐状态为j的转移个数.

下面推导在隐状态确定的条件下, 多元线性回归模型参数(β_k, σ_k²)的全条件后验分布.记所有隐状态为k的观测数据集合为D_k、因变量集合为Y_k、自变量集合为X_k, 且记$ {{\mathit{\boldsymbol{\tilde \beta }}}_k} = {\left( {\mathit{\boldsymbol{X}}_k^{\rm{T}} \cdot {\mathit{\boldsymbol{X}}_k}} \right)^{ - 1}}\mathit{\boldsymbol{X}}_k^{\rm{T}} \cdot {\mathit{\boldsymbol{Y}}_k}$, 则

Y_k的似然函数可以表示为

其中：${v_k}s_k^2 = {\left( {{\mathit{\boldsymbol{Y}}_k} - {\mathit{\boldsymbol{X}}_k}{{\mathit{\boldsymbol{\tilde \beta }}}_k}} \right)^{\rm{T}}}\left( {{\mathit{\boldsymbol{Y}}_k} - {\mathit{\boldsymbol{X}}_k}{{\mathit{\boldsymbol{\tilde \beta }}}_k}} \right) $且v_k=n_k－r, n_k为隐状态k的观测样本数, r为自变量个数.根据条件概率计算公式可得：

记${\mathit{\boldsymbol{\mu }}_n} = {\left( {\mathit{\boldsymbol{X}}_k^{\rm{T}}{\mathit{\boldsymbol{X}}_k} + {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{ok}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{X}}_k^{\rm{T}}{\mathit{\boldsymbol{X}}_k}{{\mathit{\boldsymbol{\tilde \beta }}}_k} + {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_{ok}}{\mathit{\boldsymbol{\mu }}_{ok}}} \right) $, 则

由此, 可将参数(β_k, σ_k²)的后验分布看成一个多元正态分布和一个倒伽玛分布的乘积, 具体形式为：

为简化表达, 可将该后验分布如下表示：

则参数β_k和σ_k²的后验分布为N[μ_n, σ_k²Λ_n^-1], Inv-Gamma[a_n, b_n], 即参数为μ_n和σ_k²Λ_n^-1的多元正态分布和参数为a_n和b_n的倒伽玛分布, 其中

除了模型的参数估计外, 模型选择也是统计分析的一个重要部分.在贝叶斯理论中, 常用的模型选择方法主要是贝叶斯因子和各种信息准则, 如AIC, BIC, DIC等.本文利用经典的DIC信息准则对隐马尔科夫多元线性回归模型中未知的隐状态个数进行模型选择.

由于DIC_k=D(θ_k)+d_k, 其中

故DIC_k=-4E_θk[logP(Y|θ_k, M_k)|Y]+2logP(Y|${{\mathit{\boldsymbol{\tilde \theta }}}_k} $).

本文比较的是隐状态个数不同的各个模型, 因此将m个隐状态的模型记为M_m, 则观测变量在备选模型下的对数似然函数为

其中

在MCMC迭代的过程中, E_{θk, Mm}[logP(Y|θ_m, M_m)|M₀]的具体计算方法为

其中J是MCMC迭代中计算后验均值的次数.

然而, 对于隐马尔科夫模型来说, 计算DIC的困难之处在于计算logP(Y|${{\mathit{\boldsymbol{\tilde \theta }}}_k} $).因为MCMC算法的执行过程中, 每次迭代时都会对每个观测向量D=(Y, X)重新分配隐状态, 因此, 单纯用每个参数的后验均值来计算logP(Y|${{\mathit{\boldsymbol{\tilde \theta }}}_k} $)进而计算该备选模型的DIC是没有意义的.这也是利用DIC对隐马尔科夫模型进行模型选择的困难之处.文献[19]在利用DIC信息准则对混合结构方程模型中未知的混合个数进行贝叶斯模型选择时, 在MCMC迭代的每一步, 对参数采取重复抽样取均值的方法来计算logP(Y|${{\mathit{\boldsymbol{\tilde \theta }}}_k} $), 其具体方法为：在第j次MCMC迭代时, 对参数θ_k抽样L次, 计算参数θ_k的均值, 并用该均值作为${{\mathit{\boldsymbol{\tilde \theta }}}^{\left( j \right)}} $代入计算logP(Y|${{\mathit{\boldsymbol{\tilde \theta }}}^{\left( j \right)}} $), 即

则

于是, 第m个备选模型的DIC为

为了检验本文所介绍的隐马尔科夫多元线性回归模型的贝叶斯参数估计和模型选择的效果, 将事先给定模型中的隐状态个数K, 隐状态转移概率矩阵A和每个隐状态下多元线性回归模型的参数θ_k.首先根据转移概率矩阵A生成一个隐状态序列, 再依据每个观测时间点的隐状态取值与对应于每个隐状态的多元线性回归模型参数的取值, 生成每一个观测时间点的观测向量.然后利用本文所介绍的方法, 结合观测向量集合, 对待定的隐状态个数进行贝叶斯模型选择, 并对固定隐状态个数下模型中的参数进行贝叶斯估计^[20].最后再将模型选择与参数估计的结果同真实模型进行比较, 检验该方法的可靠性.

设隐状态个数取值K=3, 则隐状态的概率转移矩阵为如下三阶方阵

假设每个隐状态下的多元线性回归模型都有如下3个自变量：

其中：ε₁：N(0, 0.5²), ε₂：N(0, 0.4²), ε₃：N(0, 0.3²).

实证模拟中, 取观测时间点的总数t=500, 则生成的观测变量集合是一个4×500的二维数组.先验分布中超参数α, μ_ok, Λ_ok的取值分别为α=1, μ_ok=(0, 0, 0)^T, Λ_ok=I₃, 其中：I₃表示三阶单位矩阵; v_ok=8, s_ok²=2.5.所有的超参数中k=1, 2, …, K, 其中K是隐状态的个数.

运用MCMC算法进行后验迭代模拟时, 总迭代次数为5 000次, 去掉前3 000次的迭代结果, 使用后2 000次结果的后验均值作为参数的估计值, 则在模型选择中, j=2 000.此外, 在利用DIC进行模型选择时, 取L=5.隐状态个数不同的各个模型下的DIC计算结果如表 1所示.

根据表 1中隐状态个数不同的各个模型的DIC计算结果可以发现, M₃的DIC值最小, 所以选择隐状态个数为3的隐马尔科夫多元线性回归模型.在3个隐状态个数下, 模型中各参数的贝叶斯估计结果如表 2所示, 表明本文的贝叶斯推断方法是有效可靠的.

本文对隐状态个数未知时的隐马尔科夫多元线性回归模型的贝叶斯推断进行了研究.首先, 运用了偏差信息准则对未知的隐状态个数进行了贝叶斯模型选择; 然后, 运用马尔科夫链蒙特卡洛算法对模型中的未知参数进行贝叶斯估计; 最后, 按照所定义的模型生成了观测向量集合, 并用实证模拟证明了文章所介绍的贝叶斯参数估计和模型选择的方法是可靠的.在进行贝叶斯模型选择时, 本文使用的是经典的DIC信息准则, 今后还将尝试使用贝叶斯因子、赤池信息量准则(AIC)、贝叶斯信息准则(BIC)等方法来对隐马尔科夫多元线性回归模型的隐状态个数进行模型选择.此外, 本文研究的是隐马尔科夫多元线性回归模型, 今后也可以研究更复杂的模型, 比如, 隐马尔科夫半参数回归模型、隐马尔科夫分位数回归模型等.