基于γ散度的稳健小域估计方法

Robust Small Area Estimation Method Based on γ-Divergence

摘要: 被称为Fay-Herriot模型的两阶段正态分层模型被广泛应用于获得小范围内基于模型的均值估计。但是，当数据中出现异常值时，经验贝叶斯估计方法的性能可能很差，估计量会受到异常值的过度影响。本文提出了一种利用密度幂散度的修正方法，并提出了一种新的稳健经验贝叶斯小域估计方法，该方法对存在异常值的数据进行估计时有不错的效果。根据模型参数鲁棒估计量的渐近性质，推导了该估计量的均方误差并估计均方误差。通过数值模拟与实证分析，研究了该方法的数值性能，相比于其他估计量在处理异常值方面表现更为稳健。

Abstract: The two-stage normal hierarchical model known as the Fay-Herriot model is widely used to obtain model-based mean estimates for small areas. However, when outliers are present in the data, the performance of empirical Bayes estimation methods may be poor, and the estimators can be excessively influenced by outliers. This paper proposes a modification method utilizing density power divergence and introduces a new robust empirical Bayes small area estimation method, which performs well when estimating data with outliers. Based on the asymptotic properties of the robust estimator of model parameters, the mean squared error (MSE) of this estimator is derived and estimated. Through numerical simulations and empirical analyses, the numerical performance of this method is studied, demonstrating its robustness in handling outliers compared to other estimators.

1. 引言

在本文研究中，通常使用“小域”来表示无法产生足够精度的直接估计的任何领域。通常，领域样本大小倾向于随着领域的总体大小而增加，但情况并非总是如此。实践中，小域估计方法在政府统计、人口统计、医学统计、农业统计、贫困率估计等领域有广泛的应用[1]-[6]。在理论层面，小域估计也受到高度重视，从Fay提出的区域水平的Fay-Herriot模型[1]、Ghosh等提出在FH模型下的稳健贝叶斯预测方法[2]、Sinha和Rao等提出的稳健经验贝叶斯方法[3]、Jiang和Rao等关于稳健小域估计的概述[4]、Morales等提出的几类混合模型小区域估计理论[5]，小域估计已经形成了较为完整的理论体系。

为解决小域估计中FH模型中估计量受异常值影响的问题，Ghosh等提出了针对区域水平模型的稳健贝叶斯预测方法，并用于克服由异常值造成的过度收缩问题[2]，但并没有考虑异常值对模型系数的影响。Sinha和Rao利用Huber-Ψ函数，通过对由异常值引起的项减小权重的方式，建立了单元水平模型的小域稳健估计量[3]。目前，该方法普遍应用于小域稳健估计。在稳健估计的其他研究方面，Smith等用稳健投影和分位数的方法研究了商业调查中的小域稳健估计[7]，Bertarelli等通过偏差修正的方式研究了小域稳健估计问题[8]，Jiang和Rao对稳健小域估计的研究历史进行系统总结后指出其他现代统计方法也可以用于小域估计[4]。

密度幂散度的稳健估计是一种在统计学中用于处理异常值或重尾现象，从而得到更加稳定和可靠的参数估计的方法。对于将密度幂散度与小域估计相结合的相关理论在Sugasawa和Kubokawa的使用混合模型进行小域估计的综述文章中进行了系统的整理[9]，实现方法与环境构建也得到整理完善[10]。在相关理论研究中，Kim和Lee开发了一种稳健的估计方法，可以最小化散度族的指数多项式以达到简化计算的效果[11]，Basak和Basu通过最小化合适的统计距离对Bregaman散度进行拓展优化[12]，Sugasawa尝试了将散度应用到小域领域来改善区域均值的估计问题[13]，Nakagwa和Hashimoto提出了基于$\gamma$ 散度的稳健贝叶斯推断，在处理严重污染的数据时有不错的估计效果[14]。调优参数的选取也至关重要，调优参数的选择会影响目标函数对异常值的敏感度。较大的参数会使得目标函数对异常值更加敏感，而较小的参数则会降低这种敏感度。而对于调优参数的选择Sugasawa和Yonekura提出一种选择调优参数的准则对于后验分布的选择[15]，Basak等人以最小距离的选择方法处理调优参数的选择问题，但该方法对于污染严重的数据处理上没有很好的效果[16]，Yonekura和Sugasawa对于所提出的选择准则在针对一般贝叶斯的背景下进行了改善[17]。本文将基于$\gamma$ 散度与小域估计中的FH模型相结合来改善数据因为样本量不足且异常值存在导致出现对区域均值估计效果不佳的情况，在下文中的数值模拟与实证分析可以证实本文的方法对于改善该问题有不错的效果。

2. 模型及结合散度的贝叶斯预测

2.1. Fay-Herriot模型

Fay-Herriot模型是一种区域层次模型，已知，如果特定区域的样本量较小，仅基于特定区域的样本数据的直接调查估计量会产生过大的标准误差。经验贝叶斯方法被广泛用于改进直接估计，通过缩小一些综合估计和借用强度。关于小域估计与Fay-Herriot模型的全面概述，请参考Fay [1]和Ghosh [2]的文章。

一个基本的区域级模型的两阶段的层次模型，称为Fay-Herriot模型，描述为

$y_i|{\theta }_i~N\left({\theta }_i,D_i\right),{\theta }_i~N\left(x_i^{\text{T}}\beta ,A\right)\left(i=1,\cdots ,m\right)$ (1)

(1)式中，$y_i$ 为小域均值${\theta }_i$ 的直接估计量，$D_i$ 为抽样方差，假设抽样方差为已知，$x_i$ 和$\beta$ 分别为协变量和回归系数的向量，A为未知方差。设$f={\left(b^{\text{T}},A\right)}^{\text{T}}$ 为(1)式中的未知参数的向量。由于$y_i~N\left(x_i^{\text{T}}\beta ,A+D_i\right)$ 在(1)的模型下，$\varphi$ 可以通过最大化对数边际似然来估计。

$\log f\left(y;\varphi \right)=-\frac{m}{2}\log \left(2\pi \right)-\frac{1}{2}{\displaystyle \sum _{i=1}^m\log \left(A+D_i\right)}-\frac{1}{2}{\displaystyle \sum _{i=1}^m\frac{{\left(y_i-x_i^{\text{T}}\beta \right)}^2}{A+D_i}}$ (2)

