半参数顺序回归的贝叶斯推断

赵焕丽, 何幼桦

(上海大学理学院, 上海 200444)

Cox[1]提出的logit 模型是二分类问题最常用的解决方案.而对于多分类或多等级问题,McCullagh[2]扩展二分类logit 模型为比例优势模型.在比例优势模型中, 设响应变量z为具有有限个等级1<2< ··· < K的顺序变量, 在自变量x ∈Rk×1处, 记pj(x) =P(z=j|x),j=1,2,··· ,K, 其累积概率

由此, 在x处,z所属等级落入{1,2,··· ,j}或{j+1,j+2,··· ,K}的概率分别为γj(x)和1−γj(x).将{1,2,··· ,j}和{j+1,j+2,··· ,K}作为两个分类, 由logit 模型[1]得到针对顺序变量的比例优势模型

式中:αj关于等级j单调递增;β ∈Rk×1为待估计参数.式(2)是一个参数回归模型.

当响应变量与解释变量之间的关系不能用有限个参数描述时, 一般引入非参数模型.当响应变量与部分解释变量满足线性关系, 而与其他解释变量的关系不能用有限个参数描述时, 则建立半参数回归模型是合适的.熊笛等[3]用一个连续非线性函数替换式(2)中的线性部分, 建立了一种半参数顺序回归模型

并构造了参数部分αj的最小二乘估计和非参数部分f(·)的局部线性估计量.

1986 年, Engle等[4]提出了一般形式的半参数模型

式中: (x,u)为解释变量;β和f(·)为待估计参数和待估计函数.

考虑到半参数模型适用范围广泛, 本工作在模型(4)的基础上, 建立了更一般形式的半参数顺序回归模型

式中:γj(x,u) =P(z≤j|x,u)为给定解释变量x,u时z所属等级不超过j的概率.模型(5)是模型(2)与模型(3)的推广.

对于半参数模型的估计, 研究人员给出了多种经典估计方法[5-7].关于贝叶斯方法,Koop等[8]对局部线性正态回归模型的半参数推断发展了一种新的贝叶斯方法, 但该方法只讨论了取共轭先验下的正态线性模型, Koop等[9]使用了类似的方法.李琪琪[10]导出了半参数回归模型中参数的贝叶斯最小风险线性无偏估计, 并讨论了相对于最小二乘加权估计(least squares weighted estimation, LSWE)的优良性.Dimitrakopoulos[11]在误差项服从Dirichlet 过程时发展了一种估计含随机波动的时变参数回归模型的贝叶斯半参数方法, 并将其应用在通货膨胀问题的分析中.Chow等[12]在非参数部分服从Dirichlet 过程时利用数值方法计算出贝叶斯估计值.Kim等[13]构造了正态先验下参数与非参数部分的一种贝叶斯估计.

本工作针对建立的半参数顺序回归模型(5), 以非参数部分随机过程的有限维分布作为先验, 构造了参数部分与非参数部分的贝叶斯估计量.基于Kim等[13]提出的方法, 在正态情形下推导出估计量的解析表达式, 最后通过仿真模拟与实证分析评价贝叶斯估计量的表现.

1 贝叶斯估计

根据样本{(xi,ui,zi),i=1,2,··· ,n}构造模型(5)中α,β,σ2,f(·)的贝叶斯估计, 从而在新状态{(xi,ui),i=n+1,n+2,··· ,n+t}下估计zi的值.下面分3 步构造α,β,σ2,f(·)的贝叶斯估计量: ①构造pj(x,u)的贝叶斯估计; ②根据式(1)计算γj(x,u); ③结合先验分布推导α,β,σ2,f(·)的后验分布, 从而构造贝叶斯估计量.

本工作仅就正态情形给出推导过程, 但估计方法适用于一般分布.在正态情形下, 可写出估计量的解析表达式; 在非正态情形下, 可利用蒙特卡罗-马尔科夫(Monte Carol-Markov chain, MCMC)等方法进行数值计算.

