郭媛媛¹, 杨雪梅², 孙志华^1,3
1. 中国科学院大学数学科学学院, 北京 100049;
2. 华北电力大学数理学院, 北京 102206;
3. 中国科学院大数据挖掘与知识管理重点实验室, 北京 100190
2019年5月14日 收稿; 2019年10月9日 收修改稿
基金项目: 国家自然科学基金（11971045）、中央高校基本科研业务费专项资金和中国科学院大数据挖掘与知识管理重点实验室开放课题资助
通信作者: 孙志华, E-mail: sunzh@amss.ac.cn

摘要: 单指标分位回归模型是一类重要的半参数模型，具有降维的优点的同时保留了非参数分位回归模型的稳健性。但现有的单指标分位回归模型的估计程序大部分都是通过内点法来实现。对单指标分位回归模型估计程序的MM（majorize-minimize）算法进行研究。首先找到目标函数的优化函数，然后通过最小化优化函数来得到估计，再逐步迭代至收敛。数值模拟和实证研究表明MM算法在单指标分位回归模型的估计中具有较好的稳定性，能够得到比较准确的估计结果，且相比于内点法，计算效率更高，耗时更短。
关键词: 单指标分位回归模型MM算法替代函数计算效率
MM algorithm of the estimation of single-index quantile regression
GUO Yuanyuan¹, YANG Xuemei², SUN Zhihua^1,3
1. School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China;
2. School of Mathematics and Physics, North China Electric Power University, Beijing 102206, China;
3. Key Laboratory of Big Mining and Knowledge Management, Chinese Academy of Sciences, Beijing 100190, China


Abstract: The single-index quantile regression model is an important semiparametric model with the merit of dimensionality reduction. Furthermore, it retains the robustness of a nonparametric model. For most existing estimating procedures of single-index quantile regression models, the estimators are obtained via minimizing the objective functions by the interior point method. In this paper, we investigate the MM (majorize-minimize) algorithm of the single index quantile regression model estimating procedure. We first construct the majorize function of the objective function and then minimize the substituted majorize function to find the estimators. Our numerical simulations and empirical study show that for the considered model, the MM algorithm has good stability and can yield more accurate estimation. Compared with the interior point method, the MM algorithm is more efficient and takes less time.
Keywords: single-index quantile regression modelMM-algorithmsurrogate functioncomputational efficiency
分位回归模型具有稳健的特点，并且能够对响应变量的分布做出精细的描述，因此获得很多****的关注。目前，分位回归方法已经成为分析数据的一个非常重要的工具，广泛地应用在金融、医疗、生物研究等领域。分位回归模型的研究和应用可参考文献[1-4]及相关文献。参数分位回归模型是在实际中应用广泛的一类回归模型^[5]。但很多时候无法对参数分位回归模型进行正确的设定。基于误定的参数分位回归的统计推断结果经常是不可信的，文献[6-7] 阐述了参数分位回归需要进行检验的必要性。非参数分位回归不存在误定的风险，但当样本量比较小且协变量比较多时，非参数分位回归方法可能会受到维数祸根的问题的困扰。
对分位回归，构建目标函数时用到的损失函数ρ_τ(r)=τr-I_{{r < 0}}具有不光滑的特点, 从而使得求解目标函数的最小值比较困难，且可能出现多个最小值点的情况，参见文献[5, 8-11]。一种解决上面问题的方法是将分位回归模型的求解问题转化为线性规划问题，再利用单纯形法或内点法进行计算。不管是单纯形法，还是内点法，运算效率都不能令人满意。2000年，Hunte和Lange^[12]提出一种新的用于求解分位回归问题的算法，即MM算法。MM算法概念简单，易于执行，且数值稳定，比内点法拥有更强的数值计算能力。文献[11]对4种求解分位回归问题的算法，即内点法、MM算法、坐标下降法和ADMM算法进行比较研究，验证了MM算法具有数值稳定和计算效率高的特点。
单指标分位回归模型具有降维的效果，同时保持了非参数分位回归的稳健性，其估计问题的研究吸引了很多研究者的兴趣。文献[13-14] 提出基于两步迭代的估计方法，文献[15-16]进一步提出不需要迭代的估计方法，文献[17-18]又提出基于贝叶斯方法的估计方法，文献[19-21]探讨单指标分位回归模型的变量选择以及加权复合单指标分位回归模型的估计。然而这些文献所提的估计方法都基于内点法来实现，内点法在计算分位回归模型时，计算效率低、耗时久，尤其在样本量较大的情况下，这种缺点更为明显。MM算法在求解分位回归模型的估计时比较高效和便捷，这在文献[22-24]中均有体现，但是没有文献研究单指标分位回归模型的MM算法，故本文研究单指标分位回归模型估计的MM算法。
我们借鉴文献[12] 的方法，对单指标分位回归模型的每一步迭代程序中目标函数构建其替代函数，从而将复杂的优化问题简单化。然后，基于优化函数再进行求解计算得到估计值。我们构建的优化函数是光滑的，并能够保证每次迭代目标函数是下降的。数值模拟和实例分析结果表明基于MM算法的估计程序具有较好的稳定性，能够得到比较准确的估计结果，并且相较于传统的内点算法具有更强的数值计算能力，用时更短。
1 单指标分位回归模型的估计介绍对于给定的分位数τ∈(0, 1)，在给定x的条件下，响应变量y的τ分位数θ_τ(x)与协变量x之间的关系如下：

$\theta_{\tau}(\boldsymbol{x})=g\left(\boldsymbol{x}^{\mathrm{T}} \boldsymbol{\gamma}\right), $

