劉鶴飛

（曲靖師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，云南 曲靖 655011）

0 引言

回歸模型是最經(jīng)典的統(tǒng)計(jì)學(xué)模型，其廣泛應(yīng)用于經(jīng)濟(jì)學(xué)、管理學(xué)、心理學(xué)、教育學(xué)、醫(yī)學(xué)等領(lǐng)域。從最簡(jiǎn)單的一元線性回歸到多元線性回歸，再到非線性回歸、廣義線性回歸，眾多研究者對(duì)回歸模型進(jìn)行了深入的研究、探索和各種改進(jìn)。

隱馬爾可夫模型是一種基于隨機(jī)過(guò)程的統(tǒng)計(jì)模型。近年來(lái)，隱馬爾可夫模型在各個(gè)領(lǐng)域都得到了廣泛的應(yīng)用和發(fā)展。在經(jīng)濟(jì)管理領(lǐng)域，利用隱馬爾可夫模型處理異質(zhì)面板數(shù)據(jù)，在生物醫(yī)學(xué)領(lǐng)域，利用隱馬爾可夫模型對(duì)DNA序列的分布進(jìn)行推斷[1]；在人工智能領(lǐng)域，利用隱馬爾可夫模型進(jìn)行語(yǔ)音識(shí)別[2]、圖像處理等。

隱馬爾可夫模型是由兩個(gè)隨機(jī)過(guò)程構(gòu)成[3]。一個(gè)是狀態(tài)轉(zhuǎn)移序列，它是一條單純的馬爾可夫鏈，另一條是與狀態(tài)對(duì)應(yīng)的觀測(cè)序列。在實(shí)際問(wèn)題中，我們只能看到觀測(cè)變量的集合，無(wú)法看到觀測(cè)狀態(tài)序列的集合。隱馬爾可夫模型的研究?jī)?nèi)容就是根據(jù)可觀測(cè)的序列集合去推斷不可觀測(cè)的狀態(tài)轉(zhuǎn)移特征以及每個(gè)狀態(tài)下的分布信息[4]。

本文將隱馬爾可夫模型與回歸模型相結(jié)合，提出隱馬爾可夫回歸模型的概念。這種模型在許多領(lǐng)域都有實(shí)際的應(yīng)用，例如，在經(jīng)濟(jì)領(lǐng)域，由于股票指數(shù)與經(jīng)濟(jì)增速、CPI指數(shù)、銀行利率等指標(biāo)都有線性相關(guān)關(guān)系，并且可以建立這幾個(gè)自變量與因變量（股票價(jià)格指數(shù)）之間的多元線性回歸模型。但是，股票市場(chǎng)處于牛市和熊市的不同狀態(tài)時(shí)，因變量和自變量之間的回歸關(guān)系是不同的。此外，不同狀態(tài)之間也是可以相互轉(zhuǎn)化的，研究者還想了解不同的狀態(tài)之間相互轉(zhuǎn)化的關(guān)系。隱馬爾可夫多元線性回歸模型就是研究不同狀態(tài)之間的相互轉(zhuǎn)化規(guī)律以及每個(gè)狀態(tài)下因變量和自變量之間回歸關(guān)系的模型。

本文將以多元線性回歸模型為例，詳細(xì)介紹隱馬爾可夫多元線性回歸模型的數(shù)學(xué)定義，推導(dǎo)用貝葉斯方法對(duì)模型的參數(shù)進(jìn)行估計(jì)的理論過(guò)程。最后利用MCMC算法模擬參數(shù)的后驗(yàn)，用后驗(yàn)均值作為參數(shù)的估計(jì)值，并與模型參數(shù)的真值進(jìn)行比較，檢驗(yàn)該方法估計(jì)的可靠性。

1 模型描述

假設(shè)隱狀態(tài)的轉(zhuǎn)移過(guò)程滿足以下馬爾可夫鏈的條件：

其中u=1，2，…，K；s=1，2，…，K；t=2，3，…，T。

這里，aus是從前一個(gè)時(shí)間點(diǎn)的隱狀態(tài)u向后一個(gè)時(shí)間點(diǎn)的隱狀態(tài)s的轉(zhuǎn)移概率。我們稱所有可能的隱狀態(tài)轉(zhuǎn)移概率構(gòu)成的矩陣為隱狀態(tài)轉(zhuǎn)移概率矩陣，記為：

