基于偏正态数据下位置、均值回归模型的参数估计

王丹璐，吴刘仓，郑桂芬

(昆明理工大学理学院，云南昆明650504)

1.引言

生活中常见的统计分布可分为两类:对称分布和偏态分布.我们所学的大部分分布如正态分布、t分布、Laplace分布等都属于对称分布，故对称分布在统计理论与方法的研究中占有重要地位.但是现实生活中收集到的大部分数据都不具有严格的对称性，几乎都存在偏斜，这时再用对称分布描述数据的性质并不科学，所以引入了偏态分布来解决这一问题.偏态分布具有偏斜、尖峰、厚尾的特征，由于偏态分布比对称分布更能捕捉到全面准确、及时有效的信息，所以近年来偏态分布引起了许多学者的研究兴趣，如郑等[1]研究了基于偏Laplace正态数据下的位置、均值回归模型的参数估计，朱等[2]研究了偏t正态数据下混合线性联合位置与尺度模型的参数估计.此外，由于偏态分布在数据拟合运用方面更广阔，所以在农业、环境、生物医学应用及金融经济等领域也有很多研究意义.

偏正态分布作为正态分布的推广，不仅具有正态分布的良好统计特性，同时也有具有偏态分布的特征，适用性更广泛，所以国内外许多学者研究了偏正态分布问题.自Azzalini[3]首次定义了偏正态分布的概念，随后关于它的性质、参数估计、统计诊断、推广以及多元偏正态分布已经被系统深入地研究.本文中我们使用的是Sahu等[4]提出的多元偏正态分布的单变量形式，该偏态分布是基于偏椭圆分布变换而来的，相比Azzalini提出的偏正态分布该函数更具有吸引力，因为该分布对于偏度参数的刻画较简单，使用EM算法更容易实现.在Sahu提出的偏正态分布模型中已有部分研究，比如Cancho等[5]在经典的贝叶斯方法的基础上，对偏正态非线性回归模型进行了深入的探讨;Bolfarine等[6]提出了多元偏正态线性模型的四种诊断措施;Lachos等[7]在多元偏正态回归模型基于极大似然估计的EM算法，研究了多元偏正态回归模型的参数估计;Arellano-Valle等[8]研究了偏正态混合线性模型;Arellano-Valle等[9]研究了其贝叶斯推断.

目前的对于该偏正态分布中，并没有对位置和均值建模的内容，但考虑到位置和均值在偏态分布中的重要地位，所以本文详细介绍了利用EM算法和梯度下降法对位置回归模型和均值回归模型的参数估计，并通过模拟和实例结果表明本文提出的模型是可行的.

本文的目的是讨论偏正态位置、均值回归模型的参数估计，全文结构安排如下:第二部分介绍了偏正态分布的部分性质;第三部分分别介绍了偏正态分布下位置和均值回归模型;第四部分利用梯度下降法和EM算法分别对位置和均值回归模型的参数进行了估计;第五部分通过Monte Carlo模拟证实本文提出方法的有效性;第六部分通过实例分析并且和多元线性回归模型比较，表明本文研究的模型和方法是科学合理的.

2.偏正态分布

根据Sahu(2003)等提出偏正态函数，若随机变量Y服从偏正态分布，那么它的概率密度函数(pdf)记做:

其中φ(·|a，b2)和Φ(·|a，b2)分别表示N(a，b2)的密度函数(pdf)和分布函数(cdf)，记做Y～SN(μ，σ2，λ)，其中μ为位置参数，σ2为尺度参数，λ为偏度参数.当λ＜0时，曲线较长的尾部在左边，称作负偏分布;当λ＞0时，曲线较长的尾部在右边，称作正偏分布;当λ=0时，偏正态分布重新退化为对称的正态分布.

I偏正态分布下的随机表示

设X0～N(0，1)，X1～N(μ，σ2)，是两个独立随机变量，则随机变量Y～SN(μ，σ2，λ)的随机表示为

从而得到偏正态分布的层次表示

即该偏正态分布可以分层为一个条件正态分布Y|(T=t)和一个半正态分布因此(Y，T)的联合密度函数为

根据贝叶斯准则有

3.位置和均值回归模型

I位置回归模型

由偏正态分布的概率密度函数(2.1)和位置回归参数模型(3.1)可以得到

其中yi为第i个响应变量，位置参数为μi，尺度参数为σ2，偏度参数为λ，xi=(xi1，xi2，···，xip)T为解释变量，β=(β1，β2，···，βp)T为p×1维的位置回归模型的未知参数.

II均值回归模型

其中α=(α1，α2，···，αp)T为p×1维的均值回归模型的未知参数.