其中x∈?^d是d维协变量，g(·)表示未知的一元联系函数。另外γ=(γ₁, …, γ_d)^T为未知的单指标向量，满足‖γ‖=1且γ₁>0, ‖·‖表示Euclidean范数。这个约束条件是为了模型的可识别性^[25]，已广泛应用在有关单指标模型的文献中。
本文采用局部线性方法对γ和g(·)进行估计，详细内容可参考文献[13]，具体算法如下：
step 1??采用文献[26]提出的平均导数法得到γ的初值，记为${\mathit{\boldsymbol{\hat \gamma }}^{\left( 0 \right)}}$, 并对其进行标准化，使其满足$\left\| {\mathit{\boldsymbol{\hat \gamma }}} \right\| = 1$且${\mathit{\boldsymbol{\hat \gamma }}_1} < 0$。
step 2??对给定的$\mathit{\boldsymbol{\hat \gamma }}$，考虑下面的最小化问题：

$\min \limits_{\left(a_{j}, b_{j}\right)} \sum\limits_{i=1}^{n} \rho_{\tau}\left(y_{i}-a_{j}-b_{j}\left(\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right)^{\mathrm{T}} \hat{\boldsymbol{\gamma}}\right) w_{i j}, $

(1)

其中ρ_τ(r)=r(τ-I(r≤0))，w_ij=$\frac{{{K_h}\left( {\mathit{\boldsymbol{x}}_i^{\rm{T}}\mathit{\boldsymbol{\gamma }} - \mathit{\boldsymbol{x}}_j^{\rm{T}}\mathit{\boldsymbol{\gamma }}} \right)}}{{\sum_{l = 1}^n {{K_h}} \left( {\mathit{\boldsymbol{x}}_l^{\rm{T}}\mathit{\boldsymbol{\gamma }} - \mathit{\boldsymbol{x}}_j^{\rm{T}}\mathit{\boldsymbol{\gamma }}} \right)}}$，K(·)为核函数且K_h(·)=$K\left( {\frac{ \cdot }{h}} \right)/h$。记所得估计为$\left\{ {{{\hat a}_j}, {{\hat b}_j}} \right\}_{j = 1}^n$。带宽的选择为h_τ=h_m{τ(1-τ)/?(Φ^-1(τ))²}^1/5，其中?(·), Φ(·) 分别为标准正态分布的概率密度函数和分布函数，h_m=${\left\{ {\frac{{\left[ {\int {{K^2}} (v){\rm{d}}v} \right]\left[ {{\mathop{\rm var}} \left( {y\mid {\mathit{\boldsymbol{x}}^{\rm{T}}}\mathit{\boldsymbol{\gamma }} = u} \right)} \right]}}{{n{{\left[ {\int {{v^2}} K(v){\rm{d}}v} \right]}^2}{{\left[ {\frac{{{{\rm{d}}^2}}}{{{\rm{d}}{u^2}}}E\left( {y\mid {\mathit{\boldsymbol{x}}^{\rm{T}}}\mathit{\boldsymbol{\gamma }} = u} \right)} \right]}^2}\left[ {{f_{{U_0}}}(u)} \right]}}} \right\}^{1/5}}$，为均值回归时的最优带宽，详细推导过程见文献[27]。
step 3??基于step 1得到的估计$\left\{ {{{\hat a}_j}, {{\hat b}_j}} \right\}_{j = 1}^n$，进一步考虑下面的问题：

$\min \limits_{\boldsymbol{\gamma}} \sum\limits_{j=1}^{n} \sum\limits_{i=1}^{n} \rho_{\tau}\left(y_{i}-\hat{a}_{j}-\hat{b}_{j}\left(\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right)^{\mathrm{T}} \boldsymbol{\gamma}\right) w_{i j}, $

(2)

其中w_ij和h_τ直接使用上一步计算的结果。记所得估计为${\mathit{\boldsymbol{\tilde \gamma }}}$。
step 4??重复step 1和step 2直至收敛。
在上述算法中，${\mathit{\boldsymbol{\hat \gamma }}}$标准化的方法为：γ=sign₁${\mathit{\boldsymbol{\tilde \gamma }}}$/‖${\mathit{\boldsymbol{\tilde \gamma }}}$‖，其中sign₁指γ的第一个分量的符号函数，只有第一步${\mathit{\boldsymbol{\hat \gamma }}}$需要标准化。最终g(·)在u处的估计结果为?(·;h; ${\mathit{\boldsymbol{\hat \gamma }}}$)=a，其中

$\begin{aligned}(\hat{a}, \hat{b})=\arg \mathop {\min }\limits_{a, b} & \sum\limits_{i=1}^{n} \rho_{\tau}\left(y_{i}-a-b\left(\boldsymbol{x}_{i}^{\mathrm{T}} \hat{\boldsymbol{\gamma}}-u\right)\right) \\& K_{h}\left(\boldsymbol{x}_{i}^{\mathrm{T}} \hat{\boldsymbol{\gamma}}-u\right).\end{aligned}$

(3)

2 估计程序的MM算法2.1 MM算法介绍下面介绍MM算法的基本思想。假设需要最小化的目标函数为L(θ): ?^p→?，θ^k为第k步的迭代值。MM算法每次迭代分两步来进行。首先，构造目标函数的优化函数Q(θ|θ^k): ?^p×?^p→? 满足

$\begin{array}{c}Q\left(\boldsymbol{\theta}^{k} \mid \boldsymbol{\theta}^{k}\right)=L\left(\boldsymbol{\theta}^{k}\right) , \\Q\left(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{k}\right) \geqslant L(\boldsymbol{\theta}) \forall \boldsymbol{\theta}.\end{array}$

