聯(lián)合均值和散度逆高斯回歸模型的參數估計

張露露，黃希芬

（云南師范大學 數學學院，昆明 650500）

0 引言

眾所周知，生產實踐中產生的數據在大多數情況下是正偏態(tài)數據。因此，對正偏態(tài)數據的分析一直是研究者們感興趣的課題，而這也符合現(xiàn)實需求。統(tǒng)計學中的許多參數模型都適用于正偏態(tài)數據的分析。例如，Gamma 分布族由于在數學上可得到顯式表達式而備受研究者們的喜愛；逆高斯分布族具有良好的概率和統(tǒng)計性質。鄧平穩(wěn)和謝治州（2022）[1]基于以上兩個回歸模型，通過對一組房價數據的分析來預測未來的房價走勢，結果表明兩種模型在預測不同問題時各有千秋。事實上，逆高斯分布因其復雜的表達式及計算過程而遠不如Gamma 分布受歡迎。然而，人們發(fā)現(xiàn)對于可靠性分析、精算學、生態(tài)學等領域中的一些數據，逆高斯分布族更符合其數據特點。

逆高斯分布起源于Brownian 運動問題。隨著人們對其了解漸漸深入，也衍生出不同的表達式[2]。人們可根據數據特點、參數估計算法不同等因素而對其進行不同的選擇和應用。例如，Chankham 等（2022）[3]基于逆高斯分布變異系數的置信區(qū)間對泰國PM2.5的擴散進行分析和研究；Amin等（2021）[4]基于響應變量服從逆高斯分布的假設，構造記憶型控制圖，將其應用于紗線制造業(yè)；黃順林等（2010）[5]基于零調整逆高斯回歸模型，對一組汽車責任險缺失數據進行了研究；莊會富等（2018）[6]提出將逆高斯分布和KI準則相結合，對影像進行非監(jiān)督變化檢測。除此之外，還有一群學者致力于發(fā)展各種算法來求解逆高斯回歸模型回歸參數的估計值。在均值的假設下，先后發(fā)展了嶺估計[7]、Liu估計[8]、KL估計[9]等方法；在均值的假設下，除了以上三種估計方法，還有迭代加權最小二乘估計法[1]；在均值μ=的假設下，發(fā)展了貝葉斯估計法[10]。顯然，在均值參數μ＞0 的必要條件下，第一種和第三種假設有明顯的局限性，即必須大于0。另一種回歸模型是聯(lián)合均值和散度模型，即假設均值參數和散度參數分別對應各自的回歸參數。具體地，假設均值參數，散度參數。在估計參數時，涉及Fisher信息矩陣[11,12]以及EM算法[13]。而對于以上兩種逆高斯回歸模型，除了貝葉斯估計法，上述其他方法都涉及矩陣求逆的問題，即算法的有效性依賴于矩陣是否奇異。

基于以上方法的局限性，本文立足于聯(lián)合均值和散度逆高斯回歸模型，將Minorization-Maximization（MM）算法[14]用于求解兩類回歸參數的極大似然估計。MM算法本質上是EM算法的一種延伸，但其在處理多元函數方面又具有EM算法不具備的優(yōu)勢，比如可以繞開多元函數求積分問題，通過合適的不等式把多元目標函數轉化為一系列低維函數之和，甚至可能把感興趣的未知參數彼此分離開，而在這之后求解函數最值問題時可繞開矩陣求逆的局限性。MM 算法由于其應用的靈活性以及計算簡便等優(yōu)點，已被廣泛應用于各種問題的研究，比如變量選擇[15]、分位數回歸[16]、半參數模型等。

1 MM算法原理簡介

假設已知l(θ|Yobs)為目標函數，其中，Yobs為數據觀測值，向量θ=(θ1，θ2，…，θq)′∈Θ，Θ是參數空間，θ是感興趣的未知參數。則參數向量θ的極大似然估計為θ^=。在極大化目標函數l(θ|Yobs)時，往往會涉及復雜的計算過程，其一階導函數在通常情況下也沒有解析解，這時常常要用到數值方法，如采用Newton-Raphson 迭代算法求其近似解，而這在多元函數情況下又涉及矩陣是否可逆的問題。因此，MM算法提供了另一種求解思路。MM算法主要分為兩個步驟：第一步是Minorization步驟，該步驟是至關重要的一環(huán)，其關鍵是尋找到合適的不等式，通過不等式放縮構造目標函數l(θ|Yobs)的一個極小化函數（或稱作替代函數）Q(θ|θ(k))，使其滿足：