因为有潜变量的存在，所以利用普通极大似然估计方法对参数估计比较困难，故我们使用基于梯度下降法的EM算法来解决含有潜变量问题的参数估计.接下来研究模型参数的极大似然估计的EM算法.

4.极大似然估计的EM算法和梯度下降法

自Dempster等[10]提出EM算法(Expectation Maximization Algorithm)来处理不完全数据以来，EM算法及其扩展应用已经被广泛的研究和应用到各个领域.EM算法是一种常用的迭代算法，每一次迭代由两部分组成:E-step，M-step.E-step是根据参数初始值或上一次迭代所得到的结果来计算出对数似然函数的期望值;M-step是通过E-step计算得到的对数似然函数的期望值，求出该似然函数的最大的参数值，记为新的参数值，用该参数值代替初始值或上一次迭代所得到的结果使得对数似然函数最大化.E-step和M-step两步不断交替计算，直到满足收敛条件时停止迭代，从而得到的参数估计值更加的逼近真实的参数值.

下面给出EM算法在偏正态数据下位置回归模型和均值回归模型的参数估计中的具体计算步骤:

I位置回归模型下极大似然估计的EM算法

由(3.2)式可得偏正态数据下位置回归模型的似然函数为

令θ1=(βT，σ2，λ)T，则L(βT，σ2，λ)=L(θ1)，两边取自然对数，得到对数似然函数

设T为潜变量，利用位置回归模型和偏正态分布的层次表达法得

其中c1为独立于θ1的常数.同时，因为该似然函数是单峰可导的，则似然方程的解存在唯一，且为极大似然估计，因而是强相合的.为了得到θ1的极大似然估计，将(4.3)式极大化，但此时极大化得到的估计依赖于潜变量.因此，为了处理潜变量，引入EM算法，使用观测数据yi求出完全数据下对数似然函数的条件期望.

E-step:利用k次迭代，得到估计

利用(2.5)和(2.6)式给出的条件期望，可以计算出(4.4)式中的条件期望E和有

M-step:两步长梯度法在构造搜索方向时，充分利用上次迭代点信息，使计算不易陷入局部最大值的陷阱，故结合梯度下降法[11]将对数似然函数极大化，获取新的参数值，给出两个初值设计以下迭代

其中

利用(4.7)式给出的目标函数，得到位置参数的Score函数为

利用梯度下降法将E-step和M-step反复迭代，直至算法收敛.

II均值回归模型下极大似然估计的EM算法

令θ2=(αT，σ2，λ)T，则L(αT，σ2，λ)=L(θ2)，两边取自然对数，得到对数似然函数

利用均值回归模型和偏正态分布的层次表达法得

用层次表达法得到完全数据下的对数似然函数为

其中c2为独立于θ2的常数.

E-step:利用k次迭代，得到估计:

利用(2.5)和(2.6)式给出的条件期望，可以计算出(4.10)的条件期望和有

皮肤表面脂质是由甘油三酯、蜡酯、角鲨烯、脂肪酸以及少量胆固醇、胆固醇酯和双甘酯组成的非极性脂类混合物[14]，其中角鲨烯具有高度的不饱和性，在UV诱导下，角鲨烯脂质过氧化形成角鲨烯过氧化物（SQOOH）以及丙二醛（MDA）[15]，增加羰基化反应、糖基化反应，从而导致皮肤泛黄。

M-step:利用梯度下降法将对数似然函数极大化，获取新的参数值，给出两个初值设计以下迭代

其中

利用(4.14)式给出的目标函数，得到均值参数的Score函数为

利用梯度下降法将E-step和M-step反复迭代，直至算法收敛.

III极大似然估计的渐近性质

同时考虑该极大似然估计的渐近正态性，取θ0为θ的真值，假定以下正则条件

(i)xi=(xi1，xi2，···，xip)T是固定的;

(ii)参数空间是紧的，真值θ0为参数空间的内点;

(iii)yi，i=1，2，···，n相互独立，Y～SN(μ，σ2，λ)，其中在位置回归下的极大似然估计时，;在均值回归下的极大似然估计时，.

在假定(i)-(iii)之下，似然方程L(θ)=0在L(1)(θ)=0在n→+∞时有相合解满足且似然方程任一相合解必为θn的BAN估计.

5.Monte Carlo模拟

I位置回归模型下的Monte Carlo模拟

为评价上述位置回归模型参数估计方法的有效性，使用该位置回归模型参数估计方法对有限样本进行Monte Carlo模拟，并用均方误差(MSE)、偏差(Bias)来评价该参数估计的精确度，其定义如下

其中θ1=(βT，σ2，λ)，θ1(0)是参数θ1的真值.根据模型(3.1)，考虑偏正态数据下位置参数的回归模型.