(4)

然后, 对优化函数Q(θ|θ^k)进行最小化, 得到下一步的迭代值θ^k+1，则有

$Q\left(\boldsymbol{\theta}^{k+1} \mid \boldsymbol{\theta}^{k}\right) \leqslant Q\left(\boldsymbol{\theta}^{k} \mid \boldsymbol{\theta}^{k}\right) .$

(5)

综合式(4)和式(5)，可知L(θ^k+1)≤L(θ^k)。这种下降趋势保证了MM算法具有显著的数值稳定性。
观察目标函数式(1)~式(3)，可以发现它们均为非光滑函数，因而不易得到最优解，故借用文献[12] 提出的MM算法的思想来处理这个问题。其主要的处理方式如下：首先给ρ_τ(r)加一扰动ε，得到其近似函数

$\rho_{\tau}^{\varepsilon}(r)=\rho_{\tau}(r)-\frac{\varepsilon}{2} \ln (\varepsilon+|r|).$

然后, 在r^k处用下式来优化ρ_τ^ε(r)

$\zeta_{\tau}^{\varepsilon}\left(r \mid r^{k}\right)=\frac{1}{4}\left[\frac{r^{2}}{\varepsilon+\left|r^{k}\right|}+(4 \tau-2) r+c\right], $

其中，c是满足ζ_τ^ε(r^k|r^k)=ρ_τ^ε(r^k)的常数。下面给出单指标分位回归模型估计程序中最小化问题(1)、(2)和(3)的基于MM算法的详细求解过程。
2.2 非参数部分的估计首先，将式(1)中的目标函数修正为

$L_{\varepsilon}^{j}(\boldsymbol{\theta})=\sum\limits_{i=1}^{n} \rho_{\tau}^{\varepsilon}\left(r_{i j}(\boldsymbol{\theta})\right) w_{i j}, j=1, 2, \cdots, n, $

其中：r_ij(θ)=y_i-a_j-b_j(x_i-x_j)^T${\mathit{\boldsymbol{\hat \gamma }}}$，θ为参数(a_j, b_j)。进一步构造近似函数L_ε^j(θ)在第k次迭代值θ^k处的优化函数：

$Q_{\tau, j}^{\varepsilon}\left(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{k}\right)=\sum\limits_{i=1}^{n} \zeta_{\tau}^{\varepsilon}\left(r_{i j} \mid r_{i j}^{k}\right) w_{i j}, j=1, 2, \cdots, n.$

(6)

下面考虑第k+1次迭代值的确定问题。MM算法是通过优化函数Q_τ^ε(θ|θ^k)的最小化来得到θ的估计：使Q_τ^ε(θ|θ^k)达到最小的θ值即为下一次迭代值θ^k+1。由于残差项r为θ的线性函数，对式(6)求导，令其为零并化简整理可得

$\begin{array}{l}\sum\limits_{i=1}^{n} w_{i j}\left(\frac{Y_{i}}{\left(\varepsilon+\left|r_{i j}^{k}\right|\right)}+(2 \tau-1)\right) \boldsymbol{X}_{i}= \\\sum\limits_{i=1}^{n} \frac{w_{i j}}{\left(\varepsilon+\left|r_{i j}^{k}\right|\right)} \boldsymbol{X}_{i} \boldsymbol{X}_{i}^{\mathrm{T}} \boldsymbol{\theta}, j=1, 2, \cdots, n, \end{array}$

其中Y_i=y_i, X_i=(1, (x_i-x_j)^T${\mathit{\boldsymbol{\hat \gamma }}}$)^T, r_ij^k=y_i-a_j^k-b_j^k(x_i-x_j)^T=Y_i-X_i^Tθ^k，故θ^k+1的表达式：

$\begin{array}{c}\boldsymbol{\theta}^{k+1}=\left(\sum\limits_{i=1}^{n} \frac{w_{i j}}{\left(\varepsilon+\left|r_{i j}^{k}\right|\right)} \boldsymbol{X}_{i} \boldsymbol{X}_{i}^{\mathrm{T}}\right)^{-1} \sum\limits_{i=1}^{n} w_{i j} \\\left(\frac{Y_{i}}{\left(\varepsilon+\left|r_{i j}^{k}\right|\right)}+(2 \tau-1) \boldsymbol{X}_{i}\right).\end{array}$

(7)

再通过对式(7)的迭代，逐渐逼近目标函数L_ε^j(θ)的最小值。
由此，可以将MM算法总结为如下步骤：
1) 选择迭代初始值θ⁰和一个较小的正常数ε，置k=0；
2) 计算r_ij^k，根据式(7)确定第k+1次迭代值θ^k+1；
3) 令k=k+1，判断是否满足收敛准则，若满足收敛准则，即：当

$\left|\frac{Q_{\tau}^{\varepsilon}\left(\boldsymbol{\theta}^{k} \mid \boldsymbol{\theta}^{k}\right)-Q_{\tau}^{\varepsilon}\left(\boldsymbol{\theta}^{k+1} \mid \boldsymbol{\theta}^{k}\right)}{Q_{\tau}^{\varepsilon}\left(\boldsymbol{\theta}^{k} \mid \boldsymbol{\theta}^{k}\right)}\right|<\delta $

$\text { 或者 } Q_{\tau}^{\varepsilon}\left(\boldsymbol{\theta}^{k} \mid \boldsymbol{\theta}^{k}\right)-Q_{\tau}^{\varepsilon}\left(\boldsymbol{\theta}^{k+1} \mid \boldsymbol{\theta}^{k}\right)<\delta$