其中$y={\left(y_i,\cdots ,y_m\right)}^{\text{T}}$ .在平方误差损失条件下，${\theta }_i$ 的贝叶斯预测量为：

${\tilde{\theta }}_i\left(y_i;\varphi \right)=y_i-\frac{D_i}{A+D_i}\left(y_i-x_i^{\text{T}}\beta \right)$ (3)

${\theta }_i$ 的经验贝叶斯估计量${\widehat{\theta }}_i\left(y_i;\varphi \right)={\tilde{\theta }}_i\left(y_i;{\widehat{\varphi }}_{\alpha }\right)$ 。

当$y_i$ 可以用辅助信息$x_i$ 很好地解释表达时，经验贝叶斯估计量可以估计出很好的效果。但是，$x_i$ 并不一定在所有区域都可以很好的解释表达直接估计量$y_i$ ；也就是说，在某些地区$y_i$ 可能离$x_i^{\text{T}}\beta$ 很远，本

文称之为异常观测。对于这类的观测，相应的${\theta }_i$ 的可以从与(1)式中与假设分布不同的分布中生成；也就是说，假设的${\theta }_i$ 的分布可能是错误的。在本文中考虑了一种异常观测存在的情况，并关注在这种情况下会出现以下两个问题。

问题1：贝叶斯预测量(3)可能将$y_i$ 向$x_i^{\text{T}}\beta$ 过度收缩。

问题2：基于(2)的估计量$\widehat{\varphi }$ 会受到异常观测值的过度影响。

这些问题已经在Fay [1]和Ghosh [2]等人的作品中得到了解决，但本文将通过使用密度幂散度扩展了这一主题的知识体系。

本文的中心思想是基于使用边际似然函数(2)的贝叶斯预测量(3)的替代表达式。贝叶斯预测量(3)可以写成：

${\tilde{\theta }}_i\left(y_i;\varphi \right)=y_i-D_i\frac{\partial }{\partial y_i}\log f\left(y_i;\varphi \right)$ (4)

如果$y_i|{\theta }_i~N\left({\theta }_i,D_i\right)$ 上述表达式成立；也就是说，只要边际似然函数$\log f\left(y_i;\varphi \right)$ 在${\theta }_i$ 分布不服从如上式(1)中的正态分布时，我们应该改变$\varphi$ 。从(4)式中可以看出，经典经验贝叶斯估计量${\theta }_i$ 可以通过最大化边际似然函数的方法然后求导数来确定。因此，本文将用密度幂散度代替对数边际似然的方法，其中包括Kullback-Leibler散度作为一种特殊情况。(1)式下的密度幂散度具有封闭形式，而本文的稳健贝叶斯预测量具有更简单的形式。本文还考虑了模型参数的稳健估计，并给出了它们的渐近性质。此外，本文构造了基于参数自举法的稳健经验贝叶斯估计法的均方误差估计量，并证明了其渐近有效性。

已有很多学者将广义似然用于贝叶斯推断，解决了观测值假设分布出现错误的情况，但本文处理的是(1)式中未直接观测到的区域均值${\theta }_i$ 的假设分布出现错误的情况。Ghosh等人[7]在Fay-Herriot模型(1)中提出了一种稳健的贝叶斯预测量，使用$\beta$ 的影响函数来处理属性1，但没有处理属性2。Sinha和Rao [3]提出利用Huber的函数推导出一般线性混合模型中的贝叶斯预测方法和参数估计的方法，从而解决了性质1和2；然而，所得的贝叶斯预测方法在处理属性1时具有局限性。

2.2. 与γ散度相结合的稳健贝叶斯估计量

对于$\gamma$ 散度，在统计学和数据分析中，$\gamma$ 散度通常用于构建稳健估计方法，以处理非正态性或存在异常值的数据。通过引入$\gamma$ 散度和$\gamma$ 似然函数，可以构建基于单元水平模型的小域稳健估计方法，从而得到更加稳健的估计结果。

设$X_1,\cdots ,X_n$ 是独立同分布的(i.i.d.)该随机变量是从概率密度函数$g\left(x\right)$ 中生成。设$f\left(x\right)$ 为基本的概率密度函数，$\delta \left(x\right)$ 为污染概率密度函数。假设$g\left(x\right)$ 是由$g\left(x\right)=\left(1-\epsilon \right)f\left(x\right)+\epsilon \delta \left(x\right)$ 给出的污染概率密度函数，其中$\epsilon$ 是污染的比值。在本文中，$\gamma$ 散度通常考虑污染密度函数$\delta \left(x\right)$ 主要位于密度函数$f\left(x\right)$ 尾部的情况。换句话说，对于一个异常值$x_0$ ，$\gamma$ 散度将认为$f\left(x_0\right)\approx 0$ 。本文考虑参数模型$f_{\theta }\left(x\right)=f\left(x;\theta \right)\left(\theta \in \Theta \right)$ 作为一个候选模型，其中$\Theta \subset {ℝ}^P$ 是$\theta ={\left({\theta }_1,\cdots ,{\theta }_n\right)}^{\text{T}}$ 的一个参数空间。

本文假设目标密度包含在候选模型$\left\{f_{\theta }|\theta \in \Theta \right\}$ 中，也就是说，$f\left(x\right)$ 可以用参数$f\left(x\right)=f_{{\theta }_0}\left(x\right)$ 来表示。在下文中，为了简单起见，本文将省略函数的参数。Nakagawa和Hashimoto所研究的关于概率密度$f_{\theta }$ 和g之间的散度如下：

$\begin{array}{c}{D}_{\gamma }\left(g,f_{\theta }\right)=\frac{1}{\gamma \left(\gamma +1\right)}\log {\displaystyle \int g^{1+\gamma }\text{d}x}-\frac{1}{\gamma }\log {\displaystyle \int g{f}_{\theta }^{\gamma }\text{d}x}+\frac{1}{\gamma +1}\log {\displaystyle \int f_{\theta }^{1+\lambda }\text{d}x}\\ =-d_{\gamma }\left(g,g\right)+d_{\gamma }\left(g,f_{\theta }\right)\end{array}$