其中，θ(k)是θ的第k次迭代值。若此時的替代函數還不盡如人意，則可把該替代函數看作新的目標優(yōu)化函數，繼續(xù)重復第一個步驟的思想。第二步是Maximization步驟，顧名思義即極大化替代函數，從而得到參數向量θ的極大似然估計值。

2 聯(lián)合均值和散度逆高斯回歸模型

2.1 模型介紹

假設響應變量Yi～Inverse Gaussian(μi，τi),且相互獨立，i=1，2，…，n，則其概率密度函數為：

其中，均值參數μi＞0 且滿足lnμi=X′i β，散度參數τi＞0 且滿足lnτi=,Xi=(xi1，xi2，…，xiP)′和Zi=(zi1，zi2，…，xiJ)′是解釋變量，β=(β1，β2，…，βP)′是均值模型的未知參數，η=(η1，η2，…，ηJ)′是散度模型的未知參數,β和η是待估計的回歸參數。

對應的對數似然函數為：

其中，c是與待估參數無關的項，可視作常數。l1(η|Yobs)，l2(β，η|Yobs)和l3(β，η|Yobs)分別為：

2.2 應用MM算法進行參數求解的步驟

本文的目的是把MM 算法應用到聯(lián)合均值和散度逆高斯回歸模型，求解兩類回歸參數的極大似然估計值。在開發(fā)MM 算法時，最大的困難在于第一步，即構造一個合適的極小化函數，而成敗的關鍵在于找到合適的不等式，這需要具體問題具體分析。在實際問題中，尤其是涉及高維模型時，MM 算法更吸引人的地方在于，可將替代函數構造為參數分離的低維函數之和，把高維優(yōu)化問題轉化為低維優(yōu)化問題，從而降低計算的復雜度。

本文希望找到對數似然函數l(β，η|Yobs)的一個替代函數，將參數向量β和η分離開，進而可以分別對其進行觀察和研究，然后再把各個分量從其向量中分離出來，即最終的替代函數形式是一系列一元函數之和，從而繞開了求解過程中的矩陣奇異性問題。

觀察式（2）可知，可把目標函數l(β，η|Yobs)重寫為三個函數（式（3）至式（5））的和。接下來，根據各個函數的特點逐一分離β和η，分別為其找到合適的替代函數。通過觀察可知，l1(η|Yobs)已與β無關，而l2(β，η|Yobs)和l3(β，η|Yobs)則需要借助MM算法分離β和η。對于，可利用算術-幾何均值不等式[17]將待估參數η和β分離，即利用不等式：

其中，c是與待估參數無關的項。

基于式（6）和式（7），可將函數l2(β，η|Yobs)和l3(β，η|Yobs)中的β和η分離開，從而得到只關于未知參數η的替代函數：

以及只關于未知參數β的替代函數：

觀察Q1(η|Yobs)和Q2(β|Yobs)可知，分別直接對其進行極大化是得不到解析解的，因此本文借助Newton-Raphson算法進行迭代求解。但在實際應用時，并不能保證矩陣總是非奇異的，為了避免遇到此類問題，接下來將繼續(xù)借用MM算法分別對參數向量η和β進行分量分離，將多元函數Q1(η|Yobs)和Q2(β|Yobs)分解為多個一元函數之和，從而避開求解極大似然估計時矩陣求逆時遇到奇異矩陣的問題。

當ωij=0 時，令。根據離散型Jensen 不等式[14]，構造如下不等式：

從而可將η1，η2，…，ηJ從函數Q1(η|Yobs)中分離開來，得到J個新的替代函數，即：

其中，Q1j(ηj|Yobs)是只關于參數ηj的一元函數。令?Q1j(ηj|Yobs)?ηj=0，利用Newton-Raphson 算法進行迭代求解，可得ηj的迭代公式為：

綜上所述，可得到如下迭代算法：

3 數值模擬