时，可终止迭代，其中δ是预先取定的足够小的数。否则返回2)继续迭代，直到满足收敛准则。
2.3 单指标向量的估计首先定义式(2)中目标函数的近似函数：

$L_{\tau}^{\varepsilon}(\boldsymbol{\gamma})=\sum\limits_{j=1}^{n} \sum\limits_{i=1}^{n} \boldsymbol{\rho}_{\tau}^{\varepsilon}\left(y_{i}-\hat{a}_{j}-\hat{b}_{j}\left(\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right)^{\mathrm{T}} \boldsymbol{\gamma}\right) w_{i j}.$

在γ^k处的优化函数可以构建为

$\begin{array}{c}R_{\tau}^{\varepsilon}\left(\boldsymbol{\gamma} \mid \boldsymbol{\gamma}^{k}\right)=\sum\limits_{j=1}^{n} \sum\limits_{i=1}^{n} \frac{1}{4}\left[\frac{r_{i j}(\boldsymbol{\gamma})^{2}}{\left.\varepsilon+\mid r_{i j}^{k}\right) \mid}+\right. \\\left.(4 \tau-2) r_{i j}(\boldsymbol{\gamma})+c\right] w_{i j}, \end{array}$

其中: r_ij(γ)=y_i-${{\hat a}_j} - {{\hat b}_j}$(x_i-x_j)^Tγ，c是满足使R_τ^ε(γ^k|γ^k)=L_τ^ε(γ^k)的常数。然后，对R_τ^ε(γ|γ^k)关于γ求导并令其为0，得到γ^k+1为

$\begin{array}{c}\boldsymbol{\gamma}^{k+1}=\left(\sum\limits_{j=1}^{n} \sum\limits_{i=1}^{n} \frac{w_{i j}}{\left(\varepsilon+\left|r_{i j}^{k}\right|\right)} \boldsymbol{X}_{i j} \boldsymbol{X}_{i j}^{\mathrm{T}}\right)^{-1} \\\sum\limits_{j=1}^{n} \sum\limits_{i=1}^{n} w_{i j}\left(\frac{Y_{i j}}{\left(\varepsilon+\left|r_{i j}^{k}\right|\right)}+(2 \tau-1)\right) \boldsymbol{X}_{i j} .\end{array}$

(8)

其中X_ij=${{\hat b}_j}$(x_i-x_j), Y_ij=y_i-a_j, r_ij^k=y_i-a_j^k-$\hat b_j^k$(x_i-x_j)^Tγ^k。按照2.2节所给出的MM算法将θ替换为γ，Q_τ^ε(·)替换为R_τ^ε(·)，对式(8) 的迭代，即可得到γ的估计。
最后，可将基于MM算法的单指标模型的估计总结为如下步骤：
1) 参考第1节step 1所提供的方法，得到γ的初始估计；
2) 对给定的${\mathit{\boldsymbol{\hat \gamma }}}$，按照2.2节计算出$\left\{ {{{\hat a}_j}, {{\hat b}_j}} \right\}_{j = 1}^n$；
3) 对于给定的$\left\{ {{{\hat a}_j}, {{\hat b}_j}} \right\}_{j = 1}^n$，按照2.3节计算出${\mathit{\boldsymbol{\hat \gamma }}}$；
4) 重复2)、3)步骤，直至收敛。
2.4 联系函数g(·)的估计对于联系函数的估计，式(6)中目标函数式的近似函数可定义为

$L_{1, \tau}^{\varepsilon}(a, b)=\sum\limits_{i=1}^{n} \rho_{\tau}^{\varepsilon}\left(y_{i}-a-b\left(\boldsymbol{x}_{i}^{\mathrm{T}} \hat{\boldsymbol{\gamma}}-u\right)\right) K_{h}\left(\boldsymbol{x}_{i}^{\mathrm{T}} \hat{\boldsymbol{\gamma}}-u\right).$

按照2.2节方法构建L_{1, τ}^ε(a, b)的优化函数，求解可得

$(\hat{a}, \hat{b})=\left(\sum\limits_{i=1}^{n} \frac{w_{i}}{\left(\varepsilon+\left|r_{i}^{k}\right|\right)} \boldsymbol{X}_{i} \boldsymbol{X}_{i}^{\mathrm{T}}\right)^{-1}$

$\sum\limits_{i=1}^{n} w_{i}\left(\frac{Y_{i}}{\left(\varepsilon+\left|r_{i}^{k}\right|\right)}+(2 \tau-1)\right) \boldsymbol{X}_{i} .$

则?(u; h, ${\mathit{\boldsymbol{\hat \gamma }}}$)=a, 其中, w_i=K_h(x_i^T${\mathit{\boldsymbol{\hat \gamma }}}$-u), X_i=(1, x_i^T${\mathit{\boldsymbol{\hat \gamma }}}$-u)^T, Y_i=y_i, r_i^k=y_i-a^k-${{\hat b}^k}$(x_i^T${\mathit{\boldsymbol{\hat \gamma }}}$-u)。
3 数值模拟对于近似函数ρ_τ^ε(r)中的扰动ε，按照文献[12]给出的建议，选取ε满足εn|lnε|=δ。
3.1 模拟1借鉴文献[13]模拟1的模型设置，考虑模型

$y = \sin \left( {\frac{{{\rm{ \mathit{ π} }}(u - A)}}{{C - A}}} \right) + 0.1Z{\rm{ , }}$

