基于Fiducial预测密度的模型选择

Model Selection Based on Fiducial Predictive Distribution

摘要: 模型选择作为统计分析的一个重要工具，它能选出候选模型集中拟合数据生成过程最好的那个模型。为了得到更好的预测标准，本文提出了一种基于Fiducial预测密度的模型选择方法，同时加入容许集来进一步压缩候选模型集，给出理论的同时也给出了MH算法去应用它。最后，对本文提出的模型选择进行了模拟研究与实例分析，结果均表明我们的方法优于其他方法。

Abstract: Model selection is an essential tool in statistical analysis, enabling the identification of the best-fitting model for the data-generating process from a set of candidate models. To achieve better predictive performance, this paper proposes a model selection method based on Fiducial predictive density, incorporating an admissible set to further reduce the candidate model space. Theoretical foundations are provided, along with a Metropolis-Hastings (MH) algorithm for its implementation. Finally, simulation studies and empirical analyses are conducted to evaluate the proposed method. The results demonstrate that our method outperforms existing approaches in terms of both predictive accuracy and efficiency.

1. 引言

随着大数据时代的到来，数据分析已经成为各行各业不可或缺的工作之一，如何有效地利用统计建模从数据中挖掘出有用的信息也越来越受到人们的关注。模型选择可以在建模过程中找到对响应变量最具有解释性的自变量集，提高模型的解释性和预测精度，因而也随之备受关注。模型选择在统计学与计量经济学中有着悠久的历史，研究者们提出了不同的模型选择方法与准则。目前可将其分为两大类：(1) 以优化理论为支撑的正则化模型选择；(2) 通过构建候选模型的概率描述来进行模型选择。

正则化变量选择因其计算速度快和对超高维度数据具有良好的适应性而得到迅速发展。常见的方法有Lasso [1]、SCAD [2]、EN [3]、MCP [4]等方法，它们具有良好的稀疏性和收敛性。但以优化理论为支撑的正则化方法并不能给出所选模型为最优模型的概率保证，而且往往会导致过拟合现象。

量化的概率保证可以给从业者直观指导，这在实际应用中是十分重要的。其中一种概率保证是贝叶斯模型选择，有基于边际似然的最高模型后验准则BIC [5]和基于偏差信息量的最高模型后验准则DIC [6]，但前者计算复杂度往往很大且难以处理超高维问题，后者不能就近区分数据生成模型和过拟合模型[7]。另一种是频率学变量选择方法，常见的方法有AIC方法[8]，TIC方法[9]、Mallows Cp方法[10]、交叉验证[11]等，但它们都过于依赖参数个数。

近年来，基于Fiducial推断理论的模型选择开始得到关注。Fiducial推断的思想是一种介于频率学派与贝叶斯学派的新思想，可以追溯到由Fisher [12]提出了单参数的Fiducial分布的概念开始，并建议用Fiducial分布代替Bayes后验分布来获取参数的区间估计。对于单参数分布族，Fisher的Fiducial置信区间与经典置信区间重合；对于多参数分布族，Fiducial置信集的覆盖率接近名义水平，但在重复抽样频率意义上不精确。随后Fiducial推断得到进一步的发展，这其中包括Dempster-Shafer理论[13]，推断模型[14]，置信分布[15]等。这些现代Fiducial理论极大地推广了Fisher的Fiducial思想。Majumder和Hannig [16]建立了广义Fiducial推断的高阶渐近定义。Long和Xu [17]在研究模型分类问题中提出了Fiducial预测密度(FPD)的概念，并且证明了FPD拥有很好的依概率收敛的性质。

Fiducial推断理论的快速发展，使得基于Fiducial推断构建模型选择准则成为可能。该准则不仅避免了贝叶斯模型选择方法依赖于先验信息的问题，还可以给出所选模型为最优模型的概率保证，并且为研究模型选择准则的大样本性质提供可能。Hannig和Lee [18]在2009年首次将Fiducial推断理论与模型选择相结合，导出了Fiducial模型选择范式，其主要思想是将模型M作为未知参数，引入结构方程中。但直接通过模型的Fiducial概率来选择最优模型往往会出现较为严重的过拟合现象，因此，Hannig在此基础上施加最小描述长度惩罚(Minimum description length, MDL)使得结果有一定的改观。李涵等[19]提出基于Fiducial推断的Kriging模型变量选择方法，给出了如何从一般Kriging和正交Kriging模型选出最优模型。张淑芹等[20]给出了拟Fiducial推断的Kriging模型选择方法，发现相比于将Lasso和EN应用进Kriging模型进行模型选择，拟Fiducial推断的模型选择方法具有更高的拟合准确性和预测精度。