本文通過模擬研究來說明MM 算法在求解聯(lián)合均值與散度逆高斯回歸模型的參數極大似然估計值時的可行性和有效性。

模擬研究的隨機數據由以下模型產生：

為了避免兩類回歸參數因分量個數不同、真值取值不同而可能引起的模擬效果不同，設定P=J=5,兩類回歸參數向量取相同的值，即β=η=(-2，-1，0.5，1，2)′。接下來，通過設計不同的解釋變量值來分析兩類解釋變量的取值是否會影響均值模型對應的回歸參數和散度模型對應的回歸參數的估計結果。具體實驗設計如下：

實驗一：兩類解釋變量之間相互獨立，且組內也相互獨立，且Xip～U(0，1)，Zij～U(0，1)。

實驗二：兩類解釋變量共線，但組內相互獨立，即Xip=Zij，Xip～U(0，1)。

除此之外，為研究樣本量對兩組實驗結果的影響，兩組實驗的樣本量n均取{50，100，200，300}。整個模擬由R軟件實現(xiàn)，每種情況模擬500次。

表1 中的和是在實驗一的條件下，基于500 次模擬樣本進行MM算法求解得到的回歸參數估計值，取其樣本均值，進而求出樣本均值與真值之間的偏差（Bias），以及估計值的樣本標準差（SD）。隨著樣本量的增加，兩類參數的估計值總體上是越來越接近真值的。具體地，對比各個均值回歸參數的偏差絕對值，(β1，β5)=(-2，2)的偏差絕對值隨著樣本量的增加而一致地減小；(β2，β3)=(-1，0.5)的偏差絕對值總體上越來越小，但β2在樣本量n=100、β3在n=300 時的偏差絕對值有所增加，但與其他四組在數值上并沒有很大的差異；相比其他四個變量，β4=1 的估計值與樣本量之間并沒有以上關系，雖然其最大的偏差絕對值不超過0.05，但相比其他變量，這個結果是不理想的。接下來，對比各個散度回歸參數的偏差絕對值，前三個回歸參數的偏差絕對值隨著樣本量的增加而一致地減小；η4=1在n=100 時，其偏差絕對值增加，其偏差最大值不超過0.01；除了(n，η5)=(200，2)這組結果波動較大外，η5的偏差絕對值總體上是越來越小的。在樣本量相同的條件下，對比兩組參數的偏差絕對值，絕大多數情況下都是均值回歸參數比散度回歸參數的估計效果更好。最后，對比樣本標準差可以看出，兩組參數都是隨著樣本量的增加而一致地減小，但依然是前者有著更小的數值。

表1 實驗一中兩組參數向量的每個參數分量的模擬結果

表2是在實驗二的條件下對500次模擬樣本進行MM算法求解得到的數值結果。隨著樣本量的增加，兩類參數的估計值總體上是越來越接近真值的。首先，β1、β2、β3在n=300 以及β5在n=200 時對應的偏差絕對值有所增加，但總體趨勢還是樣本量的增加有益于模擬效果的優(yōu)化；β4=1 的估計值與樣本量之間并沒有以上的趨勢，其最大的偏差絕對值不超過0.04。前三個散度回歸參數的偏差絕對值隨著樣本量的增加而一致地減小；η4=1 在n=200 時偏差絕對值增加，在n=300 時顯著減小；η5的偏差絕對值隨著樣本量的增加而減小，但在n=300 時有所增加。然后，對比兩組參數的偏差絕對值，從n=100 開始，隨著樣本量的增加，均值回歸參數比散度回歸參數的估計效果一致地更好。最后，對比樣本標準差可以看出，兩組參數隨著樣本量的增加而一致地減小，但前者模擬效果更好。

表2 實驗二中兩組參數向量的每個參數分量的模擬結果

表1 和表2 展示了對每個分量進行對比研究的結果，而表3 則通過均方誤差（MSE）來研究每組回歸參數的整體模擬效果。結果表明：首先，隨著樣本量的增加，每組實驗的模擬效果都越來越好，但顯然均值回歸參數比散度回歸參數具有更好的效果；其次，實驗一的結果整體上比實驗二的效果好，即數據之間獨立比共線具有更好的結果，這也符合以往的理論和實踐；最后，不管數據之間的相關性如何，當樣本量從50 增加到100 時，整體模擬效果都顯著變好，因此在實際應用中，當數據量達到100 時就可以得到一個很接近實際情況的結果。