在給定隱狀態(tài)的條件下，描述自變量與因變量關(guān)系的多元線性回歸模型定義為：

其中，β是P維回歸系數(shù)向量，εk是模型的誤差，且：

2 貝葉斯推斷原理

將隱狀態(tài)轉(zhuǎn)移概率矩陣記為A；將多元線性回歸模型中的參數(shù)(βk，)記為θ；將隱狀態(tài)向量記為Z。則該模型的貝葉斯推斷問(wèn)題為，其中D=(Y，X)。

馬爾可夫鏈蒙特卡洛算法模擬過(guò)程的具體步驟如下：

（1）更新隱狀態(tài)Z；

（2）更新潛變量A；

（3）更新模型參數(shù)θ。

其中，更新模型參數(shù)θ又可以分為兩小步：

（1）更新正態(tài)分布的方差；

（2）更新多元線性回歸系數(shù)βk。

每一步更新參數(shù)，都是借助于Gibbs抽樣[6]和MH算法[7]。這需要對(duì)每個(gè)參數(shù)設(shè)定其先驗(yàn)分布，并推導(dǎo)出后驗(yàn)分布。

3 先驗(yàn)分布

貝葉斯理論認(rèn)為，每一個(gè)參數(shù)都是一個(gè)隨機(jī)變量[8]。因此，在進(jìn)行貝葉斯推斷時(shí)，必須事先為每一個(gè)參數(shù)都選擇一個(gè)合適的先驗(yàn)分布。Robert[9]在研究隱馬爾可夫正態(tài)模型時(shí)，選擇對(duì)稱的狄尼克萊分布作為轉(zhuǎn)移概率矩陣每一行的先驗(yàn)分布。Minka[10]在研究線性回歸模型參數(shù)的貝葉斯估計(jì)時(shí)，選擇多元正態(tài)分布和倒伽瑪分布作為線性回歸模型參數(shù)的先驗(yàn)分布。借鑒以上經(jīng)驗(yàn)，本文為模型中所有參數(shù)選擇的先驗(yàn)分布如下：

其中，Ak表示轉(zhuǎn)移概率矩陣的第k行，D表示狄尼克萊分布，α是狄尼克萊分布的超參數(shù)；μ0k，Λ0k是多元正態(tài)分布中的超參數(shù)，Inv-Gamma表示倒伽瑪分布，是倒伽瑪分布中的超參數(shù)。

4 后驗(yàn)推斷

貝葉斯后驗(yàn)推斷的主要任務(wù)是利用樣本的觀測(cè)信息和參數(shù)的先驗(yàn)信息推導(dǎo)出參數(shù)的后驗(yàn)信息[11]。對(duì)于本文的隱馬爾可夫多元線性回歸模型來(lái)說(shuō)，即也就是要根據(jù)觀測(cè)變量集合和模型參數(shù)的先驗(yàn)分布去推斷模型的隱狀態(tài)、概率轉(zhuǎn)移矩陣、多元回歸模型的參數(shù)。其中，轉(zhuǎn)移概率矩陣A、多元回歸模型的系數(shù)β是隱馬爾可夫多元線性回歸模型中的參數(shù)，是需要進(jìn)行估計(jì)的對(duì)象。由于這個(gè)分布的復(fù)雜性，本文使用了MCMC方法來(lái)進(jìn)行后驗(yàn)?zāi)M，用后驗(yàn)均值作為相關(guān)參數(shù)的估計(jì)。這就需要相關(guān)參數(shù)的全條件后驗(yàn)分布。

其中，πk是在轉(zhuǎn)移概率矩陣A的作用下，隱狀態(tài)達(dá)到穩(wěn)定時(shí)，隱狀態(tài)k的穩(wěn)定概率；是觀測(cè)變量在隱狀態(tài)k下的似然函數(shù)。

下面推導(dǎo)在隱狀態(tài)確定的條件下，多元線性回歸模型參數(shù)(βk，)的全條件后驗(yàn)分布。

記所有隱狀態(tài)為k的觀測(cè)數(shù)據(jù)集合為Dk，所有隱狀態(tài)為k的因變量集合記為Yk，所有隱狀態(tài)為k的自變量的集合為Xk。