其中u=x^Tγ，x=(x₁, x₂, x₃)^T，γ=$\frac{1}{{\sqrt 3 }}$(1, 1, 1)^T，$A = \frac{{\sqrt 3 }}{2} - \frac{{1.645}}{{\sqrt {12} }}$，$C = \frac{{\sqrt 3 }}{2} + \frac{{1.645}}{{\sqrt {12} }}$；x_i i.i.d.~U(0, 1), i=1, 2, 3；Z~N(0, 1)，且x_i与Z相互独立。取样本容量n=200，在分位数τ分别取0.25、0.5、0.75的情况下估计γ，分别模拟100。选择高斯核函数$K(\nu ) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} }}{{\rm{e}}^{ - {\nu ^2}/2}}$，迭代精度δ取1×10^-6，带宽h_τ的选取如前所述。
类似参考文献[13, 19, 28]所用的评价标准，用Bias(样本偏差)、SE(样本标准差)、MSE(样本均方误差)、CI (95%的置信区间)、CP(真值落在置信区间的比例)评价单指标向量部分估计${\mathit{\boldsymbol{\hat \gamma }}}$的估计精度，用ASE和AAE评价联系函数部分?(·) 的估计精度, 并给出ASE与AAE的Mean(样本均值)与SE(样本均方误差)，其中ASE和AAE的定义如下：

$\mathrm{ASE}=n^{-1} \sum\limits_{i=1}^{n}\left[\hat{g}\left(u_{i}\right)-g\left(u_{i}\right)\right]^{2}, $

$\mathrm{AAE}=n^{-1} \sum\limits_{i=1}^{n}\left|\hat{g}\left(u_{i}\right)-g\left(u_{i}\right)\right| .$

在不同的分位数下，计算上述评价指标值，所得结果列于表 1和表 2，可以看出，用MM算法计算单指标分位回归模型，无论是单指标向量的估计还是联系函数的估计，都有比较好的结果，且在较小和较大的分位数下依然有良好的表现，这表明本文所提出的计算方法是有效的。将该估计结果与文献[13] 中关于该模型的模拟结果进行对比，可以看出，用MM算法计算单指标分位回归模型，与内点法相比，估计的精度是相似的，估计系数的偏差大小都在10^-3~10^-2，估计系数的样本标准差数量级均为10^-2，但在计算效率上，本文所提出的方法大大优于内点法，这将在模拟3中展示和说明。
Table 1

表 1 模拟1中不同的τ下，估计${\mathit{\boldsymbol{\hat \gamma }}}$的Bias、SE、MSE、CI、CPTable 1 The Bias、SE、MSE、CI、CP of ${\mathit{\boldsymbol{\hat \gamma }}}$ under different choices of τ in simulation 1 评价 标准 Bias SE MSE CI CP ${{\mathit{\boldsymbol{\hat \gamma }}}_1}$ 0.001 9 0.075 5 0.005 7 (0.431 1, 0.727 4) 0.94 τ=0.25 ${{\mathit{\boldsymbol{\hat \gamma }}}_2}$ -0.007 5 0.079 4 0.006 3 (0.414 1, 0.725 5) 0.93 ${{\mathit{\boldsymbol{\hat \gamma }}}_3}$ -0.003 6 0.095 2 0.009 1 (0.387 0, 0.760 3) 0.96 ${{\mathit{\boldsymbol{\hat \gamma }}}_1}$ 0.001 7 0.038 8 0.001 5 (0.503 0, 0.655 1) 0.94 τ=0.5 ${{\mathit{\boldsymbol{\hat \gamma }}}_2}$ 0.001 0 0.037 8 0.001 4 (0.502 0, 0.650 4) 0.94 ${{\mathit{\boldsymbol{\hat \gamma }}}_3}$ 0.005 9 0.034 9 0.001 2 (0.514 7, 0.651 8) 0.94 ${{\mathit{\boldsymbol{\hat \gamma }}}_1}$ -0.006 6 0.025 8 0.000 7 (0.520 1, 0.621 2) 0.95 τ=0.75 ${{\mathit{\boldsymbol{\hat \gamma }}}_2}$ 0.002 7 0.024 6 0.000 6 (0.531 8, 0.628 3) 0.93 ${{\mathit{\boldsymbol{\hat \gamma }}}_3}$ 0.006 1 0.006 0 0.000 6 (0.536 6, 0.630 2) 0.93

表 1 模拟1中不同的τ下，估计${\mathit{\boldsymbol{\hat \gamma }}}$的Bias、SE、MSE、CI、CPTable 1 The Bias、SE、MSE、CI、CP of ${\mathit{\boldsymbol{\hat \gamma }}}$ under different choices of τ in simulation 1

Table 2

表 2 模拟1中不同τ的选择下，?(u)的ASE、AAE的Mean、SETable 2 Outcomes of ?(u) for the models under different choices of τ in simulation 1 ?(u) τ=0.25 τ=0.5 τ=0.75 Mean SE Mean SE Mean SE ASE 0.142 1 0.020 1 0.043 6 0.004 8 0.060 7 0.006 1 AAE 0.340 9 0.021 1 0.175 1 0.009 1 0.168 6 0.009 3

表 2 模拟1中不同τ的选择下，?(u)的ASE、AAE的Mean、SETable 2 Outcomes of ?(u) for the models under different choices of τ in simulation 1

3.2 模拟2借鉴文献[13]模拟2的模型设置，考虑模型

$y=10 \sin (0.75 u)+\sqrt{1+\sin (u)} Z, $

其中u=x^Tγ，x=(x₁, x₂)^T，γ=(1, 2)^T/${\sqrt 5 }$，x_ii.i.d~N(0, 0.25²), i=1, 2；Z~N(0, 1)，且x_i与Z相互独立。取样本容量n=400，在分位数τ分别取0.25、0.5、0.75的情况下估计γ，分别模拟100次。其他设置同模拟1，所得结果列于表 3和表 4，可以看出新提出的算法在方差的分布不相同时，依然具有良好的估计结果，具有稳健性。将该估计结果与文献[13] 中对应模型的模拟结果对比，依然可以发现，用MM算法计算单指标分位回归模型，与内点法相比，估计的精度是相似的，再次表明我们所提出的方法是有效的。
Table 3