为了减少候选模型的个数，从而在增强真实模型所反映的信号的同时，进一步减轻计算负担，Williams等[21]引入了ε-容许集的概念。该方法极大地压缩了候选模型集合的数量。这项工作带来了一种新的Fiducial模型选择方法，被称为EAS (ε-admissible set，ε-容许集)。Williams等[22]进一步将EAS方法推广到了向量自回归模型。赵勇超等[23]提出了一种基于Fiducial边际似然函数的模型不确定性度量新方法FML，并将所提出的模型不确定性度量整合到EAS中，提出了一种Fiducial变量选择准则FMC，该准则可以很好地处理高维甚至超高维的情形。

本文将基于FPD函数提出一个新的模型不确定性度量方法，该方法在进行模型选择时同时考虑了参数的Fiducial分布和未来观测值的Fiducial预测密度，同时考虑进ε-容许集来加强对候选模型集的压缩，并记该模型选择方法为FPMC，同时给出了MH算法去实现它。本文结构如下：第二节介绍了Fiducial预测模型选择准则FPMC，第三节通过数值模拟及实例分析来说明本文方法的优越性。第四节是对本文的总结与展望。

2. Fiducial预测密度模型选择

2.1. 高维线性回归模型

在本节中给出本文所考虑的模型及问题，并对部分符号作出说明。考虑$n$ 个独立的观测$\left(x_i,y_i\right)$ ，$i=1,\cdots ,n$ 其中$x_i\in {ℝ}^p$ ，$y_i\in ℝ$ 。记$y={\left(y_1,\cdots ,y_n\right)}^{\prime }$ 。

令$S=\left\{M|M\subset \left\{1,2,\cdots ,p\right\}\right\}$ 是所有候选模型所构成的集合，则集合$S$ 中总共包括了$2^p$ 个候选模型。$\left|M\right|$ 是模型$M$ 的模型长度。对于候选模型$M$ ，本文考虑回归系数为$\beta$ ，误差项为${\epsilon }_M\sim N_n\left(0,{\sigma }_M^2{I}_n\right)$ 的多元正态线性回归模型，即

$y=X_M{\beta }_M+{\epsilon }_M$ (1)

其中设计阵$X_M$ 定义为仅由索引集$M$ 对应的$X$ 的那些列组成的矩阵，${\beta }_M$ 是由$M$ 索引的$\beta$ 对应的$\left|M\right|$ 维向量。令${\theta }_M=\left({\beta }_M,{\sigma }_M^2\right)\subset {\Theta }_M$ ，其中${\Theta }_M$ 是第$M$ 个模型的参数空间。

在本文的后续研究中，允许协变量个数大于观测值个数，即$p>n$ ，但此时回归模型的回归系数并不存在唯一的最小二乘估计。事实上，现有文献的各种高维变量选择问题的研究主要都是施加一定的稀疏性假设，使得高维问题在数学上可处理。尽管稀疏性假设有时并不一定正确，但就目前而言这是必要的。为了方便在数学框架下讨论高维变量选择问题，本文研究也是在一定的稀疏性假设下进行的。

令$M_o$ 为真实模型，假设$M_o\in S$ ，即候选模型集合中包含真实模型，则模型选择的目标分为两种，第一种是识别，即在$2^p$ 个候选模型中识别出真实模型$M_o$ 。第二种是预测，即并不需要识别出真实模型，只需要找到可以达到最优预测效果的模型即可，该模型可能只是真实模型的一个有效替代。相比较后者，前者的要求更高，因为完全识别出真实模型是很难的，事实上在实际案例中无法明确哪一个模型是真实模型，即$M_o\notin S$ 。但后者所选的替代模型可能是个非常复杂的模型，即存在较多的冗余信息。为此，本文综合的来考虑两个目的，即在识别效果可比的前提下，是否可以达到较优的预测效果。

2.2. 参数的Fiducial密度

Fiducial推断的基本思想源于对Fisher提出的Fiducial思想的理解。对于模型$M$ ，首先将数据$y$ 和参数${\theta }_M$ 的关系表示为

$y=G\left(U,{\theta }_M\right)$

其中$G$ 是一个确定的函数，称为数据生成方程，$U$ 是分布已知的随机部分，并且与参数${\theta }_M$ 独立。给定$y$ 的情况下，如果对于任意的${\theta }_M$ ，$G^{-1}$ 均存在，则

${\theta }_M=G^{-1}\left(y,u\right)$

为了简单起见，本文只考虑简单情况，即数据生成函数$G$ 可逆，对于不可逆的情况，Hannig [24]也给出了相应的解决方法，在这里不多赘述。在已知$U$ 的分布下，多次重复抽样$U$ ，通过结构方程的逆可以得到一组关于${\theta }_M$ 的随机样本，我们称之为${\theta }_M$ 的Fiducial样本，与之对应的密度称为${\theta }_M$ 的Fiducial密度，记为$r\left({\theta }_M|y\right)$ 。