在此記法下，Yk的似然函數(shù)可以表示為：

根據(jù)條件概率的計(jì)算公式可得：

于是，可以把參數(shù)(βk，的后驗(yàn)分布看成一個(gè)多元正態(tài)分布和一個(gè)倒伽瑪分布的乘積，它的具體形式為：

為簡(jiǎn)單起見(jiàn)，可將該后驗(yàn)分布看成如下兩部分：

則βk和這兩個(gè)參數(shù)的后驗(yàn)分布分別為N，即參數(shù)為μn和的多元正態(tài)分布和參數(shù)為an和bn的倒伽瑪分布。其中：

5 實(shí)證模擬

為了檢驗(yàn)上文所推導(dǎo)的模型參數(shù)的貝葉斯估計(jì)方法，這里將事先給定模型參數(shù)A和θ的真實(shí)值。根據(jù)概率轉(zhuǎn)移矩陣A生成一個(gè)隱狀態(tài)序列集合，再根據(jù)每個(gè)觀測(cè)點(diǎn)的隱狀態(tài)取值和對(duì)應(yīng)于該隱狀態(tài)的多元線性模型參數(shù)的值，生成每一個(gè)觀測(cè)時(shí)間點(diǎn)的觀測(cè)向量。然后利用所有的觀測(cè)變量集合，根據(jù)所推導(dǎo)的后驗(yàn)分布，用MCMC算法對(duì)模型的參數(shù)進(jìn)行后驗(yàn)?zāi)M。取后驗(yàn)均值作為模型中參數(shù)的估計(jì)值。最后將參數(shù)的貝葉斯估計(jì)結(jié)果與事先給定的真實(shí)值進(jìn)行比較，觀察估計(jì)的效果。

取隱狀態(tài)的個(gè)數(shù)K=2，則隱狀態(tài)的概率轉(zhuǎn)移矩陣是一個(gè)二階方陣，令：

設(shè)每個(gè)隱狀態(tài)下的多元線性回歸模型都有3個(gè)自變量，具體為：

模型中，取觀測(cè)時(shí)間點(diǎn)總數(shù)T=200，則生成的觀測(cè)變量集合是一個(gè)4×200的二維數(shù)組。先驗(yàn)分布中超參數(shù)的取值分別為：，其中I3表示三階單位矩陣

使用MCMC算法進(jìn)行后驗(yàn)?zāi)M時(shí)，取迭代總次數(shù)為5000，去掉前面3000次的迭代結(jié)果，用后面2000次的后驗(yàn)均值作為參數(shù)的估計(jì)。各參數(shù)的真實(shí)值與估計(jì)值如表1所示。

表1 模型參數(shù)的真值與估計(jì)值

6 結(jié)論

本文介紹了隱馬爾可夫線性回歸模型，推導(dǎo)了隱馬爾可夫多元線性回歸模型參數(shù)的貝葉斯估計(jì)方法。并且通過(guò)實(shí)證模擬，將模型參數(shù)的貝葉斯估計(jì)結(jié)果與事先設(shè)定的模型參數(shù)的真實(shí)值進(jìn)行比較，發(fā)現(xiàn)估計(jì)效果良好，這說(shuō)明本文給出的模型參數(shù)的貝葉斯估計(jì)方法是可靠的?？梢杂秒[馬爾可夫多元線性回歸模型來(lái)研究含有多個(gè)狀態(tài)的自變量與因變量的關(guān)系模型，分析不同狀態(tài)之間的轉(zhuǎn)化關(guān)系，以及每個(gè)狀態(tài)下多個(gè)自變量與因變量之間的線性回歸關(guān)系。

本文僅以簡(jiǎn)單的多元線性回歸模型為例介紹了隱馬爾可夫線性回歸模型，未來(lái)還可以研究隱馬爾可夫非線性回歸模型等更復(fù)雜的模型。此外，還可以研究隱狀態(tài)個(gè)數(shù)未知情形下的隱馬爾可夫回歸模型，利用貝葉斯因子、可逆跳躍MCMC算法等方法對(duì)隱狀態(tài)個(gè)數(shù)進(jìn)行模型選擇。