其中$\gamma >0$ 是一个关于稳健性的调优参数，而$d_{\gamma }\left(g,f_{\theta }\right)$ 被称为$\gamma$ -交叉熵。这个散度也被称为“$\gamma$ -发散”。为了推导出$\theta$ 的最小$\gamma$ 散度估计量，本文可以考虑最小化问题：

$\underset{\theta \in \Theta }{\min }{d}_{\gamma }\left(g,f_{\theta }\right)=\underset{\theta \in \Theta }{\min }\left\{-\frac{1}{\gamma }\log {\displaystyle \int g{f}_{\theta }^{\gamma }\text{d}x}+\frac{1}{\gamma +1}\log {\displaystyle \int f_{\theta }^{1+\gamma }\text{d}x}\right\}$

对于$\theta$ 。虽然真实的密度函数$g\left(x\right)$ 是未知的，但本文注意到$\gamma$ -交叉熵可以通过对$d_{\gamma }\left(\overline{g},f_{\theta }\right)$进行经验估计，其中$\overline{g}$ 是$x_1,\cdots ,x_n$ 的经验概率密度。然后，通过$\arg {\min }_{\theta \in \Theta }{d}_{\gamma }\left(\overline{g},f_{\theta }\right)$ 得到$\theta$ 的稳健估计量。

现在，本文考虑如下的$\gamma$ -交叉熵的单调变换：

${\tilde{d}}_{\gamma }\left(g,f_{\theta }\right)=-\frac{1}{\gamma }\left\{\exp \left(-\gamma d_{\gamma }\left(g,f_{\theta }\right)\right)-1\right\}=-\frac{\frac{1}{\gamma }{\displaystyle \int g{f}_{\theta }^{\gamma }\text{d}x}}{{\left\{{\displaystyle \int f_{\theta }^{1+\gamma }\text{d}x}\right\}}^{\gamma /\left(1+\gamma \right)}}+\frac{1}{\gamma }$ (5)

注意到，(5)等式右侧的第二项，即“$+1/\gamma$ ”是下文证明命题1中的${\tilde{d}}_{\gamma }\left(g,f_{\theta }\right)$ 的必要条件。但在本文计算后验分布时，这一项在分母和分子中将被抵消不会对计算造成影响。这种转换在本文中是必要的。然后本文将给出了变换后的$\gamma$ -交叉熵的一些性质。

命题1. 设g和f是概率密度函数且设${\kappa }_1,{\kappa }_2$ 和$\kappa$ 是常数。那么${\tilde{d}}_{\gamma }\left(g,f\right)$ 具有以下性质：

(i) ${\tilde{d}}_{\gamma }\left({\kappa }_{1}g,{\kappa }_{2}f\right)={\kappa }_1{\tilde{d}}_{\gamma }\left(g,f\right)+\left(1-{\kappa }_1\right)/\gamma$。

(ii) ${\tilde{d}}_{\gamma }\left(g,f\right)={\tilde{d}}_{\gamma }\left(g,g\right)$ ，当且仅当$f=\kappa g$ 时成立。

(iii) ${\lim }_{\gamma \to 0}{\tilde{d}}_{\gamma }\left(g,f\right)=-{\displaystyle \int g\log f\text{d}x}=d_{KL}\left(g,f\right)$，其中$d_{KL}\left(g,f\right)$ 是g和f之间的Kullback-Leibler交叉熵。

证明(i)和(ii)的证明在这里省略，因为这两条性质是由$\gamma$ -交叉熵的定义直接证明。本文主要对性质(iii)进行证明，(iii)的证明可以通过泰勒展开式$f^{\gamma }=1+\gamma \log f+O\left({\gamma }^2\right)\left(\gamma \to 0\right)$ 来证明。

然后

$\begin{array}{l}\frac{1}{\gamma }{\displaystyle \int g{f}^{\gamma }\text{d}x}-1=\frac{1}{\gamma }{\displaystyle \int g\left(f^{\gamma }-1\right)\text{d}x}={\displaystyle \int g\log f\text{d}x}+O\left(\gamma \right)\left(\gamma \to 0\right),\\ {\left({\displaystyle \int f^{1+\gamma }\text{d}x}\right)}^{\gamma /\left(1+\gamma \right)}=1+O\left({\gamma }^2\right)\left(\gamma \to 0\right).\end{array}$

因此，有：

${\tilde{d}}_{\gamma }\left(g,f\right)=-\frac{1}{\gamma }\left[\frac{{\displaystyle \int g{f}^{\gamma }\text{d}x}}{{\left\{{\displaystyle \int f^{1+\gamma }\text{d}x}\right\}}^{\gamma /1+\gamma }}-1\right]\to -{\displaystyle \int g\log f\text{d}x}\left(\gamma \to 0\right)$

这样就完成了证明。

由经验密度函数$\overline{g}$ 代替真实密度g，于是可以得到：

$\begin{array}{c}-n{\tilde{d}}_{\gamma }\left(\overline{g},f_{\theta }\right)={\displaystyle \sum _{i-1}^n\frac{f_{\theta }{\left(X_i\right)}^{\gamma }}{\gamma {\left\{{\displaystyle \int f_{\theta }{\left(X_i\right)}^{1+\gamma }\text{d}x}\right\}}^{\gamma /\left(1+\gamma \right)}}}-\frac{n}{\gamma }\\ ={\displaystyle \sum _{i=1}^n{q}_{\theta }^{\left(\gamma \right)}\left(X_i\right)}-\frac{n}{\gamma }=Q_n^{\left(\gamma \right)}\left(\theta \right)\end{array}$ (6)

其中

$q_{\theta }^{\left(\gamma \right)}\left(x\right)=q^{\left(\gamma \right)}\left(\theta ;x\right)=\frac{1}{\gamma }{f}_{\theta }{\left(x\right)}^{\gamma }{\left\{{\displaystyle \int f_{\theta }{\left(x\right)}^{1+\gamma }\text{d}x}\right\}}^{-\gamma /\left(1+\gamma \right)}$