对于线性回归模型$M$ ，参数${\theta }_M=\left({\beta }_M,{\sigma }_M^2\right)$ 的充分统计量为

$\left({\widehat{\beta }}_M,{\widehat{\sigma }}_M^2\right)=\left({\left({X^{\prime }}_M{X}_M\right)}^{-1}{X^{\prime }}_{M}Y,\frac{RS{S}_M}{n-\left|M\right|}\right)$

由抽样分布定理和基本的数理运算可得

${\widehat{\beta }}_M={\beta }_M+{\sigma }_M{\left({X^{\prime }}_M{X}_M\right)}^{-\frac{1}{2}}Z$

$\left(n-\left|M\right|\right){\widehat{\sigma }}_M^2={\sigma }_M^{2}U$

其中$Z\sim N_{\left|M\right|}\left(0,I_{\left|M\right|}\right)$ ，$U\sim {\chi }_{n-\left|M\right|}^2$ ，并且两者相互独立。当给定观测值$y$ 和$\left(z,u\right)$ 时，该结构方程的逆存在唯一解：

${\beta }_M={\widehat{\beta }}_M+{\sigma }_M{\left({X^{\prime }}_M{X}_M\right)}^{-\frac{1}{2}}z$ ，${\sigma }_M^2=\frac{\left(n-\left|M\right|\right){\widehat{\sigma }}_M^2}{u}$

因此${\theta }_M=\left({\beta }_M,{\sigma }_M^2\right)$ 的联合Fiducial密度为

$r\left({\theta }_M|y\right)=\frac{{\pi }^{-\frac{\left|M\right|}{2}}}{2^{\frac{n}{2}}\Gamma \left(\frac{n-\left|M\right|}{2}\right)}RS{S}_M{\left|{X^{\prime }}_M{X}_M\right|}^{\frac{1}{2}}{\left({\sigma }_M^2\right)}^{-\frac{n}{2}-1}\exp \left(-\frac{\Vert y-X_M{\beta }_M\Vert }{2{\sigma }_M^2}\right)$

而且容易求得${\sigma }_M^2$ 的边际密度函数是$IGa\left(\frac{n-\left|M\right|}{2},\frac{RS{S}_M}{2}\right)$ 的概率密度函数，为

$r\left({\sigma }_M^2|y\right)=\frac{{\left(\frac{RS{S}_M}{2}\right)}^{\frac{n-\left|M\right|}{2}}}{\Gamma \left(\frac{n-\left|M\right|}{2}\right)}{\left({\sigma }_M^2\right)}^{-\frac{n-\left|M\right|}{2}-1}\exp \left(\frac{RS{S}_M}{2{\sigma }_M^2}\right)$

${\beta }_M$ 的边际密度函数为

$\begin{array}{l}r\left({\beta }_M|y\right)={\displaystyle {\int }_0^{+\infty }r}\left({\theta }_M|y\right)\text{d}{\sigma }_M^2\\ \text{ =}\frac{\Gamma \left(\frac{n}{2}\right)}{\Gamma \left(\frac{n-\left|M\right|}{2}\right)}\frac{RS{S}_M^{\frac{n-\left|M\right|}{2}}}{{\pi }^{\frac{\left|M\right|}{2}}}{\left|{X^{\prime }}_M{X}_M\right|}^{\frac{1}{2}}{\left({\Vert y-X_M{\beta }_M\Vert }_2^2\right)}^{-\frac{n}{2}}\\ \text{ =}\frac{\Gamma \left(\frac{n}{2}\right)}{\Gamma \left(\frac{n-\left|M\right|}{2}\right)}\frac{RS{S}_M^{\frac{-\left|M\right|}{2}}}{{\pi }^{\frac{\left|M\right|}{2}}}{\left|{X^{\prime }}_M{X}_M\right|}^{\frac{1}{2}}{\left[1+\frac{RS{S}_M{\Vert X_M\left({\beta }_M-{\widehat{\beta }}_M\right)\Vert }_2^2}{{\left(n-\left|M\right|\right)}^2}\right]}^{-\frac{n}{2}}\end{array}$ (2)

即${\beta }_M$ 服从多元$t$ 分布$t_{n-\left|M\right|}\left({\widehat{\beta }}_M,\frac{RS{S}_M}{n-\left|M\right|}{\left({X^{\prime }}_M{X}_M\right)}^{-1}\right)$ 。

2.3. 候选模型集压缩方法