当(xi,ui)处的观察值较少时,pij=pj(xi,ui)的极大似然估计值可能会取不合理的0 或1,因此首先对pij进行贝叶斯估计.取pi= (pi1,pi2,··· ,piK)的先验为Dirichlet 分布, 响应变量服从多项分布, 即

则后验

因此,pij的贝叶斯估计[14]为

在(xi,ui)下, 模型(5)满足因此, 在不计一个常数的情况下,{αj,1 ≤j

当n=2,t=1,K=3,k=1 时, 模型(5)的样本模型为

T=为n阶单位矩阵的每一行重复K−1 行,f(u)=(f(u1),f(u2),··· ,f(un+t))′,ε=(ε11,··· ,ε1,K−1,ε21,··· ,εn,K−1)′, 则

为对t个新样本预测其分类, 需估计参数β⋆和函数f(·).但f(·)是一个随机过程, 有无穷多个参数, 因此无法直接估计函数f(·), 但可估计f(·)在un+1,un+2,··· ,un+t处的函数值下面对正态情形构造参数部分与非参数部分的贝叶斯估计量.

在扰动项ε服从正态分布的情况下, 有

为了构造参数和非参数部分的贝叶斯估计量, 设参数部分β⋆的先验服从多维正态分布,函数f(·)服从高斯过程(Gaussian process, GP), 则f(·)在u处的函数值服从高斯过程的有限维分布N(·,·).σ的先验分布取Jeffreys 先验, 即

以下推导过程中将f(u),f0(u)分别简记为f,f0.

定理1记的后验分布为正态-逆伽马分布.

证明 根据贝叶斯理论, 在先验(15), (16)和(17)满足的情况下, (β⋆,f,σ)的后验分布为

式中:

式(18)的推导主要用到矩阵二次型的化简.

注意到, 对任意v ̸=0∈RK−1+k+n+t,

式(23)中的“>”号由v1,v2不同时为0 保证.因此, 矩阵A正定, 式(19)中的A−1存在.

由式(18)可知,θ,σ|Y,x,u服从正态-逆伽马分布(normal-inverse-gamma distribution,NIG), 即分别对σ2,θ积分, 可求得θ,σ2的边际分布.

推论1后验分布

在二次损失下,θ的贝叶斯估计为后验期望估计即可得^θ∗,f.

推论2后验分布σ2|Y,x,u服从逆伽马分布(inverse-gamma distribution, IGa), 即

在二次损失下,σ2的贝叶斯估计为后验期望估计

对新的样本(x,u), 根据

预测

2 数值模拟

观察在不同先验分布、样本点处的观察次数下^β⋆,f,^σ2的表现.设回归模型样本量n= 70, 预测t= 30 步, 研究不同先验分布, 不同观察次数下的表现.先验设置中, 有

表1 N =1 000 时的参数估计Table 1 Estimation of parameters when N =1 000

(2) 为探究先验分布对估计量效果的影响, 以f(u)为例, 取先验均值f0(u)=16u(0.5−u),重复观察次数m=3(随机种子为123), 由推论1 估计得到的与真实的f(u) = sin(2πu)的比较如图1(b)所示.

图1 非参数部分f 与估计的比较Fig.1 Comparing of nonparametric f and estimation

图1(a)表明, 在f(u)的先验分布的均值f0(u) = 0 时,会对先验分布有一定修正, 估计值趋向于真实的f(u).每个样本点处的重复观察次数m越多,的表现越好.表1 中与真实的α2,α3差异不大,与真实的β,σ2有一定差异, 但会随着重复观察次数m的增大越来越接近真实值.从,的1 000 次模拟的方差来看, 每个样本点处的重复观察次数m越多,,,,的估计值越稳定.

对比图1(b)所示的f0(u)=16u(0.5−u)与f(u)=sin(2πu)的曲线,在区间[0,0.7]内f0(u)接近f(u),估计得到的接近真实值;在区间[0.7,1]内f0(u)偏离f(u)很大,端点处f0(1)=−8,估计得到的很大程度上修正了先验分布与真实值之间的偏离,端点处接近−2.3.因此在先验分布接近真实值时,的估计效果较理想, 当先验分布偏离真实值时,也有不错的表现.同理, 可探究先验分布的选取对的影响.

3 实证分析

对收入等级预测问题建立一个半参数顺序回归模型, 以食品、衣着、居住、家庭设备及用品、交通通信、文教娱乐、医疗保健和其他共8 项主要生活性消费支出[15]的占比作为解释变量, 根据消费结构预测收入等级.在很多经济调查中, 真实收入的收集非常困难, 本工作根据消费习惯来预测收入等级, 有助于基于收入进行数据分析.

根据国家统计局2002—2012 年人均消费支出占比的统计数据(http://data.stats.gov.cn/easyquery.htm?cn=C01), 对不同收入等级居民的消费结构进行实证分析.8 项消费支出占比用每项的消费支出/总消费支出计算得到.农村居民家庭的收入等级采取国家统计局的五等份划分法, 分为低收入、中等偏下收入、中等收入、中等偏上收入、高收入家庭(分别用1, 2, 3, 4,5 表示); 城镇居民家庭的收入等级在国家统计局的划分法基础上稍作改动, 将低收入与较低收入家庭划分为低收入家庭, 较高收入与最高收入家庭划分为高收入家庭, 这样城镇居民家庭由原来的7 个收入等级改为5 个收入等级.

选取2002—2011 年的数据(共120 组)作为训练样本, 用半参数顺序回归模型(5)拟合8 项人均消费支出占比与家庭收入等级之间的关系.令u为8 项消费支出占比,

则由模型

从表2 的实验结果可以看到, 对收入五等级问题的12 组样本预测中, 预测准确率为58.33%.预测错误的5 组样本中, 有4 组样本的预测等级与实际等级只相差1 个等级.

表2 半参数顺序回归模型的外推收入等级Table 2 Extrapolation income level of semiparanetric ordinal regression mode

已有研究表明, 在诸多影响消费需求的因素中, 收入水平始终是影响消费需求的最重要因素.而本工作实证分析的结果表明, 以居民家庭的消费结构为解释变量可以较准确地预测收入等级, 因此消费结构反过来也反映了家庭的收入情况.以8 项生活性消费支出考虑家庭的收入等级的方式, 比采用恩格尔系数更为全面.综上可知, 在小样本情况下, 贝叶斯估计由于利用了先验信息, 往往有更好的估计效果.

4 结束语

本工作在正态情形下构造了半参数回归模型中参数与非参数部分的贝叶斯估计, 多次模拟结果表明, 在先验均值均取0 时仍然有不错的估计效果, 在先验分布接近真实值时, 估计效果会更理想.

相比模型(3), 本工作建立的模型是基于比例优势模型和半参数模型(4)所形成的更一般的半参数顺序回归模型.模型(3)只考虑了被解释变量与非参数部分变量u的关系, 无法反映参数部分变量x的影响, 而本模型建立了解释变量与被解释变量的半参数关系, 适用范围更广泛.例如, 当x为外生变量时, 本模型可反映不同情况下u对被解释变量的影响.在根据消费结构预测家庭收入等级的实例中, 本工作将线性部分取为哑变量, 使模型对城镇居民和农村居民的收入等级预测问题都适用, 同时线性部分的参数反映了相同收入等级的城镇居民家庭与农村居民家庭的消费结构差异.以2002—2011 年间的数据作为样本, 对2012 年的数据作预测,即在外推情形下, 对五等级预测问题的预测准确率达到58.33%, 预测错误的样本中, 预测等级与实际等级大多只相差1 个等级.