本文称成$Q_n^{\left(\gamma \right)}\left(\theta \right)$ 为$\gamma$ -似然，这是一种拟似然(或称加权似然)。根据命题1的性质(iii)可以得到：

$\underset{\gamma \to 0}{\lim }{Q}_n^{\left(\gamma \right)}\left(\theta \right)={\displaystyle \sum _{i=1}^n\log f_{\theta }\left(x_i\right)}$

所以，$\gamma$ -似然是对数似然的一种泛化。本文将用$\gamma$ -似然代替FH模型区域均值估计的对数边际似然。因为在确定模型的前提下使用$\gamma$ -似然所以可以得到：

$\underset{\gamma \to 0}{\lim }{Q}^{\gamma }\left(y;\varphi \right)=\log f\left(y;\varphi \right)$

在模型(1)下，$y_i$ 是独立的且服从$y_i~N\left(x_i^{\text{T}}\beta ,A+D_i\right)$ 的正态分布，因此(6)式可以用该形式表达

$Q^{\gamma }\left(y;\varphi \right)={\displaystyle \sum _{i=1}^m\left\{\frac{s_i\left(y_i;\varphi \right)}{\gamma }-\frac{i}{\gamma }\right\}}$ (7)

其中$V_i={\left\{2\pi \left(A+D_i\right)\right\}}^{-1/2}$ 且

$s_i\left(y_i;\varphi \right)={\left\{V_i^{-\gamma }{\left(1+\gamma \right)}^{1/2}\right\}}^{\gamma /\left(1+\gamma \right)}{V}_i^{\gamma }\exp \left(-\frac{\gamma {\left(y_i-x_i^{\text{T}}\beta \right)}^2}{2\left(A+D_i\right)}\right)$ 。

于是本文将$\gamma$ -似然函数(7)式中的$Q^{\gamma }\left(y;\varphi \right)$ 来代替(4)中对数边际似然$\log f\left(y;\varphi \right)$ 。来定义${\theta }_i$ 的鲁棒贝叶斯预测因子${\tilde{\theta }}_i^R$ 。可以得到：

$\frac{\partial }{\partial y_i}{Q}^{\gamma }\left(y;\varphi \right)=\frac{\partial }{\partial y_i}\frac{s_i\left(y_i;\varphi \right)}{\gamma }=\frac{1}{A+D_i}\left(y_i-x_i^{\text{T}}\beta \right)s_i\left(y_i;\varphi \right)$

则稳健贝叶斯预测为：

${\tilde{\theta }}_i^R=y_i-\frac{D_i}{A+D_i}\left(y_i-x_i^{\text{T}}\beta \right)s_i\left(y_i;\varphi \right)$ (8)

在(7)式中的收缩因子为$s_i\left(y_i;\varphi \right)/\left(A+D_i\right)$ ，它依赖于$y_i$ ，而经典的贝叶斯预测(3)中的收缩因子是$D_i/\left(A+D_i\right)$ ，它不依赖于$y_i$ 。而且，当$\gamma =0$ 时，$s_i\left(y_i;\varphi \right)=1$ 则${\tilde{\theta }}_i^R$ 将变为${\theta }_i$ 。此外，当$D_i\to 0$ ，稳健贝叶斯预测器${\tilde{\theta }}_i^R$ 将简化为直接估计量$y_i$ ，与经典的贝叶斯与测量${\theta }_i$ 一样。

3. 稳健经验贝叶斯估计量与均方误差

3.1. 稳健的参数估计

本文定义$\varphi$ 的稳健估计量${\widehat{\varphi }}_{\gamma }$ ，${\widehat{\varphi }}_{\gamma }=\arg \max Q^{\gamma }\left(y;\varphi \right)$，其中$Q^{\gamma }\left(y;\varphi \right)$ 在(7)中给出。则稳健估计量${\widehat{\varphi }}_{\gamma }$ 满足

$\begin{array}{l}\frac{\partial Q^{\gamma }}{\partial \beta }={\displaystyle \sum _{i=1}^m\frac{x_i{s}_i\left(y_i;\varphi \right)\left(y_i-x_i^{\text{T}}\beta \right)}{A+D_I}}=0,\\ \frac{\partial Q^{\gamma }}{\partial A}={\displaystyle \sum _{i=1}^m\left\{\frac{\gamma {\left(y_i-x_i^{\text{T}}\beta \right)}^2{s}_i\left(y_i;\varphi \right)-\gamma s_i\left(y_i;\varphi \right)/\sqrt{2\pi }}{2{\left(A+D_i\right)}^2}-\frac{{\gamma }^2{s}_i\left(y_i;\varphi \right)/\sqrt{2\pi }}{2\left(1+\gamma \right){\left(A+D_i\right)}^{3/2}}\right\}}.\end{array}$ (9)

本文使用Newton-Raphson来求解这些方程，求解过程在此省略详细求导过程请参考小域估计的文章[1]。本文将稳健估计量${\widehat{\varphi }}_{\gamma }$ 替换(7)中的贝叶斯估计量$\varphi$ ，最终得到稳健经验贝叶斯估计量${\widehat{\theta }}_i^R={\tilde{\theta }}_i^R\left(y_i;{\widehat{\varphi }}_{\gamma }\right)$。

3.2. 稳健的参数估计

调优参数$\gamma$ 与稳健性有关，但很难给出选择过程的解释。根据Ghosh [2]的文章，本文考虑基于稳健贝叶斯预测量的均方误差来选择(7)中的$\gamma$ ，这使下文能够以可解释的方式指定调优参数$\gamma$ 。均方误差的公式由以下定理给出。

定理1. 在模型(1)中，$E\left\{{\left({\tilde{\theta }}_i^R-{\theta }_i\right)}^2\right\}=g_{1i}\left(A\right)+g_{2i}\left(A\right)$ ，其中

$g_{1i}\left(A\right)=\frac{A{D}_i}{A+D_i},g_{2i}\left(A\right)=\frac{D_i^2}{A+D_i}\left\{\frac{V_i^{2\gamma }}{{\left(2\gamma +1\right)}^{3/2}}-\frac{2V_i^{\gamma }}{{\left(\gamma +1\right)}^{3/2}}+1\right\}$