根据模型(5.1)产生模拟数据，其中解释变量x的分量独立产生于均匀分布U(-1，1)，其余参数的真值设定为-2.

II均值回归模型下的Monte Carlo模拟

为评价上述均值回归模型参数估计方法的有效性，使用该位置均值模型参数估计方法对有限样本进行Monte Carlo模拟，并用均方误差(MSE)、偏差(Bias)来评价该参数估计的精确度，其定义如下

其中θ2=(αT，σ2，λ)，θ2(0)是参数θ2的真值.根据模型(3.3)，考虑偏正态数据下均值参数的回归模型.

分别取样本量n=100，200，300，400，重复模拟1000次.模拟结果见表1，评价结果见表2、表3.

表1 位置、均值回归模型的参数估计模拟结果

表2 位置回归模型的参数估计评价结果

表3 均值回归模型的参数估计评价结果

从上表可以得到，无论样本是左偏还是右偏，即参数λ是负值还是正值，随着样本量n的增大，所有参数的估计值越来越接近真值，而且通过估计量的均方误差(MSE)和偏差(Bias)的判断也可以得到当样本量不断增大，估计效果越好.以上结论表明，本文提出的偏正态数据下位置和均值回归模型及所使用的EM算法对参数的极大似然估计取得了较理想的效果.

6.实例分析

根据2017年国际糖尿病联名协会公布的一组数据，目前全球有4.7亿糖尿病患者，其中中国糖尿病患者高达1.14亿，位居世界第一，所以对糖尿病患者数据的研究是必要的.本文选取UCI机器学习数据库的一组糖尿病数据进行分析，数据包括了478名糖尿病患者的身体质量指数BMI和其他相关指标，该数据中包含一个响应变量Y-糖尿病患者BMI，7个解释变量X，分别是怀孕次数、血糖(mmol/L)、血压(mmHg)、皮质厚度(mm)、胰岛素(pmol/ml)、糖尿病遗传函数、年龄.利用R软件计算得y偏度为0.5529，峰度为3.0376，从峰度值和偏度值，可以知道，糖尿病患者的BMI数据具有明显的偏斜，而正态分布的偏度值为0，峰度值为3，直方图如下:

图1 糖尿病患者BMI的直方图

根据图1和偏度系数说明糖尿病患者的身体质量指数BMI数据近似的服从偏正态分布，所以可以利用该组数据做偏正态分布的位置和均值回归模型的参数估计，考虑Y与X之间的模型如下:

同时，我们加入多元线性回归模型做比较，多元线性回归模型如下:

糖尿病患者BMI的位置、均值回归模型参数估计和多元线性回归参数估计结果:

表4 糖尿病患者BMI的位置、均值回归模型参数估计和多元线性回归参数估计结果

由于在同一组糖尿病患者BMI的数据中，尺度参数和偏度参数应为一致，从表格中分析到，位置回归模型和均值回归模型的尺度参数σ2和偏度参数λ基本一致，但是多元线性回归模型的尺度参数σ2与其他两个模型差别较大，并且不存在偏度参数λ.这是由于使用多元线性回归模型需要忽略数据的偏斜特征.从解释变量的系数来看，多元线性回归模型中的怀孕次数γ1、血糖γ2、血压γ3、皮脂厚度γ4、年龄γ7都不显著，并且系数值与其余两个模型相差过大，胰岛素γ5呈现正相关，这都与实际情况不符.以上，我们可以判断出该数据用多元线性回归模型做参数估计导致估计效果不合理甚至错误.这是由于这组数据中并不是严格对称的，具有一定的偏斜，但多元线性回归模型针对的是对称分布数据进行统计推.所以在数据呈现偏正态情形下，不应该再选用多元线性回归模型做参数估计，可以选择本文提出的偏正态数据下的位置回归模型和均值回归模型做数据分析.

在均值和位置回归模型中，解释变量中血糖X2、皮脂厚度X4、糖尿病遗传函数X6呈现正相关且系数较大，这与实际情况相符合.血糖升高是判断是否患糖尿病重要指标，肥胖是二型糖尿病的主要诱因，并且如果父母患糖尿病那么子女患糖尿病的风险比正常人要高.解释变量中胰岛素X5呈负相关，这与患有糖尿病的患者缺少胰岛素实际情况相符合.

7.结论

Monte Carlo模拟表明了本文提出的EM算法对位置和均值回归模型的未知参数进行了较好的估计，与现有的模型和常用估计方法相比较，本文提出的偏正态数据下位置和均值回归模型具有更好的运用性，更适合实际生活中的偏态数据的研究.除此以外，在实例分析中，对糖尿病患者BMI数据的应用研究也表明了本文提出的模型和方法能更准确分析解释数据，因此表明本文提出的模型和方法是可行的.