表3 兩組參數向量分別在兩組實驗中的MSE

綜上所述：首先，均值回歸參數的模擬效果總體上比散度回歸參數的模擬效果更好；其次，兩組解釋變量之間是否共線對結果的影響并不顯著；最后，樣本量的增加可以改善模型的結果，且在實際應用時，當數據量達到100時就可以得到一個很好的估計結果。

4 實例分析

本文用一組實際數據檢驗該方法的實用性。該數據集[12,18]共收集了125 組數據，涉及5個變量（見表4）。具體地，先記錄下汽油罐的初始罐溫以及初始氣壓，當該油罐裝滿汽油后，再記錄下其溫度和氣壓，以及設備回收的碳氫化合物重量。將以上實驗重復125 次，則可收集到125組數據。

表4 變量說明

把設備回收的碳氫化合物重量設為響應變量Y,且Yi～Inverse Gaussian(μi，τi)，i=1，2，…，125，再將剩下4個變量設為解釋變量，具體的變量說明見表4。已有研究[12,18]通常只取其中的32組數據進行建模求解，但本文的數值模擬結果表明，當數據集達到100組時得到的結果更接近實際情況，所以對全部125組數據進行建模求解。本文的目的是研究一組解釋變量對逆高斯分布的均值參數和散度參數是否有影響，以及哪些變量有顯著影響，故模型假設X=Z。

為了檢驗參數估計值的顯著性，本文采用Bootstrap法[19]對數據集進行抽樣，得到一組Bootstrap 樣本集，再用MM算法估計該組樣本集下的回歸參數，重復該過程B次得到B組參數結果。記。若每個參數的B個估計值具有正態(tài)性，則可構造如下的t分布置信區(qū)間：

若每個參數的B個估計值不具有正態(tài)性，則可構造如下的百分位置信區(qū)間：

本文基于逆高斯聯(lián)合均值和散度回歸模型，先運用MM算法求解得到參數的極大似然估計值（MLE），再利用上述Bootstrap 方法得到B=20000 組參數估計值，最后基于顯著性水平α=0.05 分別計算得到t分布置信區(qū)間以及百分位置信區(qū)間，具體結果見表5。對比估計值，發(fā)現(xiàn)兩組回歸參數有著不同的估計值，其中解釋變量(X3，X4)有著相反的符號。觀察t分布置信區(qū)間，發(fā)現(xiàn)汽油溫度X2、初始氣壓X3及汽油氣壓X4都對逆高斯分布的均值有顯著影響；初始罐溫X1和汽油氣壓X4對散度有顯著影響，該結論與已有的研究結果[19]一致。再觀察百分位置信區(qū)間，發(fā)現(xiàn)汽油溫度X2以及汽油氣壓X4都對逆高斯分布的均值有顯著影響；初始罐溫X1、初始氣壓X3以及汽油氣壓X4對散度有顯著影響。

表5 參數估計值及其兩種Bootstrap置信區(qū)間

5 結束語

響應變量服從逆高斯模型的數據在實踐中普遍存在，許多學者也致力于這方面的研究。基于不同的現(xiàn)實背景，人們會著眼于不同的均值回歸參數形式，也有研究散度回歸模型的，但同時對均值和散度參數進行回歸分析的研究相對較少。本文對聯(lián)合均值和散度逆高斯回歸模型進行了探討，在求解模型參數的極大似然估計時，由于逆高斯概率密度函數的表達式復雜，沒有顯式解，故需借助Newton-Raphson 算法等數值方法，這時求解多元函數值時總是涉及矩陣求逆問題，而要求現(xiàn)實數據總是滿足矩陣可逆是不現(xiàn)實的。本文基于MM 算法可以進行參數分離的優(yōu)良性質，繞開了矩陣奇異問題，為求解回歸參數估計值提供了一種新的可能。本文的MM 算法也可以解決只考慮均值或散度具有回歸參數形式的逆高斯回歸模型。數值模擬及實驗分析驗證了所提方法的可行性。