其中$g_{2i}\left(A\right)$ 随$\gamma$ 的增大而增大。

经典贝叶斯预测量(3)的均方误差对应于$g_{1i}\left(A\right)$ ，因此${\tilde{\theta }}_i^R$ 相对于${\theta }_i$ 的超额均方误差为$g_{2i}\left(A\right)$ ，当$\gamma =0$ 时$g_{2i}\left(A\right)$ 接近于零。因此，${\tilde{\theta }}_i^R$ 的稳健性与均方误差之间模型(1)下存在着一个权衡。为了得到更好的调优参数$\gamma$ ，在此本文定义$Exc\left(\gamma \right)=100\%\times {\displaystyle {\sum }_{i=1}^m{g}_{2i}\left({\widehat{A}}_{\gamma }\right)}/{\displaystyle {\sum }_{i=1}^m{g}_{1i}\left({\widehat{A}}_{\gamma }\right)}$ 作为相对超额均方误差的一个百分比，其中${\widehat{A}}_{\gamma }$ 是A在式(11)中的稳健估计所得。本文对调优参数$\gamma$ 的选择遵循使$Exc\left(\gamma \right)$ 不超过用户指定的$c\%$ ；换句话说我们寻找$\gamma$ 满足$Exc\left({\gamma }^{*}\right)=c\%$ 。该等式存在唯一解，因为根据定理1可知$Exc\left(\gamma \right)$ 随$\gamma$ 的增大而增大。在实证分析中，我们首先计算指定$c\%$ 的调优参数$\gamma$ ，并将所有估计过程中的$\gamma$ 替换为$\gamma$ 进行的。而在下文中为了理论证明的简单性，我们假设$\gamma$ 是已知的，而在数值模拟与实证分析中我们将使用该调优参数的选择方法来研究它可能造成的影响。

3.3. 稳健估计量的渐进性质

我们考虑在模型(1)下稳健估计量的渐近性质。为此，下文需要假设以下正则性条件。

条件1：$0<D_{*}\le {\min }_{1\le i\le m}{D}_i\le {\max }_{1\le i\le m}{D}_i\le D^{*}<\infty$ ，其中$D_{*}$ 与$D_{*}$ 不依赖于m。

条件2：${\max }_{1\le i\le m}{x}_i^{\text{T}}{\left(X^{\text{T}}X\right)}^{-1}{x}_i=O\left(m^{-1}\right)$ ，其中$X={\left(x_1,\cdots ,x_m\right)}^{\text{T}}$ 。

条件3：当$m\to \infty$ 时，$X^{\text{T}}X/m$ 将收敛于一个正定矩阵。

Rao的文章中也假设了类似的条件。因为(9)式的导数在模型(1)下期望为0，我们将得到以下的结果。

定理2. 在条件1~3情况下，${\widehat{\beta }}_{\gamma }$ 与${\widehat{A}}_{\gamma }$ 是渐进独立的且分布分别为$N\left(\beta ,m^{-1}{J}_{\beta }^{-1}{K}_{\beta }{J}_{\beta }^{-1}\right)$ 与$N\left(A,K_A/m{J}_A^2\right)$ ，其中

$J_{\beta }=\frac{1}{m{\left(\gamma +1\right)}^{3/2}}{\displaystyle \sum _{i=1}^m\frac{V_i^{\beta }{x}_i{x}_i^{\text{T}}}{A+D_i}}$ ，$J_A=\frac{1}{2m}{\displaystyle \sum _{i=1}^m\frac{V_i^{\gamma }\left({\gamma }^2+2\right)}{{\left(A+D_i\right)}^2{\left(\gamma +1\right)}^{5/2}}}$ ，

$K_{\beta }=\frac{1}{m{\left(2\gamma +1\right)}^{3/2}}{\displaystyle \sum _{i=1}^m\frac{V_i^{2\beta }{x}_i{x}_i^{\text{T}}}{A+D_i}}$ ，${{\rm K}}_A=\frac{1}{m}{\displaystyle \sum _{i=1}^m\frac{V_i^{2\gamma }}{{\left(A+D_i\right)}^2}}\left\{\frac{2\left(2{\gamma }^2+1\right)}{{\left(2\gamma +1\right)}^{5/2}}-\frac{{\gamma }^2}{{\left(\gamma +1\right)}^3}\right\}$ 。

当$\gamma =0$ ，${\widehat{\beta }}_{\gamma }$ 的渐进协方差与${\widehat{A}}_{\gamma }$ 的

渐进方差将简化为Datta与Lahiri文章中所给出的$\beta$ 与A的极大似然估计量的渐进协方差和方差，因为(7)将简化为对数似然(2)。

3.4. 稳健经验贝叶斯估计量的均方误差

为了评估估计量${\widehat{\theta }}_i^R$ 的风险，下文考虑均方误差$M_i=E\left\{{\left({\widehat{\theta }}_i^R-{\theta }_i\right)}^2\right\}$ ，其中的期望是关于${\theta }_i$ 和

$y_i\left(i=1,\cdots ,m\right)$ 的联合分布的并遵循假设模型(1)。均方误差可视为综合贝叶斯风险，是小域估计中风险的标准度量。

由于${\widehat{\theta }}_i^R$ 依赖于估计量${\widehat{\varphi }}_{\gamma }$ ，均方误差$M_i$ 考虑了由于${\widehat{\varphi }}_{\gamma }$ 引起的额外波动。因此很难对$M_i$ 进行直接计算，所以使用了$M_i$ 的二阶近似。按照这个思路，下文用以下的定理给出了$M_i$ 的近似。

定理3. 在条件1~3下，

$M_i=g_{1i}\left(A\right)+g_{2i}\left(A\right)+\frac{g_{3i}\left(A\right)}{m}+\frac{g_{4i}\left(A\right)}{m}+\frac{2g_{5i}\left(A\right)}{m}+o\left(m^{-1}\right)$ (10)

其中$g_{1i}\left(A\right)$ 和$g_{2i}\left(A\right)$ 已经在定理1中给出且