表 3 模拟2中不同的τ下，估计${\mathit{\boldsymbol{\hat \gamma }}}$的Bias、SE、MSE、CI、CPTable 3 The Bias、SE、MSE、CI、CP of ${\mathit{\boldsymbol{\hat \gamma }}}$ under different choices of τ in simulation 2 评价 标准 Bias SE MSE CI CP τ=0.25 ${{\mathit{\boldsymbol{\hat \gamma }}}_1}$ -0.014 1 0.118 1 0.014 1 (0.201 5, 0.664 5) 0.97 ${{\mathit{\boldsymbol{\hat \gamma }}}_2}$ -0.000 2 0.060 8 0.003 7 (0.775 3, 1.013 9) 0.93 τ=0.5 ${{\mathit{\boldsymbol{\hat \gamma }}}_1}$ -0.000 5 0.114 9 0.013 2 (0.221 3, 0.671 9) 0.93 ${{\mathit{\boldsymbol{\hat \gamma }}}_2}$ -0.006 2 0.057 3 0.003 3 (0.775 8, 1.000 3) 0.95 τ=0.75 ${{\mathit{\boldsymbol{\hat \gamma }}}_1}$ -0.004 4 0.108 4 0.011 7 (0.230 2, 0.655 3) 0.95 ${{\mathit{\boldsymbol{\hat \gamma }}}_2}$ -0.003 1 0.056 3 0.003 1 (0.780 9, 1.001 7) 0.92

表 3 模拟2中不同的τ下，估计${\mathit{\boldsymbol{\hat \gamma }}}$的Bias、SE、MSE、CI、CPTable 3 The Bias、SE、MSE、CI、CP of ${\mathit{\boldsymbol{\hat \gamma }}}$ under different choices of τ in simulation 2

Table 4

表 4 模拟2中不同τ的选择下，?(u)的ASE、AAE的Mean、SETable 4 Outcomes of ?(u) for the models under different choices of τ in simulation 2 ?(u) τ=0.25 τ=0.5 τ=0.75 Mean SE Mean SE Mean SE ASE 0.054 4 0.047 1 0.029 7 0.027 1 0.051 3 0.050 4 AAE 0.385 3 0.176 6 0.269 6 0.126 9 0.367 1 0.204 0

表 4 模拟2中不同τ的选择下，?(u)的ASE、AAE的Mean、SETable 4 Outcomes of ?(u) for the models under different choices of τ in simulation 2

3.3 与内点法比较由前两个模拟可以发现，MM算法在单指标分位回归估计的计算问题中表现良好，接下来比较MM算法与内点法的计算效率，考虑如下3个模型：
模型1: Y=exp(X′₁γ₁)+ε₁;
模型2: Y=5cos(X′₂γ₂)+ε₂;
模型3: Y=5sin(X′₃γ₃)+exp (-16(X′₃γ₃)²)+ε₃.
其中γ₁=(1, 1, 1)^T/${\sqrt 3 }$, X_{1, i}~N(0, 1), i=1, 2, 3, ε₁~N(0, 1);γ₂=(2, 1, 2)^T/3, X_{2, i}~U(0, 1), i=1, 2, 3, ε₂~N(0, 1);γ₃=(1, 2)^T/${\sqrt 5 }$, X_{3, i}~N(0, 1), i=1, 2, ε₃~Exp(1)。对每个模型，取样本量分别为100、200，模拟100次，用MM算法和内点法进行估计计算，并用100次模拟的总体运行时间与平均迭代次数比较算法的运行效率，所得结果列于表 5、表 6。
Table 5

表 5 n=100，模型1、2、3的估计结果比较Table 5 Estimation comparison among models 1, 2, and 3 with n=100 总体运行时间/h 平均迭代次数 model 1 MM 0.43 66.45 IP 2.20 62.26 model 2 MM 0.30 54.83 IP 0.65 17.27 model 3 MM 0.26 41.32 IP 0.54 17.36 注：IP指内点算法。

表 5 n=100，模型1、2、3的估计结果比较Table 5 Estimation comparison among models 1, 2, and 3 with n=100

Table 6

表 6 n=200，模型1、2、3的估计结果比较Table 6 Estimation comparison among models 1, 2, and 3 with n=200 总体运行时间/h 平均迭代次数 model 1 MM 1.51 56.97 IP 25.16 52.99 model 2 MM 1.09 50.73 IP 10.98 22.37 model 3 MM 0.95 37.98 IP 8.01 20.03

表 6 n=200，模型1、2、3的估计结果比较Table 6 Estimation comparison among models 1, 2, and 3 with n=200