在高维模型下经常会出现$p\gg n$ 的情况，包括Fiducial模型选择在内的概率描述类模型择方法无法实现。其原因一方面是遍历整个候选模型集合来分别计算各自的模型概率度量具有极高的计算复杂度；另一方面，模型概率度量在整个候选模型集合上取值为1，当候选模型集中包含了较多模型时，包括真实模型在内的每个模型所分摊的概率将会很小，这导致真实模型所反映出的信号变的很弱，从而无法达到模型选择的目的。为此合理且有效的压缩候选模型集是必要的。

Williams和Hannig [21]提出如下的ε-容许集的概念。

给定$\epsilon >0$ ，一个带有指标集$M\subset \left\{1,2,\cdots ,p\right\}$ 的系数向量${\beta }_M$ 被称为是ε-容许的，当且仅当$h\left({\beta }_M\right)=1$ ，其中

$h\left({\beta }_M\right):=I\left(\frac{1}{2}{\Vert X^{\prime }\left(X_M{\beta }_M-X{b}_{\min }\right)\Vert }_2^2\ge \epsilon \right)$

其中$b_{\min }$ 是在${\Vert b\Vert }_0\le \left|M\right|-1$ 的条件下，下述优化问题的解

$\underset{b\in {ℝ}^p}{\min }\frac{1}{2}{\Vert X^{\prime }\left(X_M{\beta }_M-Xb\right)\Vert }_2^2$

候选模型集的压缩效果依赖于$\epsilon$ ，即一个较小的$\epsilon$ 会产生一个相对稀疏的容许集合，但过大的$\epsilon$ 可能会导致真实模型都变为不被容许。$\epsilon$ 偏大或偏小都会降低模型选择的精确度。在本文中，建议使用Williams和Hannig [21]相同的默认规则来确定$\epsilon$ ：

${\epsilon }_M={\Lambda }_M{\widehat{\sigma }}_M^2{\left(\frac{n^{0.51}}{9}+\left|M\right|\frac{\log {\left(p\pi \right)}^{1.1}}{9}+1-p_o\right)}_{+}$

其中${\Lambda }_M:={\Vert H_{M}X\Vert }_{Frobenius}^2$ ，$H_M:=X_M{\left({X^{\prime }}_M{X}_M\right)}^{-1}{X^{\prime }}_M$ 以及${\widehat{\sigma }}_M^2:={\text{RSS}}_M/\left(n-\left|M\right|\right)$。$p_o$ 反映的是真实模型$M_o$

的稀疏程度，我们用Zhu等人[25]提出的自适应最佳子集选择(Adaptive Best Subset Selection, ABESS)算法来获得真实模型$M_o$ 的稀疏度估计。不同数据集下的线性回归都可以适用该$\epsilon$ 规则。

2.4. Fiducial预测密度模型选择准则

类似于Bayes预测密度的定义，Long和Xu [17]将参数的Fiducial密度替代Bayes后验密度，给出了Fiducial预测密度定义。

对于线性回归模型(1)，如果我们已知$\left|M\right|$ 维向量$x$ 的数值，且$x$ 独立于设计矩阵$X_M$ ，则预测变量$z$ 与$y$ 都相互独立，并且$z$ 的分布为$N\left(x^{\prime }{\beta }_M,{\sigma }_M^2\right)$ ，即$z$ 的概率密度函数为

$f\left(z|{\theta }_M,M\right)=\frac{1}{\sqrt{2\pi }{\sigma }_M}\exp \left(-\frac{{\left(z-x^{\prime }{\beta }_M\right)}^2}{2{\sigma }_M^2}\right)$

将${\theta }_M$ 的Fiducial密度$r\left({\theta }_M|y\right)$ 限制到那些被容许的集合上，得

$r_{\epsilon }\left({\theta }_M|y\right)=r\left({\theta }_M|y\right)\cdot h\left({\beta }_M\right)$

则$z$ 的Fiducial预测密度(简记为FPD)为

