貝葉斯推斷下部分線性模型中非參數部分的估計基于重構參數化方法

中圖分類號：0212.8 文獻標志碼：A

Estimation of the Nonparametric Component of Partial Linear Models under Bayesian Inference

Based on Reconstructed Parametric Methods

ZHAO Bo-han，YANG Jian-kui （School of Science，Beijing University of Posts and Telecommunications，Beijing lOo876，China）

Abstract： When dealing with the nonparametric component of partial linear models，traditional Gaussian process prior methods，while capable of estimating the nonparametric component，have low computational efficiency and are not suitable for handling high-dimensional data with large sample sizes. To address this issue，a reparameterized method was employed to interpolate and reconstruct the nonparametric component of partial linear models. New prior distributions were assigned to the parameters obtained after reconstruction. Based on these new priors，Bayesian inference was then utilized to derive the closed-form posterior distributions of the parameters. Numerical simulation results demonstrate that the method proposed in this article can effectively reduce computational costs andachieve satisfactoryinferential outcomes.

Keywords： partially linear model; Bayesian inference； Gaussian process; reconstruction parameterization

簡單的線性模型往往不能很好地擬合各領域中的復雜數據，而非參數回歸模型雖然可以較好地擬合數據但缺乏可解釋性，因此在實際的應用中經常受阻。為解決這一問題，1986年Engle等在研究電力與氣候環境之間關系的問題時提出了部分線性模型[1]，包含參數和非參數兩部分，使模型兼具可解釋性和靈活性，并避免了非參數回歸模型中一些棘手的問題，例如“維度災難”。對于部分線性模型，特別是對其中非參數部分的估計問題，已有諸多文獻做了研究。Hardle提出利用局部線性估計方法得到非參數回歸分量的估計[2]，并且得到了漸進分布，利用樣條平滑方法可得到非參數分量以及線性參數的估計[3」，對非參數分量采用參數化形式逼近的方法研究了估計量的漸進性質[4」，此后提出了對參數部分最小二乘估計以及對非參數分量的 N-W核估計方法[5]和投影樣條估計非參數函數的方法[6]，基于核方法通過密度估計和數值反演構造了單調函數估計部分線性模型的非參數部分[7]。有研究基于核估計并結合小型蒙特卡洛探討了部分線性模型[8]，比如基于LASSO的部分線性單指標模型局部懲罰樣條估計方法，并利用蒙特卡洛模擬得到非參數的估計[9]，在研究部分線性變系數模型的估計時，可利用B樣條函數近似變系數函數方法[10]。貝葉斯推斷方法的不斷發展，為解決統計問題提供了新視角，為推理部分線性模型非參數分量開拓了新方法，比如在參數估計時使用約束貝葉斯方法[1]。高斯過程回歸（Gaussian Process Regression，GPR）首先在 1978 年提出，是基于高斯過程先驗及其協方差函數的貝葉斯非參數曲線擬合和分類的一種替代方法[12]，后有學者給出了GPR 模型的基本介紹[13]。高斯過程通過調節已知輸入和目標，推斷觀測位置的未知目標值[14]，是定義靈活貝葉斯回歸和分類模型的最簡單和最明顯的方法[15]。GPR 結合貝葉斯的方法在概率論和統計領域得到廣泛應用，實驗表明在極限的線性模型中，該方法的實際優點遠遠超出線性模型[16]。給予非參數分量高斯過程先驗并結合貝葉斯推斷得到的模型后驗推斷[17]，為解決部分線性模型中非參數問題提供了新思路，但該方法并不能很好地解決海量數據中非參數的估計。為解決算法的后驗計算問題，使用高斯過程先驗的貝葉斯推斷分析變異系數模型，將數據分解成數據子樣本，并將子樣本的后驗抽樣作為替代方案[18]，但在實驗方面比較復雜。作為統計方法中重要的函數逼近方法，插值方法在非參數回歸問題上有著廣泛的應用。插值重構參數化是將插值節點處的函數值參數化，從而將非參數回歸中的未知函數轉化為具有有限未知參數的形式，即將非參數問題轉化為參數問題[19]。插值方法產生的誤差非常小，并且可以通過選擇插值結點方式避免因為數據量龐大而造成的大量計算問題。本文使用高斯過程結合貝葉斯推斷方法和基于插值重構參數化結合貝葉斯推斷這兩種方法估計部分線性模型的非參數部分，計算得到部分線性模型中非參數部分以及參數部分的后驗分布，并利用數值實驗得到兩種方法中部分線性模型的非參數部分的估計結果。

1 部分線性模型

部分線性模型（PartiallyLinearModel，PLM）一般形式為

y_i=x_i^Tβ+f（u_i）+ε_i，i=1，…，n

其中， y_i 是對應的響應變量， x_i=（x_i1，…，x_ip）^T 是 p 維的解釋變量， u_i=（u_i1，…，u_ik）^T 是 k 維解釋變量，{u_i}_i=1ⁿ∈[0，1]^k 和 {x_i}_i=1ⁿ 均是固定并且已知的協變量， {x_i}_i=1ⁿ 與 {u_i}_i=1ⁿ 無關； f（u_i） 是未知函數， β= （β₁，…，β_p）^T 是未知線性回歸參數。在部分線性模型中， {ε_i}_i=1ⁿ 是獨立同分布的隨機誤差，每個 ε_i 都服正態分布 N（0，σ²） ， σ² 是未知參數。

假設 f_n=（f（u₁），…，f（u_n））^?T 是基于 u₁，u₂，…，u_n 的一組函數值，令 X=（x₁，x₂，…x_n）^T ， ε=（ε₁ …，ε_n）^T ， 是響應值的 n 維向量，即 y=（y₁，y₂，…，y_n）^T 。式（1）的矩陣形式為

y=Xβ+f_n+ε

2 高斯過程先驗方法

在不指定輸入條件下，部分線性模型中的未知函數為 f（?） 。給予未知函數高斯過程先驗，并且高斯過程由均值函數 m（t） 以及協方差函數 C（?，?） 唯一決定，即 f（?）～GP（m（t），C（?，?） ）。由于兩個參數β 和 σ² 的所有先驗信息未知，因此，可以使用無信息先驗分布的方式推理兩個參數。計算得到部分線性模型中參數的條件后驗分布以及非參數分量的后驗推斷，最后結合Gibbs抽樣得到非參數分量以及參數的后驗樣本。

2.1 先驗分布

假設 f（?） 具有零均值的高斯過程先驗，并且任意的 u_i∈{u_i}_i=1ⁿ 都有對應的 f（u_i） ）隨機變量序列，則

f（?）～GP（0，C（?，?））

根據定義有 C（u₁，u₂）≡Cov（f（u₁），f（u₂）） ，一般由超參數 θ=（σ_f²，η₁，…，η_k）^T 參數化且具有一些已知的非線性形式。在實際應用中，高斯協方差函數具有平滑性并且能夠很好地擬合非線性數據，在研究時間序列分析、機器學習算法以及圖像識別技術等領域時大多使用高斯協方差函數[20]

由 f（?） 的先驗假設， f_n 的先驗分布為 n 維的多元正態分布， C_n 是應用協方差函數 C（?，?） 得到的協方差矩陣

[f_n]～N_n（0，C_n）

對 β 采用的無信息先驗分布方法時，假設 β 的先驗分布服從均勻分布，即 β 的概率密度函數 π[β]∝1 。假設參數 σ² 具有Jefferys無信息先驗分布，即 σ² 的概率密度函數 π[σ²]∝1/σ² 。 β 和 σ² 互相獨立，則 β 和σ² 的聯合先驗概率密度函數

π[β，σ²]∝1/σ²

其中， β 和 σ² 都與 f（?） 無關。令 I_n 代表 n×n 的單位矩陣，結合先驗假設，得到在給定 f_n?β，σ² 的條件下 的條件分布 為

[y|f_n，σ²，β]～N_n（Xβ+f_n，σ²I_n）

2.2 后驗推斷

首先推導 f_n 后驗條件分布。由式（5）和式（7）得到 和 f_n 的聯合分布服從 n+n 維的多元正態分布

其中 ，Σ_ff=C_n，Σ_fy=C_n，Σ_yf=C_n^T，Σ_yy=σ²I_n+C_n， 。

在推斷 f_n 的后驗條件分布時，引入引理1。

引理1 假設 W₁ 和 W₂ 是隨機向量，且 W₁ 和 W₂ 的聯合分布服從 m+n 維的多元正態分布，記作 。其中， μ₁ 和 μ₂ 分別是 W₁ 和 W₂ 的均值向量， Σ₁₁ 是 m×m 的協方差矩陣， Σ₁₂=Σ₂₁^T 是 m×n 的協方差矩陣， Σ₂₂ 是 n×n 的協方差矩陣，則有條件分布 [W₁|W₂]～N_m（μ₁+Σ₁₂ （204號 。

根據引理1， f_n 的后驗條件分布為

[f_n|y，β，σ²]～N_n（μ_n，Σ_n）

推斷 β 和 σ² 兩個參數的后驗條件分布。式（6）和式（7）相乘，可得 β 和 σ² 的聯合后驗分布

根據式（10）， π[σ²∣y，β，f_n] 有

式（11）中，根據逆伽馬分布概率密度函數的表達式并利用逆伽馬分布（Inverse Gamma，IG）具有歸一性，得

弋人式（10）的分母中，得到 σ² 服從形狀參數為 0.5n ，尺度參數為 0.5（y-Xβ-f_n）^T（y-Xβ-f_n）^T 的】

[σ²σ|y，βγ，f_n]～IG（0，5n，0，5（y-Xβ-f_n）^T（y-Xβ-f_n））

求解參數 β 的后驗條件分布時，引入 ，其中， f_n ）是 β 的最小二乘估計。整理式（10）有

兩邊積分，得

由于 β 的后驗概率密度函數滿足

將式（13）和式（14）帶入式（15），整理得到 β 的后驗分布為

在超參數已知條件下，通過馬爾科夫鏈蒙特卡羅算法（Markov ChainMonteCarlo，MCMC）中的Gibbs抽樣方法獲得各個參數以及 f_n 的后驗樣本。Gibbs抽樣方法的實施依賴于各參數的完整條件分布，這些分布由式（9）、式（12）和式（16）給出。為了獲取有效樣本，在正式開始收集樣本之前，通常需要一個“燒入期（burn）\"以確保馬爾可夫鏈能穩定并收斂到目標分布。取總抽樣次數 M=1 000 ，燒入期 。

Gibbs抽樣算法采樣步驟：

初始化：隨機選擇 β^（0），（σ²）^（0），f_n^（0） 作為初始值迭代更新：令 k=M+burn

Step1將 （σ²）^（k-1） β^（k-1） 代入式（9），計算得到 μ_n^（k-1） 及 Σ_n^（k-1） ，由此得 f_n^（k）

Step2將 代人 ，由式（16），得 β^（k）

Step3將 與 β^（k） 代人式（12），得 采樣：當 kgt;burn 開始第一次采樣；循環：若 k 1 β^（k） 和 （204號

2.3 超參數估計

當協方差函數中的超參數未知時，應用極大似然法（MaximumLikelihood Estimation，MLE）估計協方差函數中超參數0。由式（5）和式（7），有 [y∣β，σ²，θ]～N_n（Xβ，Σ_GP） ！ 。除去可加的有限個常數項后，對數似然滿足正比關系后，得到超參數 θ 的極大似然估計，結合Gibbs抽樣得即可得 β，σ² 和 f_n 的后驗樣本。然而，給予未知函數高斯先驗結合Gibbs 抽樣的方法雖然可以解決部分線性模型中非參數部分的估計問題，但是這種方法計算成本較高。在數據維度較高，數據量較大的情況下，使用抽樣算法的效率非常低。針對該問題，本文提出基于插值重構參數化的貝葉斯推斷方法。

3基于重構參數化的貝葉斯推斷方法

利用重構參數化方法插值重構部分線性模型當中的未知函數，應用插值器對未知函數 f（?） 建模，使其具有可以進行貝葉斯推斷的新參數向量，給予各個參數先驗分布，經過計算得到各參數封閉形式的后驗分布，由此解決部分線性模型中的估計問題。

3.1重構參數化方法

非參數回歸是統計學和機器學習中的重要問題之一，重構法是多元非參數回歸的一種通用方法。基于插值重構參數化方法[19]利用插值器將有限插值結點處未知函數的值參數化，然后通過最小二乘估計得到這些值。由于結點 Σ_m 的選擇具有靈活性，當 m?n 時插值收斂通常比統計估計收斂到基本函數的速度快得多，因此適用于大規模的復雜數據集，可大大減少計算負擔。基于相關研究[19]，利用插值重構法對部分線性模型中的未知函數建模，將新參數賦予先驗并結合貝葉斯推斷解決未知函數的估計問題。

給定一組訓練數據 {（x₁，y₁），…，（x_n，y_n）} 以及定義在 D?R^d 上的光滑未知函數 f ，利用訓練數據估計 f 。在定義域內選擇合適的插值器 I（?;A，γ） ，其中， A 是插值結點集，即 A={a₁，…，a_m}?D ，對應的函數值為 γ=（γ₁，…，γ_m）∈R_m 。令 ，則 ，其中， （204號 ）是 γ 基于訓練數據的估計值，可由最小二乘法得到： P（γ，λ），P 是懲罰函數， λ 是調優參數，待估計參數就是選定結點處的參數值 γ 。在重構參數化過程中，使用高斯核（相關）函數 φ（?，?）

此時，核插值器又稱為高斯核插值器。指定一組已知的回歸函數 g（x）=（g₁（x），…，g_k（x））^? ，為便于計算，使用的回歸函數為 g（x）=（1，x^T）^T 。高斯核插值器： I（x;A，γ）=γ^Tb（x） ，其中， b（x）=Pg（x）+ Qr（x） γ;P=R_A^-1G_A（G_A^?R_A^-1G_A）^-1，Q=[I_n-R_A^-1（G_AG_A^?R_A^-1G_A）^-1G_A^?]R_A^-1;G_A=（g（a₁），…，g（a_n））^? ，rs（x）=（＄（x-a）），，＼（x-am））T，R=（a;-aj）i，j=1.，m。

在數據量龐大的情況下，選擇的結點個數 m?n 可以大大減少計算負擔。在數據點的數量 Ω_n 趨近于無窮大時，需要結點數量 Ψ_m 也趨近于無窮大，從而可以確保插值器收斂。通常取 m=max{10d，∣logn∣} ，其中， d 代表輸人數據的維度。選取結點時，可以選擇在所有 Ψ_m 子集中通過最小化準則的結點子集： c（A）=

3.2 先驗分布

根據上述方法對部分線性模型中的未知函數進行重構參數化，并構建用于貝葉斯推斷的重構參數先驗方法。重構參數化的表達式為： f（?）=γ^Tb（?） i=1，…，m ，其中， γ 和 b（?） 同式（21）， γ=（f（a₁） ，…，（20 f（a_m）^T 是未知函數插值結點處的函數值向量;插值結點集 A={a₁，…，a_m} 選擇最小化準則下的 m 子集。為了得到各參數封閉形式的后驗分布， γ 的先驗分布為： γ～N_m（μ_γ，σ²V） ，其中， ?V 是 m×m 的對角矩陣，對角線的值為式（3）定義的高斯協方差函數 C（a_i，a_j）_i，j=1^m 矩陣中對角線上的值。

令 τ=（σ²）^-1 ，服從伽馬先驗分布： [τ]～G（a，q） ？ agt;0，qgt;0 ，由 γ，τ 相互獨立可得兩個參數聯合先驗分布為正態—伽馬分布（Normal一Gamma）[21]： [γ，τ]～NG（μ_γ，V，a，q） ？ agt;0，qgt;0 。

根據正態一伽馬的概率密度函數公式，滿足：

給定未知參數 β 的先驗概率密度函數有 。

3.3后驗分布以及預測分布

輸人協變量 u=（u₁，…，u_n）^? ，得到 b（u） 是 m×n 的矩陣，令 b=b（u）^? ，對于 n 維列向量 y=（_y1，… y_n）^T ，有 y=Xβ+bγ+ε 。

在超參數 θ 和 μ_γ 已知條件下， 的條件分布滿足

[y|β，τ，γ]～N_n（Xβ+bγ，τ^-1I_n）

所有參數的聯合后驗密度函數為

同樣引入計算方法： 是 β 的最小二乘估計，整理式（18），得

對式（19）積分，得 ，且 （204 整理，得

令 P₁=I_n-X（X^TX）^-1X^T ，并將 代人式（20），得

其中 Y一μ₀^TΣ_γ|y^-1μ₀ ，式（21）中對 γ 積分，得 τ 的邊緣密度函數為：

因此， τ 的后驗分布服從伽馬分布： 。對式（21）中 τ 積分，得 γ 的邊緣密度函數： ，于是有[γ丨y]

計算 β 的后驗分布時，引入計算方法 ，其中 是 γ 的最小二乘估計 ，代入式（18），得

其中 和 SS_P 與 γ 相互獨立。 V_π=

對式（22）中 γ 積分，得

代入 到式（23）并對 τ 積分，令 H=（b^?b）^-1b^? ， A=I_n-b（b^Tb）^-1b^T ，得到 β 的后驗分布為

其中 ]-1，（2 。

3.4超參數估計

當超參數未知時，使用極大似然法得到超參數 μ_γ?，θ 的估計。由式（17）似然函數有

[y|σ²，θ，β，μ_γ]～N_n（Xβ+bμ_γ，Σ_?NGP）

其中， Σ_NGP=σ²（bVb^T+I_n） ，不考慮有限可加常數項，對數似然滿足

式（26）中超參數不具備顯式解，可以采用標準的非線性優化算法解決。本文使用L-BFGS-B算法結合R包中的optim函數[22]即可得到超參數 μ_γ?O 的估計，再結合各參數封閉形式的后驗分布可得到各個參數的貝葉斯估計。基于重構參數化的貝葉斯方法通過給定參數先驗的方式得到各參數封閉式的后驗和預測分布，節省了計算成本。特別是在數據量較大，數據維度較高的情況下，該方法通過選擇合適的 m 結點子集估計超參數可以降低計算時間。

4數值模擬

生成數據集時， y_i=x_i^Tβ+f（u_i）+ε_i，i=1，…，n，{ε_i}_i=1ⁿ 中的每個 ε_i 獨立同分布并且服從均值為0、參數 σ²=0.01 的正態分布;協變量 x_i=（x_i1，x_i2）^T，x_i1，x_i2 是獨立同分布的隨機數，均由標準正態分布N（0，1） 生成，選擇線性部分參數向量 β=（1.2，1.3）^T 。非參數協變量選擇 k=3 的三維協變量 u_i=（u_i1 ，u_i2，u_i3）^T ，其中， u_i1，u_i2，u_i3 是獨立同分布的隨機數且服從 [0，1]^k 上均勻分布。在數值實驗中，未知函數f（u_i） 分別選擇兩個測試函數： 。

比較基于高斯過程先驗（GP）的貝葉斯推斷方法以及基于插值的正態伽馬先驗（NGP）的貝葉斯推斷方法。樣本量分別取 n=200 和500，實驗重復30次，根據模擬實驗得到的部分線性模型中參數的真實值和后驗均值之間的均方根誤差（RMSE）評估兩種方法的性能，結果見表1。在測試函數相同的條件下，當樣本量從200增加到50O時，通過重構參數化方法得到的后驗樣本RMSE更小。樣本量增加到 n=1 000 時，兩種方法在不同測試函數下RMSE見表2，GP方法計算機無法得到結果，而NGP方法則可以得到非參數部分以及參數部分的后驗樣本估計;隨著樣本量的增加，NGP方法得到的后驗樣本RMSE有下降的趨勢。

表1 模型均方根誤差比較結果

表2樣本量增加后模型比較結果 （n=1000）

5 結論

對于部分線性模型中非參數分量估計問題，本文基于插值重構方法對重構后得到的新參數向量提出了正態伽馬分布先驗的方法（NGP），并據此得到了部分線性模型中非參數分量以及參數分量的估計。該方法與基于高斯過程先驗的貝葉斯推斷方法（GP）相比，可以得到參數的封閉式后驗推斷，減少了計算的復雜度。數值實驗結果表明，本文提出的NGP方法在解決部分線性模型的估計問題上是有效的，適合分析數量較大的數據，并且對復雜化的數據具有應用價值。雖然重構參數化方法在處理高維復雜數據時的可以得到比較好的結果，但是當數據集稀疏時插值方法的有效性會大大降低。在利用部分線性模型分析和處理數據集時，除了必要的數據預處理，線性和非線性協變量的選擇也是十分重要的，這一問題或許可以借助其他的方法比如似然比檢驗等解決。NGP方法也可以應用于其他復雜模型，比如廣義加性模型或者部分線性變系數模型等。后續研究可以對NGP方法的漸進性質進行了解和拓展，進一步優化和完善估計理論。

參考文獻

[1]ENGLERF，GRANGERCWJ，ICEJ，etal.SemiparametricestimatesofterelationbetweenweatherandelectricitysalesJornal of the American Statistical Association，1986，81（394）：3l0-320.

[2]HARDLE W，LIANG H，GAO J. Partially linear models[M]. Heidelberg：Springer，2000.

[3]HECKMANE.SplinesmoothinginapartlylinearmodelJ]JournaloftheRoyalStatisticalSociety：SeriesB（Methodological），986，48（2）：244-248.

[4]SPECKAPKerelsmotginpartialineaodels]oualofteoyalStatisticalSocietyeriesB（Methodological）988，50（3）：413-436.

[5]ROBINSON P M. Root-n-consistent semiparametric regression[J]. Econometrica，1988，56 （4）： 931-954.

[6]LIANH，ZHOKF，LUSGProjetedsplineestiationofteonarametricfunctioninghimesionalprtalllinearodelsformassive data[J]. The Annals of Statistics，2019，47（5）：2922-2949.

[7]ZHANGY，ANGSMonotonefunctionestimationinpartialylinearmodels]StatiticsandItsInterface，2018，）：19-9.

[8]GENGX，AHOCAOF.EstimationofapartialylinaraditivemodelwithgeneratedcovariatesJ]Joualftitcal Planningand Inference，2020，208：94-118.

[9]趙靜．基于LASSO的部分線性單指標模型局部懲罰樣條估計[J]．統計與決策，2021，37（21）：19-23.

[10]王鵬鵬，肖燕婷．縱向數據下非參數帶測量誤差的部分線性變系數模型的估計[J]．應用數學，2022，35（2）：291.