可以发现MM算法所用的时间远远少于内点法，且随着样本量的增大，这种计算效率上的优势更加明显。这是由于用内点法解决分位回归问题，是将目标函数及约束条件转化为线性规划问题，再用内点法来求解该问题，但转化之后的线性规划问题，协变量维数与样本量的大小有着正相关的关系，这种方法极大地增加了算法的计算量与所用时间。本文第3.2节中的问题转化为线性规划后，协变量的维数为2n+p，第3.3节中的问题转化为线性规划后，维数为2n²+p，具体转化方法及维数的增加量可见文献[11]。而MM算法只需对p维矩阵做运算，故两种方法的计算效率随样本量的增加会产生越来越大的差距。
4 实例分析在这一节，应用本文所提的估计方法对Boston房屋数据集进行分析，Boston房屋数据集包含506个样本。该数据集及其详细介绍均可在R包mlbench(http://cran.r-project.org/) 找到。该数据集的响应变量MEDV(自住房屋房价的中位数)为右删失数据，当其大于50 000时，按50 000来记录，故比较适合用分位回归来分析。该数据集已被很多文献分析过，借鉴文献[13-14, 19] 的设置，选取MEDV(自住房屋房价的中位数)为响应变量，同时选取RM(每间住宅的房屋数目)，TAX(每一万美元的不动产税率)，PTRATIO(城镇中学生与教师的比例)，LSTAT(该地区房东收入较低的百分比)4个协变量，对TAX与LSTAT进行对数变换，并对响应变量做中心化处理，构造如下模型

$\begin{array}{c}\theta_{\tau}(\text { MEDV | RM, TAX, PTRATIO, LSTAT) = } \\g\left(\gamma_{1} \mathrm{RM}+\gamma_{2} \log (\mathrm{TAX})+\gamma_{3} \text { PTRATIO }\right)+ \\\gamma_{4} \log (\text { LSTAT }).\end{array}$

用本文提出的方法对该问题进行估计，计算在不同分位数下系数的估计值，并采用bootstrap方法估计标准差，方法如下，具体细节可参考文献[13, 29]。
1) 计算每个样本的误差, ${{\hat \in }_{\tau , i}} = {y_i} - {{\hat g}_\tau }\left( {\mathit{\boldsymbol{x}}_i^{\rm{T}}\mathit{\boldsymbol{\widehat \gamma }}} \right)$.再对误差进行中心化；
2) 对{${{\hat \in }_{\tau , i}}$}进行重抽样得到{$\hat \in _{\tau , i}^*$}；
3) 产生新的响应变量y_i^*=${{\hat g}_\tau }\left( {\mathit{\boldsymbol{x}}_i^{\rm{T}}\mathit{\boldsymbol{\widehat \gamma }}} \right)$+∈_{τ, i}^*；
4) 根据新的样本(x_i, y_i^*)，得到新的估计。
重复模拟100次计算标准差，所得结果列于表 7。从表 7可以发现，RM的系数在不同的分位数下皆为正，这表明每栋房屋的房间数量越多，房价就越高且收入越多的家庭更加在意每栋房屋的房间数量；log(TAX)的系数为负且随分位数逐渐变大，这表明不动产的税率越高，房价越低且收入较低的家庭更加在意不动产税率的大小；PTRATIO的系数为负且随分位数变化较小，这表明学生与教师的比例越大，即教师资源越匮乏，房价越低且低收入家庭与高收入家庭对教育的重视程度是同样大的；log(LSTAT) 系数为负且随分位数逐渐变小，这表明一个地区低收入人群所占的百分比越高，房价越低且收入较高的家庭更加在意一个地区的低收入人群比例。
Table 7

表 7 波士顿数据集在单指标分位回归模型下的系数估计及标准差估计Table 7 Coefficient estimation and standard deviation estimation of Boston data set under the single-index quantile regression mode τ RM log(TAX) PTRATIO log(LSTAT) 0.1 0.304 0 -0.610 7 -0.064 4 -0.567 2 (0.038 2) (0.051 2) (0.010 0) (0.052 3) 0.25 0.390 4 -0.508 3 -0.060 8 -0.533 3 (0.042 0) (0.045 1) (0.009 7) (0.048 7) 0.5 0.484 3 -0.446 7 -0.071 3 -0.614 6 (0.038 4) (0.044 6) (0.008 6) (0.043 8) 0.75 0.639 7 -0.338 8 -0.069 5 -0.592 0 (0.052 4) (0.067 1) (0.010 2) (0.061 9) 0.9 0.672 8 -0.078 0 -0.055 2 -0.810 6 (0.071 8) (0.089 6) (0.012 9) (0.076 8) 注：括号中数据为标准差估计。

表 7 波士顿数据集在单指标分位回归模型下的系数估计及标准差估计Table 7 Coefficient estimation and standard deviation estimation of Boston data set under the single-index quantile regression mode

图 1展示了不同分位数下g(·) 的估计，其中，横坐标为估计所得的${\mathit{\boldsymbol{\hat \gamma }}}$与协变量X的线性组合，纵坐标为对应的联系函数的估计值?(·)，散点为响应变量Y的实际观测值。观察图 1，可以发现散点(真实观测值)都紧紧围绕在估计函数?(·)的周围且估计曲线随分位数的增大而上移，这表明本文提出的方法应用在Boston数据集上是有效的。此外，还可以发现，在不同分位数下，所估计单指标x^Tγ_τ的取值范围不同，这说明采用分位回归分析该问题是非常有必要的。
Fig. 1

	Download: JPG larger image
图 1 联系函数g(u)及在不同的τ下, ?τ(u)的估计 Fig. 1 The link function g(u) and the estimation of ?τ(u) under different choices of τ 图 1 联系函数g(u)及在不同的τ下, ?τ(u)的估计 Fig. 1 The link function g(u) and the estimation of ?τ(u) under different choices of τ -->

5 结论本文研究单指标分位回归模型估计方法的MM算法。相比于内点法，MM算法极大地缩短了计算时间，提高了运算效率。此外，本文给出单指标分位回归模型在MM算法下的参数估计公式，在每次迭代过程中，将协变量与响应变量的观测值直接代入公式，即可得到参数的估计值，避免了每次迭代都要优化目标函数的麻烦。

参考文献