$\begin{array}{c}\widehat{f}\left(z|y,M\right)={\displaystyle \int f\left(z|{\theta }_M,M\right)}{r}_{\epsilon }\left({\theta }_M|y\right)\text{d}{\theta }_M\\ ={\displaystyle {\int }_{{ℝ}^{\left|M\right|}}{\displaystyle {\int }_0^{+\infty }h\left({\beta }_M\right)\frac{{\pi }^{-\frac{\left|M\right|+1}{2}}RS{S}_M^{\frac{n-\left|M\right|}{2}}{\left|{X^{\prime }}_M{X}_M\right|}^{\frac{1}{2}}}{2^{\frac{n+1}{2}}\Gamma \left(\frac{n-\left|M\right|}{2}\right)}}}{\left({\sigma }_M^2\right)}^{-\frac{n+3}{2}}\exp \left\{-\frac{{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2}{2{\sigma }_M^2}\right\}\text{d}{\sigma }_M^2\text{d}{\beta }_M\end{array}$

为了简化起见，先计算下面的积分

${\displaystyle {\int }_0^{+\infty }{\left({\sigma }_M^2\right)}^{-\frac{n+3}{2}}\exp \left\{-\frac{{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2}{2{\sigma }_M^2}\right\}\text{d}{\sigma }_M^2}$ (3)

令$t=\frac{{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2}{2{\sigma }_M^2}$ ，则

$\left(3\right)式=2^{\frac{n+1}{2}}\Gamma \left(\frac{n+1}{2}\right){\left[{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2\right]}^{-\frac{n+1}{2}}$

因此

$\begin{array}{c}\widehat{f}\left(z|y,M\right)={\displaystyle {\int }_{{ℝ}^{\left|M\right|}}h\left({\beta }_M\right)\frac{RS{S}_M^{\frac{n-\left|M\right|}{2}}{\left|{X^{\prime }}_M{X}_M\right|}^{\frac{1}{2}}\Gamma \left(\frac{n+1}{2}\right)}{{\pi }^{\frac{\left|M\right|+1}{2}}\Gamma \left(\frac{n-\left|M\right|}{2}\right)}}{\left[{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2\right]}^{-\frac{n+1}{2}}\text{d}{\beta }_M\\ =\frac{\Gamma \left(\frac{n+1}{2}\right)}{\sqrt{\pi }\Gamma \left(\frac{n}{2}\right)}{\displaystyle {\int }_{{ℝ}^{\left|M\right|}}h\left({\beta }_M\right)\frac{1}{\sqrt{{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2}}}{\left[\frac{{\left(z-x^{\prime }{\beta }_M\right)}^2}{{\Vert y-X_M{\beta }_M\Vert }_2^2}+1\right]}^{-\frac{n}{2}}r\left({\beta }_M|y\right)\text{d}{\beta }_M\\ \propto E\left\{h\left({\beta }_M\right){\left({\Vert y-X_M{\beta }_M\Vert }_2^2\right)}^{-\frac{1}{2}}{\left[\frac{{\Vert y-X_M{\beta }_M\Vert }_2^2}{{\left(z-x^{\prime }{\beta }_M\right)}^2+{\Vert y-X_M{\beta }_M\Vert }_2^2}\right]}^{\frac{n+1}{2}}\right\}\end{array}$

其中最后一行是关于${\beta }_M$ 的期望。

Vehtari等人[26]给出了一种对数逐点预测密度期望的贝叶斯留一估计，我们基于此提出对数逐点Fiducial预测密度估计，即对于模型$M$ 和数据集$y$ ，

$elfp{d}_{loo}={\displaystyle \sum _{i=1}^n\log \widehat{f}\left(y_i|y_{-i},M\right)}\triangleq h\left(y,M\right)$

进而有

$h\left(y,M\right)\propto {\displaystyle \sum _{i=1}^{n}E\left\{h\left({\beta }_M\right){\left({\Vert y_{-i}-X_{M,-i}{\beta }_M\Vert }_2^2\right)}^{-\frac{1}{2}}{\left[\frac{{\Vert y_{-i}-X_{M,-i}{\beta }_M\Vert }_2^2}{{\Vert y-X_M{\beta }_M\Vert }_2^2}\right]}^{\frac{n+1}{2}}\right\}}$

假设所有模型的先验分布$\pi \left(M\right)$ 为均匀分布，由此我们定义一种新的模型不确定性度量

$P_r\left(M|y\right)=\frac{h\left(y,M\right)\pi \left(M\right)}{{\displaystyle \sum _{M^{\prime }\in S}h\left(y,M^{\prime }\right)\pi \left(M^{\prime }\right)}}\propto h\left(y,M\right)$

于是基于Fiducial预测密度的选择准则给出，即

$\widehat{M}=\underset{M\in S}{\arg \max }{P}_r\left(M|y\right)=\underset{M\in S}{\arg \max }h\left(y,M\right)$

将该模型选择方法成为Fiducial预测密度模型选择方法，简记为FPMC。

该方法也具有相合性，即在部分正则化条件下，真实模型$M_o$ 满足

$P_r\left(M|y\right)=\frac{h\left(y,M\right)\pi \left(M\right)}{{\displaystyle \sum _{M^{\prime }\in S}h\left(y,M^{\prime }\right)\pi \left(M^{\prime }\right)}}\propto h\left(y,M\right)$

其中$P_y$ 是指与观测值$Y$ 的抽样分布有关的概率度量。该性质的证明与William和Hannig [21]文中证明GFI概率下的相合性的思路近乎相同，所以证明过程此处不展开讨论。

3. FPMC的实现

仿照赵勇超[23]文中采用分组独立的Metropolis-Hastings算法实现FMC的思路，我们也采取该算法并将其应用进FPMC。基本思路是首先通过(2)式引入潜变量从而得到一个解析表达式。然后抽取一个模型$M$ ，并从$M$ 对应的多元$t$ 分布中抽取$N_{MC}$ 样本来计算期望的蒙特卡罗估计，记为$B$ 。经过上述过程可以得到一个联合的马尔可夫链$\left(M,B\right)$ ，通过仅关注$M$ 部分，便可从$M$ 的伪边际分布中获取样本。通过计算模型在稳定后的马尔可夫链中的后验概率，选择后验概率最大的模型$M$ 作为最优模型。算法的具体实现步骤见表1：

Table 1. Distribution-independent MH algorithm

表1. 分组独立的MH算法

输入：$y$ ，$X$ ，$N_{\text{step}}$ ，$N_{\text{burnin}}$ ，$N_{\text{MC}}$ ，一个初始的模型$M^0$ ，一个提议概率函数$q\left(\tilde{M}|M^0\right)$ ， 一个马尔科夫链$S_{\text{MCMC}}=\varphi$ 以及默认的$\epsilon$ 。

步骤一：执行如下循环，对于$k=1,\cdots ,N_{\text{step}}$ ， (1) 自$M^0$ 开始通过随机游走的方式搜索下一个模型$M\in S$ ，记为$\tilde{M}$ 。 $\tilde{M}=\{\begin{array}{l}{M}^0\cup \left\{一个新的协变量\right\}\hfill \\ M^0/\left\{一个已存在的协变量\right\}\hfill \\ M^0/\left\{一个已存在的协变量\right\}\cup \left\{一个新的协变量\right\}\hfill \end{array}\begin{array}{c}\\ \\ \end{array}\begin{array}{c}\omega =1/3\\ \omega =1/3\\ \omega =1/3\end{array}$ (2) 通过(4) 对$\left({\beta }_{M^0},{\sigma }_{M^0}^2\right)$ 和$\left({\beta }_{\tilde{M}},{\sigma }_{\tilde{M}}^2\right)$ 进行估计，分别记为$\left({\widehat{\beta }}_{M^0,-i},{\widehat{\sigma }}_{M^0,-i}^2\right)$ 和$\left({\widehat{\beta }}_{\tilde{M},-i},{\widehat{\sigma }}_{\tilde{M},-i}^2\right)$ 。 (3) 分别从分布

$t_{n-\left|M\right|}\left({\widehat{\beta }}_M,{\widehat{\sigma }}_M^2{\left({X^{\prime }}_M{X}_M\right)}^{-1}{X^{\prime }}_{M}Y\right)$ ，$M\in \left\{M^0,\tilde{M}\right\}$ 抽取样本${\beta }_{M^0}^j$ 和${\beta }_{\tilde{M}}^j$ ，$j=1,2,\cdots ,N_{MC}$ 。 (4) 对于所有$j=1,2,\cdots ,N_{MC}$ ，计算 ${{\rm H}}_M^j\triangleq {\displaystyle \sum _{i=1}^{n}h\left({\beta }_M^j\right){\left({\Vert y_{-i}-X_{M,-i}{\beta }_M^j\Vert }_2^2\right)}^{-\frac{1}{2}}{\left[\frac{{\Vert y_{-i}-X_{M,-i}{\beta }_M^j\Vert }_2^2}{{\Vert y-X_M{\beta }_M^j\Vert }_2^2}\right]}^{\frac{n+1}{2}}}$ ，$M\in \left\{M^0,\tilde{M}\right\}$ 其中${\beta }_M^j$ 可由显式$L_0$ 最小化算法计算。 (5) 计算(17)式的蒙特卡洛估计 ${\widehat{P}}_r\left(M|y\right)\propto \frac{1}{N_{MC}}{\displaystyle \sum _{j=1}^{N_{MC}}{{\rm H}}_M^j}$ ，$M\in \left\{M^0,\tilde{M}\right\}$

(6) 如果$\min \left\{\frac{{\widehat{P}}_r\left(\tilde{M}|y\right)q\left(M^0|\tilde{M}\right)}{{\widehat{P}}_r\left(M^0|y\right)q\left(\tilde{M}|M^0\right)},1\right\}>u^{\ast },u^{\ast }\sim U\left(0,1\right)$ ， 那么$M^0=\tilde{M}$ 并且$S_{\text{MCMC}}=\left\{S_{\text{MCMC}},\tilde{M}\right\}$； 否则$M^0=M^0$ 并且$S_{\text{MCMC}}=\left\{S_{\text{MCMC}},\varphi \right\}$ 。 终止循环。

步骤二：舍弃马尔可夫链$S_{\text{MCMC}}$ 上前$N_{\text{burnin}}$ 个样本。计算所有模型在剩余$N_{\text{step}}-N_{\text{burnin}}$ 个样本上的后验概率。选择具有最高后验概率的模型作为最优模型$M_{\text{best}}$ 。

步骤三：输出最优模型$M_{\text{best}}$ 。

4. 数值模拟

本节对提出的FPMC进行数值模拟，并将其模拟结果与Lasso和EN进行比较。比较的指标主要包括两方面：一方面是变量识别的准确性，主要有积极变量识别率的平均(AEIR)；消极变量识别率的平均(IEIR)；稀疏度估计误差(SLE)。另一方面是预测精度，用到的是均方根误差(RMSE)。设$M_o$ 和$M_o^c$ 为真实模型$M_o$ 中包含的积极效应和消极效应所构成的集合，$M_i$ 为第$i$ 次试验中选择的最优模型，$N$ 为重复实验次数。则各个指标的定义如下

$SLE=\frac{1}{N}{\displaystyle \sum _{i=1}^N\left|M_i\right|-\left|M_o\right|}$ , $RMSE=\sqrt{\frac{1}{n_{test}}{\displaystyle \sum _{i=1}^{n_{test}}{\left(y_i-{\widehat{y}}_i\right)}^2}}$

$AEIR=\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{\left|M_i\cap M_o\right|}{\left|M_o\right|}}\times 100\%$ , $IEIR=\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{\left|M_i\cap M_o^c\right|}{\left|M_o^c\right|}}\times 100\%$

AEIR越大，说明所选模型中包含的积极变量越多，故AEIR越大越好(越接近1越好)；IEIR，RMSE越小越好，SLE越接近0越好。

4.1. 模拟1

本例与Williams和Hannig [21]的研究中的模拟设定1保持一致。从(1)式中生成数据，其中设计阵$X_{n\times p}$ 服从多元正态分布$N_p\left(0_p,\Sigma \right)$ 。真实模型的参数向量为

${\beta }_{M_o}=\left\{-1.5,-1,-0.8,-0.6,0.6,0.8,1,1.5,0_{p-8}\right\}$

其中$p\in \left\{40,80,120\right\}$ 。$\Sigma \left(={\left({\sigma }_{ij}\right)}_{p\times p}\right)$ 被设定为单位矩阵，训练集大小为$n=100$ ，测试集大小为$n_{\text{test}}=30$ ，误差项的方差${\sigma }^2$ 固定为1。在赵勇超[23]的模拟1中已通过ABESS算法得出了此时$p_o$ 的估计值取$p_o=8$ 最合适。下面的表2、表3和表4给出了FPMC、LASSO和EN的表现。

Table 2. The performances of three methods when p = 40

表2. p = 40时三种方法的表现

	SLE	RMSE	AEIR	IEIR
FPMC	1.19	3.6589	0.936	0.053
Lasso	0.53	3.6743	0.828	0.0597
EN	1.23	3.6740	0.875	0.0698

Table 3. The performances of three methods when p = 80

表3. p = 80时三种方法的表现

	SLE	RMSE	AEIR	IEIR
FPMC	1.10	3.7422	0.930	0.023
Lasso	0.62	3.7501	0.821	0.028
EN	1.52	3.7592	0.884	0.034

Table 4. The performances of three methods when p = 120

表4. p = 120时三种方法的表现

	SLE	RMSE	AEIR	IEIR
FPMC	0.95	3.6385	0.927	0.014
Lasso	1.00	3.6736	0.845	0.020
EN	1.87	3.6772	0.896	0.024

模拟1结果表明，无论变量个数p取40，80还是120，FPMC的RMSE总是最小的，即无论p小于n还是大于n，我们方法的预测准确度都是最高的。随着p值的增加新方法的SLE数值越来越接近0，并且当p大于n时FPMC是三种方法SLE值最小的。从被选模型的识别效果和可解释性方面来看，FPMC的AEIR值在三种情形下都是最大值、IEIR值都是最小值，即FPMC方法能更精确地识别积极效应、更精确地压缩消极效应。综合来说，本文提出的FPMC在预测和识别方面都更胜一筹。

4.2. 模拟2

采用与Williams和Hannig [21]的研究中模拟设定2相同的数据生成方式，其中$n=30$ 。数据生成过程为

$Y\sim N_n\left(1\cdot x^{\left(1\right)}+1\cdot x^{(2)}+\cdots +1\cdot x^{\left(9\right)},I_n\right)$

其中$x^{\left(1\right)},x^{(2)},x^{(3)}\stackrel{i.i.d}{~}{N}_n\left(0,I_n\right)$ ，$x^{\left(4\right)}\sim N_n\left(0.25\cdot x^{\left(1\right)},{0.1}^2{I}_n\right)$ ，$x^{\left(5\right)}\sim N_n\left(0.50\cdot x^{\left(2\right)},{0.1}^2{I}_n\right)$ ，$x^{\left(6\right)}\sim N_n\left(-0.75\cdot x^{\left(3\right)},{0.1}^2{I}_n\right)$ ，$x^{\left(7\right)}\sim N_n\left(1\cdot x^{\left(1\right)}+1\cdot x^{\left(3\right)},{0.1}^2{I}_n\right)$ ，$x^{\left(8\right)}\sim N_n\left(1\cdot x^{\left(2\right)}-1\cdot x^{\left(3\right)},{0.1}^2{I}_n\right)$ ，$x^{\left(9\right)}\sim N_n\left(1\cdot x^{\left(1\right)}+1\cdot x^{\left(2\right)}+1\cdot x^{\left(3\right)},{0.1}^2{I}_n\right)$ 。

由此我们建立线性模型，经过10折交叉验证我们取FPMC方法中的$p_o$ 估计值为6。进行1000次重复模拟实验，三种方法RMSE和模型长度的数值如下：

Table 5. The performances of three methods in Simulation 2

表5. 模拟2中三种方法的表现

	FPMC	Lasso	EN
RMSE	4.6854	5.3463	5.6291
Length of model	6.1	7.0	6.9

由表5可知FPMC方法的RMSE数值最小，即它的预测准确度最高。而且它在保证准确度的情况下没有过多的选择变量，控制了所选模型的复杂程度。

5. 实例分析

5.1. 实例一

活塞拍击的声音是一种由活塞二次运动引起的发动机噪音。为了降低活塞拍击的噪音，选取了6种对噪音影响较大的因素进行分析，希望通过改变这6种因素达到减少噪音的目的。6种因素分别为活塞和气缸套之间的设定间隙$x_1$ ，峰值压力位置$x_2$ ，裙部长度$x_3$ ，裙部轮廓形状$x_4$ ，裙部椭圆度$x_5$ ，活塞销偏置$x_6$ 。数据集来源于Huang等[27]并进行了归一化处理，共包含100个观测样本，每个样本有6个输入变量，模型中可能包含所有的线性主效应、二次主效应以及正交多项式编码下的所有双因素相互作用，因此共有72个基变量。我们从100个样本中取80个作为训练集，剩下的20个做测试集，模拟结果在表6中给出。

Table 6. Data simulation results of a piston slap noise example

表6. 活塞拍击噪声实例的数据模拟结果

	FPMC	Lasso	EN
RMSE	0.5088	0.6462	0.6368
Length of model	16	39	32

由表6可以得出，FPMC所选的模型长度不仅比Lasso和EN的要小，RMSE数值也是三种方法中最小的。总体来看，FPMC方法能够在有效简化模型同时具有较好的预测效果。

5.2. 实例2

关于影响维生素B生产速率的基因位点识别，是一个典型的基于超高维线性回归模型的模型选择问题。该数据集首次公开于Buhlmann等人的研究，其中包含71个观测值和4088个基因表达，其中响应变量为枯草芽孢杆菌的维生素B生产速率的对数。我们采用线性回归模型来探索数据集中响应变量与协变量之间的关系，因此该数据集符合模型2-(1)的设定。首先用ABESS算法得到该真实数据集的$p_o$ 的估计值为${\widehat{p}}_o=6$ 。然后用BESS (Best Subset Selection)算法进行初筛并保留$p=71$ 个基因表达进行分析。需要注意的是，BESS和ABESS算法都源自Zhu等人[25]的研究。区别在于前者需要指定稀疏度(在这里稀疏度指定为$p=71$ )，而后者是依靠数据驱动的自适应版本。对于最终保留的71个基因表达，分别采用FPMC、Lasso、EN来进行变量选择，基本的设置与第四节保持一致。

Table 7. Data simulation results of a gene recognition of vitamin B example

表7. 维生素B基因识别实例的数据模拟结果

	FPMC	Lasso	EN
RMSE	0.292	0.331	0.324
Length of model	8	10	9

由表7可以得出，FPMC所选的模型长度不仅比Lasso和EN的要小，RMSE数值也是三种方法中最小的。总体来看，FPMC方法表现更好。

6. 结论

本文提出了高斯回归模型的Fiducial预测模型选择方法FPMC，并同时给出它的理论步骤与实践算法。之后其与Lasso和EN相比较，通过两个模拟和两个实例我们得到三种方法中最好的是FPMC，不仅预测误差总是最小的，而且模型拟合精度和复杂度都表现良好。