$g_{3i}\left(A\right)=\frac{D_i^2{V}_i^{2\gamma }}{B_i^2{\left(2\gamma +1\right)}^{3/2}}{x}_i^{\text{T}}{J}_{\beta }^{-1}{K}_{\beta }{J}_{\beta }^{-1}{x}_i$ ，$g_{4i}\left(A\right)=\frac{D_i^2{V}_i^{2\gamma }{K}_A}{B_i^3{\left(2\gamma +1\right)}^{7/2}{J}_A^2}\left({\gamma }^4-\frac{1}{2}{\gamma }^2+1\right)$ ，

$\begin{array}{l}{g}_{5i}\left(A\right)=\frac{\gamma D_i^2{x}_i^{\text{T}}{J}_{\beta }^{-1}{K}_{\beta }{J}_{\beta }^{-1}{x}_i}{2B_i^4}\left(3B_i{C}_{11}-\gamma C_{21}\right)\\ +\frac{D_i^2{K}_A}{24B_i^6{J}_A^2}\left\{3\gamma B_i^2{C}_{21}+\left(\gamma -2\right)\left(3\gamma +8\right)C_{11}\right\}\\ +\frac{D_i^2{x}_i^T{J}_{\beta }^{-1}{x}_i}{B_i^4}\left(B_i{C}_{12}-\gamma C_{22}\right)+\frac{D_i^2}{2B_i^6{J}_A}\left\{\gamma C_{32}-2B_i{C}_{22}+\left(2-\gamma \right)B_i^2{C}_{12}\right\}\\ +\frac{D_i^2}{2B_i^4}\left\{b_A-\frac{\gamma V_i^{\gamma }}{{\left(\gamma +1\right)}^{3/2}{B}_i{J}_A}\right\}\left\{\left(2-\gamma \right)B_i{C}_{11}-\gamma C_{21}\right\}\end{array}$ 。

其中$B_i=A+D_i,b_A={\lim }_{m\to \infty }mE\left({\widehat{A}}_{\gamma }-A\right)$ 是${\widehat{A}}_{\gamma }$ 的一阶偏差，且

$C_{jk}=\left(2j-1\right)!!B_i^j\left\{V_i^{k\gamma }{\left(k\gamma +1\right)}^{-j-1/2}-V_i^{k\gamma +\gamma }{\left(k\gamma +\gamma +1\right)}^{-j-1/2}\right\}$ 。

其中$\left(2j-1\right)!!=\left(2j-1\right)\left(2j-3\right)\cdots \left(1\right)$ 。

这个定理的证明在Datta和Lahiri的补充材料中有类似的推导。近似公式(10)是基于已知的调优参数$\gamma$ ，因此在估计得到$\gamma$ 的情况下可能会有所不同。(10)中的表达式简化为Datta和Lahiri文章中给出的经

典经验贝叶斯估计量的均方误差，当$\gamma =0$ 时$C_{jk}=0$ ，$g_{5i}\left(A\right)=0$ 。

3.5. 均方误差的估计

由于定理3给出的均方误差的近似取决于未知参数A的选择，因此不能在实证分析中使用；所以使用均方误差的二阶无偏估计量。如果$E\left(\widehat{T}\right)=T+o\left(m^{-1}\right)$ ，则称估计量$\widehat{T}$ 是二阶无偏的。由定理3可知，$g_{3i}\left(A\right)$ ，$g_{4i}\left(A\right)$ 和$g_{5i}\left(A\right)$ 是光滑的，因此$g_{3i}\left(A\right)$ ，$g_{4i}\left(A\right)$ 和$g_{5i}\left(A\right)$ 是二阶无偏的。相比之下,$g_{1i}\left({\widehat{A}}_{\gamma }\right)$ 和$g_{2i}\left({\widehat{A}}_{\gamma }\right)$ 可能有相对大的偏差,因为$g_{1i}\left({\widehat{A}}_{\gamma }\right)$ 和$g_{2i}\left({\widehat{A}}_{\gamma }\right)$ 是$O\left(1\right)$ 。由于这些项的偏差和偏差校正估计的推导需要比较繁琐的代数计算，因此下文使用参数自举的方法，方法与Butar和Lahiri的文章中的类似。本文定义了自举估计量

${\widehat{M}}_i=2g_{12i}\left({\widehat{A}}_{\gamma }\right)-\frac{1}{B}{\displaystyle \sum _{b=1}^B{g}_{12i}\left({\widehat{A}}_{\gamma }^{\left(b\right)}\right)+\frac{g_{3i}\left({\widehat{A}}_{\gamma }\right)}{m}}+\frac{g_{4i}\left({\widehat{A}}_{\gamma }\right)}{m}+\frac{2g_{5i}\left({\widehat{A}}_{\gamma }\right)}{m}$ (11)

其中$g_{12i}\left({\widehat{A}}_{\gamma }\right)=g_{1i}\left({\widehat{A}}_{\gamma }\right)+g_{i2}\left({\widehat{A}}_{\gamma }\right)$ 且${\widehat{A}}_{\gamma }^{\left(b\right)}$ 是基于自举样本$y_i^{\left(b\right)},\cdots ,y_m^{\left(b\right)}$ 的自举估计量，自举样本是由下式产生$y_i^{\left(b\right)}=x_i^{\text{T}}{\widehat{\beta }}_{\gamma }+v_i^{\left(b\right)}+{\epsilon }_i^{\left(b\right)}$，$v_i^{\left(b\right)}~N\left(0,{\widehat{A}}_{\gamma }\right)$ ，${\epsilon }_i^{\left(b\right)}~N\left(0,D_i\right)$ 。根据Chang和Hall的文章，得到以下定理：

定理4. 在条件1~3的情况下设${\widehat{M}}_i^{†}$ 为${\widehat{M}}_i^{†}$ 在$B=0$ 时得到的理想版本。定义${\widehat{M}}_i={\widehat{M}}_i^{†}+U_i$ ，其中$U_i$ 表示仅执行有限数量的引导复制而产生的错误项。那么$E\left({\widehat{M}}_i^{†}\right)-M_i=o\left(m^{-1}\right)$ 且$U_i=O_p\left\{{\left(mB\right)}^{-1/2}\right\}$ 。