[1]	Zhao X, Wang W, Liu L, et al. A flexible quantile regression model for medical costs with application to Medical Expenditure Panel Survey Study[J]. Statistics in Medicine, 2018, 37(17): 2645-2666. Doi:10.1002/sim.7670
[2]	Petrella L, Raponi V. Joint estimation of conditional quantiles in multivariate linear regression models with an application to financial distress[J]. Journal of Multivariate Analysis, 2019, 173: 70-84. Doi:10.1016/j.jmva.2019.02.008
[3]	Tan X Z, Gan T Y, Chen S, et al. Modeling distributional changes in winter precipitation of Canada using Bayesian spatiotemporal quantile regression subjected to different teleconnections[J]. Climate Dynamics, 2019, 52: 2105-2124. Doi:10.1007/s00382-018-4241-0
[4]	Xu B, Lin B. Investigating the differences in CO2 emissions in the transport sector across Chinese provinces: evidence from a quantile regrssion model(Article)[J]. Journal of Cleaner Production, 2018, 27(175): 109-122.
[5]	姜成飞. 分位数回归方法综述[J]. 科技信息, 2013, 30(25): 185-240.
[6]	He X, Zhu L X. A lack-of-fit test for quantile regression[J]. Journal of the American Statistical Association, 2003, 98(464): 1013-1022. Doi:10.1198/016214503000000963
[7]	Mercedes C A, César S S, Wenceslao G M. A lack-of-fit test for quantile regression models with high-dimensional covariates[J]. Computational Statistics and Data Analysis, 2015, 33(88): 128-138.
[8]	关静. 分位数回归理论及其应用[D]. 天津: 天津大学, 2009.
[9]	朱平芳, 张征宇. 无条件分位数回归: 文献综述与应用实例[J]. 统计研究, 2012, 29(3): 88-96.
[10]	陈建宝, 丁军军. 分位数回归技术综述[J]. 统计与信息论坛, 2008, 23(3): 89-96.
[11]	Gao J Y. Computation in quantile and composite quantile regression models with or without regularization[D]. Edmonton: University of Alberta, 2015.
[12]	Hunter D R, Lange K. Quantile regression via an MM algorithm[J]. Journal of Computational and Graphical Statistics, 2000, 9(1): 60-77.
[13]	Wu T Z, Yu K M, Yu Y. Single-index quantitle regression[J]. Jornal of Multivariate Analysis, 2010, 40(101): 1607-1621.
[14]	Kong E F, Xia Y C. A single-index quantile regression model and its estimation[J]. Econometric Theory, 2012, 28(4): 730-768. Doi:10.1017/S0266466611000788
[15]	Kuruwita C N. Non-iterative estimation and variable selection in the single-index quantile regression model[J]. Communications in Statistics-Simulation and Computation, 2016, 45(10): 3615-3628. Doi:10.1080/03610918.2014.992542
[16]	Christou E. A Non-iterative method for fitting the single index quantile regression model with uncensored and censored data[D]. Pennsylvania: The Pennsylvania State University, 2016.
[17]	Hu Y A, Gramacy R B, Lian H. Bayesian quantile regression for single-index models[J]. Statistics and Computing, 2013, 23(4): 437-454. Doi:10.1007/s11222-012-9321-0
[18]	Alshaybawee T, Midi H, Alhamzawi R. Bayesian elastic net single index quantile regression[J]. Journal of Applied Statistics, 2017, 44(5): 853-871. Doi:10.1080/02664763.2016.1189515
[19]	Jiang R, Qian W M, Zhou Z G. Weighted composite quantile regression for single-index models[J]. Jornal of Multivariate Analysis, 2016, 148: 34-48. Doi:10.1016/j.jmva.2016.02.015
[20]	Liu H L, Yang H. Estimation and variable selection in single-index composite quantile regression[J]. Communications in Statistics-Simulation and Computation, 2017, 46(9): 7022-7039. Doi:10.1080/03610918.2016.1222424
[21]	Zhong W, Zhu L P, Li R Z, et al. Regularized quantile regression and robust feature screening for single index models[J]. Statistica Sinica, 2016, 26(1): 69-95.
[22]	Jiang Y L. Nonparametric quantile regression models via majorization minimization-algorithm[J]. Statistics and Its Interface, 2014, 7(2): 235-240. Doi:10.4310/SII.2014.v7.n2.a8
[23]	Xie S Y, Wan A T K, Zhou Y. Quantile regression methods with varying-coefficient models for censored data[J]. Computational Statistics and Data Analysis, 2015, 88: 154-172. Doi:10.1016/j.csda.2015.02.011
[24]	Chen Y L, Tang M L, Tian M Z. Semiparametric hierarchical composite quantile regression[J]. Communications in statistics: theory and methods, 2015, 50(4-6): 996-1012.
[25]	Lin W, Kulasekera K B. Indentifiability of single-index models and additive-index models[J]. Biometrica Jornal, 2007, 31(94): 496-501.
[26]	Chaudhuri P, Doksum K, Samarov A. On average derivative quantile regression[J]. Annals of Statistics, 1997, 25(2): 715-744.
[27]	Yu K M, Jones M C. Local linear quantile regression[J]. Journal of the American Statistic Association, 1998, 93(441): 228-237.
[28]	Lv Y Z, Zhang R Q, Zhao W H, et al. Quantile regression and variable selection of partial linear single-index model[J]. Ann Inst Stat Math, 2015, 67(2): 375-409. Doi:10.1007/s10463-014-0457-x
[29]	De Gooijer J G, Zerom D. On additive conditional quantiles with high-dimensional covariates[J]. Journal of the American Statistic Association, 2003, 98(116): 135-146.