定理4表示，估计量(11)的理想版本${\widehat{M}}_i^{†}$ 是二阶无偏的。此外，该定理说明如果$B=O\left(m^{1+\delta }\right)$ 对于较小的$\delta >0$ ，则由有限次数的自举迭代引起的误差为$o_p\left(m^{-1}\right)$ ，因此估计量${\widehat{M}}_i^{†}$ 的表现应与理想估计量${\widehat{M}}_i^{†}$

相似。然而，理论情况下是基于已知的$\gamma$ ，因此自举估计量(11)不一定在估计所得的$\gamma$ 下被证明。

尽管在(10)中采用了一种附加形式的偏差校正，但在以Butar和Lahiri为基础也提出了其他形式的偏差校正。对于$g_{5i}\left(A\right)$ ，本文必须计算${\widehat{A}}_{\gamma }$ 的一阶偏差$b_A$ 的估计，它可以从参数自举样本中计算得到。换句

话说，本文可以使用参数自举法计算$g_{5i}\left(A\right)$ 。如定理3所示，$E\left\{\left({\widehat{\theta }}_i^R-{\tilde{\theta }}_i^R\right)\left({\tilde{\theta }}_i^R-{\tilde{\theta }}_i\right)\right\}=m^{-1}{g}_{5i}\left(A\right)+o\left(m^{-1}\right)$ ，所以下文可以用

$B^{-1}{\displaystyle \sum _{b=1}^B\left\{{\widehat{\theta }}_i^R\left(y_i^{\left(b\right)},{\widehat{\varphi }}_{\gamma }^{\left(b\right)}\right)-{\tilde{\theta }}_i^R\left(y_i^{\left(b\right)},{\widehat{\varphi }}_{\gamma }\right)\right\}}\left\{{\tilde{\theta }}_i^R\left(y_i^{\left(b\right)},{\widehat{\varphi }}_{\gamma }\right)-{\tilde{\theta }}_i\left(y_i^{\left(b\right)},{\widehat{\varphi }}_{\gamma }\right)\right\}$，

替代了$g_{5i}\left(A\right)$ ，其中${\widehat{\varphi }}_{\gamma }^{\left(b\right)}$ 是参数自举法估计量。

4. 数值模拟与实证分析

4.1. 数值模拟

首先本文评估了所提出的稳健估计方法的估计精度，并将其与现有估计方法进行了比较。本文考虑Fay-Herriot模型：

$y_i={\theta }_i+{\epsilon }_i$ ，$y_i={\theta }_i+{\epsilon }_i{\theta }_i={\beta }_0+{\beta }_1{x}_i+A^{1/2}{u}_i\left(i=1,\cdots ,m\right)$ 。

当$m=30$ ，${\beta }_0=0$ ，${\beta }_1=2$ ，$A=0.5$ ，且${\epsilon }_i~N\left(0,D_i\right)$ 。辅助变量$x_i$ 由$\left(0,1\right)$ 上的均匀分布生成。接下来，将m个区域分成5组，每组包含相同数量的区域，并在每组中设置相同的$D_i$ 值。各组的$D_i$ 值分别为$\left(0.2,0.4,0.6,0.8,1.0\right)$ 。对于$u_i$ 的分布，下文假设结构为$u_i~\left(1-\xi \right)N\left(0,\text{1}\right)+\xi N\left(0,{10}^2\right)$ ，其中$\xi$ 决定了假设分布受污染的程度。本文考虑以下三种情况：(I) $\xi =0$ ，(II) $\xi =0.15$ ，(III) $\xi =0.3$ 。在情景(II)和 (III)中，一些观测值有非常大的残差，辅助信息$x_i$ 对于这些异常观测值是有用的。

本文使用提出的具有密度幂散度的稳健经验贝叶斯估计方法估计${\theta }_i$ 。将使用两个不同的膨胀率分别为$c=1$ 和$c=5$ ，并根据上文中描述的选择准则来选择$\gamma$ 。比较了几种可供选择的方法：经典经验贝叶斯估计量、稳健贝叶斯估计量与使用Sinha和Rao提出的稳健估计方程使用极大似然方法估计的模型参数，以及Ghosh提出的稳健经验贝叶斯估计量使用模型参数的极大似然估计量。

表1报告了同一组内平均的均方误差值。由于经典经验贝叶斯方法中的正态性假设在情景(I)中是正确的，因此经验贝叶斯方法产生的均方误差比其他方法小这是自然的；然而，包括本文提出的方法在内的一些稳健方法的估计效果与经典方法相当。另一方面，在情景(II)和(III)中，当存在异常观测时，本文所提出的方法产生的均方误差往往小于其他方法，特别是对于抽样方差较大的群体。在这些场景中，当$c=5$ 时所提出的方法的性能要优于$c=1$ 时的方法，因为c越大，方法的稳健性越强。然而，本文所提出的方法的估计效果在情景(I)中，$c=5$ 比$c=1$ 更差，这可以被视为为情景(II)和(III)中估计效果有更强的稳健性所付出的合理代价。

4.2. 实证分析

我们考虑应用美国劳工统计局的鲜奶消费数据，该数据在Arora和Lahiri的文章与You和Chapman的文章中都进行了研究。在该数据集中，有43个地区图1显示了$y_i$ 的散点图以及${\beta }_1,\cdots ,{\beta }_4$ 的最大似然估计量。表明在$R_1$ 和$R_2$ 区域存在一些异常值。为了证实这一点，计算本文了标准化残差

$r_i={\left(\widehat{A}+D_i\right)}^{-1/2}\left\{y_i-{\displaystyle \sum _{j=1}^4{\widehat{\beta }}_{j}I\left(i\in M_j\right)}\right\}\left(i=1,\cdots ,m\right)$。

当正确指定Fay-Herriot模型(1)时，$r_i$ 的分布接近标准正态分布。然而，如表2所示，在某些地区，$r_i$ 的绝对值很高。本文使用Sinha和Rao的稳健估计方法和提出的密度幂散度方法作为(9)在1%和5%膨胀率下的解来估计参数。

从表2可以看出，在$r_i$ 绝对值较大的区域，${\widehat{\theta }}_i^{EB}$ 与${\widehat{\theta }}_i^R$ 的差异较大。${\widehat{M}}_i^{EB}$ 和${\widehat{M}}_i$ 之间的关系也出现了类似的情况。相反，${\widehat{\theta }}_i^{SR}$ 和$M_i^{SR}$ 的值与其他值不同，这可能是由于表2中A的估计值较低。

Table 1. The average mean square error of different estimation methods within the same group; All values are multiplied by 1000

表1. 不同估计方法的均方误差在同一组内的平均值；所有的值都乘以1000

情景	组别	DEB1	DEB2	EB	REB1	REB2	GEB
(I)	1	141	142	156	174	158	158
	2	218	225	251	281	259	306
	3	269	282	316	338	325	368
	4	307	326	350	376	359	393
	5	337	364	387	404	388	417
(II)	1	190	186	193	329	189	189
	2	360	345	372	1021	362	364
	3	513	482	604	1836	528	533
	4	662	620	879	2710	668	684
	5	798	745	860	3598	813	828
(III)	1	197	195	210	227	195	195
	2	389	383	391	515	385	385
	3	575	561	579	953	565	569
	4	758	737	773	1522	770	767
	5	937	905	947	2249	911	913

DEB1为膨胀率为1%时密度幂散度，即$c=1$ ；DEB2为膨胀率为5%时密度幂散度，即$c=5$ ；EB为经验贝叶斯方法；REB1为Sinha和Rao的稳健经验贝叶斯方法；REB2，为Sinha和Rao的最大似然稳健贝叶斯预测方法；GEB，为Ghosh等人的稳健经验贝叶斯方法。

Figure 1. Scatter plot of average expenditure on fresh milk in 43 regions, along with the maximum likelihood estimation value

图1. 43个地区鲜奶平均支出的散点图，以及最大似然估计值

Table 2. Empirical Bayesian estimators, the robust empirical Bayesian estimators proposed in this paper, the robust empirical Bayesian estimators of Sinha and Rao, and the mean squared error estimates of these three methods; Multiply the values of${\widehat{M}}_i^{EB}$ ，${\widehat{M}}_i$ and ${\widehat{M}}_i^{SR}$ by 100

表2. 经验贝叶斯估计量，本文提出的稳健经验贝叶斯估计量与Sinha和Rao的稳健经验贝叶斯估计量，以及该三种方法对均方误差的估计；将${\widehat{M}}_i^{EB}$ ，${\widehat{M}}_i$ 和${\widehat{M}}_i^{SR}$ 的值乘以100

地区	区域	$y_i$	$r_i$	${\widehat{\theta }}_i^{EB}$	${\widehat{\theta }}_i^R$	${\widehat{\theta }}_i^{SR}$	${\widehat{M}}_i^{EB}$	${\widehat{M}}_i$	${\widehat{M}}_i^{SR}$
3	1	1.10	0.63	1.02	1.02	1.02	1.35	1.35	0.81
4	1	0.63	−2.05	0.78	0.76	0.91	0.85	0.85	4.70
5	1	0.75	−1.25	0.86	0.87	0.92	0.96	0.96	4.44
10	2	1.36	1.38	1.19	1.24	1.23	1.39	1.37	3.79
11	2	0.62	−3.01	0.80	0.73	1.07	0.77	0.78	5.06
12	2	1.46	1.54	1.20	1.24	1.23	1.63	1.62	5.66
19	3	1.28	0.47	1.22	1.22	1.22	1.30	1.32	0.77
22	3	1.18	−0.01	1.19	1.19	1.19	0.80	0.84	0.56
31	4	0.89	0.63	0.76	0.76	0.75	1.54	1.63	0.78
37	4	0.44	−1.84	0.54	0.54	0.61	0.64	0.65	3.59

1989年鲜奶平均支出的估计值$y_i$ ，以及抽样方差$D_i$ 。继Arora和Lahiri的研究之后，我们考虑如果第i个区域属于第j个区域的话，我们有$x_i^{\text{T}}\beta ={\beta }_j\left(j=1,\cdot \cdot \cdot ,4\right)$ 的Fay-Herriot模型(1)。地区分成四个区域分别是$R_1=\left\{1,\cdot \cdot \cdot ,7\right\}$ ，$R_2=\left\{8,\cdot \cdot \cdot ,14\right\}$ ，$R_3=\left\{15,\cdot \cdot \cdot ,25\right\}$ 和$R_4=\left\{26,\cdot \cdot \cdot ,43\right\}$ 。

5. 结论与展望

通过数值模拟与实证分析表明，与现有方法相比，本文提出的方法在数据中存在异常观测的情况下特别适用。本研究聚焦于稳健小域估计问题，主要的包含稳健小域估计方法的比较、密度幂散度在小域估计模型中的应用的研究。提出了一种针对存在数据存在异常观测值的Fay-Herriot模型的稳健小域估计方法。本文通过引入$\gamma$ 散度，给出了解决具有异常观测的稳健估计方法。对主要结论进行归纳，包括如下两个方面：

第一，对基本的小域估计方法、小域稳健估计方法进行了比较研究。通过对现有小域估计方法的分析发现，小域估计量对区域随机效应的方差较为敏感。当观测数据存在异常值或者随机效应是有偏分布时，现有估计量具有较大的偏差，且容易受异常值影响。而现有的稳健估计方法中，一部分方法是针对于特殊情形下的小域模型进行估计的方法，不具有普适性。另一部分稳健估计方法容易出现过度收缩的情形。因此本文中引入了密度幂散度族的估计方法，给出了$\gamma$ 散度的定义、相关的性质以及如何进行估计等。

第二，本文研究了区域水平模型的稳健小域估计。首先将DPD散度引入了区域水平模型，通过推广似然函数的表达，给出了稳健估计量的具体公式。结合估计量的表达式，给出了估计量的MSE及其估计。其次，将$\gamma$ 散度应用于区域水平模型，给出了估计量的表达式和MSE。最后，给出了模型调整参数的选择算法和估计量MSE估计的程序。通过模拟数据和实际数据，验证了本研究中提出方法的优越性。结果一致表明在适当的调整参数中，本文提出的方法对于模型的系数的估计和小域目标的估计均具有较小的MSE，表现出了稳健性。