多任务高斯过程回归的模型平均

Model Averaging for Multi-Task Gaussian Process Regression

摘要: 模型平均因其稳健性好、预测精度高等优点在当代统计学和计量经济学界有广泛应用。多任务高斯过程回归同时考虑了输入变量之间和输出变量之间的相关性，能够得到更为准确的预测值。考虑到输出之间的相关性，本文先构建多任务高斯过程回归模型(MGPR)，再进行模型平均及选择；不考虑输出之间的相关性时，分别构建单任务高斯过程回归(GPR)，再进行模型平均及选择。数值模拟结果表明，输出噪声为同方差时，MGPR的模型平均比GPR具有更小的损失风险，且MGPR的计算效率相对较高。输出噪声同方差时，采用中国少数民族地区城市化率、非农业人口比例及受教育年限的实际数据；输出噪声异方差时，采用苄醇转化率和转频率的实际数据，均表明MGPR的模型平均具有更小的损失风险。

Abstract: Model averaging, owing to its advantages of robustness and high predictive accuracy, finds exten-sive application in contemporary statistics and econometrics. Multi-task Gaussian process regres-sion, considering both the interdependencies among input variables and the correlations among output variables, enables more precise predictions. Taking into account the inter-output correla-tions, this paper constructs a Multi-Task Gaussian Process Regression model (MGPR) followed by model averaging and selection. Conversely, when disregarding inter-output correlations, single-task Gaussian Process Regression (GPR) models are individually constructed, followed by the same pro-cess of model averaging and selection. Numerical simulation results demonstrate that when the output noise exhibits homoscedastic variance, MGPR’s model averaging exhibits lower loss risk compared to GPR, with relatively higher computational efficiency. For the case of homoscedastic output noise, actual data on urbanization rates in China’s ethnic minority regions, non- agricultural population ratios, and years of education are employed. In the case of heteroscedastic output noise, real data on benzyl alcohol conversion rates and turnover frequencies are utilized, both indicating that MGPR’s model averaging yields lower loss risk.

1. 引言

多任务高斯过程能够处理多个相关输出变量，提高模型的灵活性和准确性。其通过协方差函数描述输出间的相关性，并共享相关性结构来获取信息。然而，单一的多任务高斯过程模型无法覆盖所有数据特征，因此引入模型平均方法可以减少模型的不确定性，提高预测性能。模型平均方法通过组合多个基本模型，综合优点并改善预测稳定性。这种方法在多个领域得到广泛应用，为解决复杂问题提供了强大而可靠的工具。

当前学术研究主要集中在贝叶斯模型平均和频率模型平均这两个方向。Zhang J.和Taflanidis A. A. [1] 提出一种高效的BMA克里金回归模型，通过将全局趋势的回归与提供局部校正的高斯过程近似相结合来制定预测。除了贝叶斯模型平均，学术界还关注了频率模型平均方法。权重计算准则是频率模型平均方法的关键问题。Wang Miaomiao等人 [2] 提到基于AIC和BIC的平滑AIC (S-AIC)和平滑BIC (S-BIC)准则，通过极小化这些信息准则来确定权重。Priyam Mitra等人 [3] 在没有对候选模型集设置限制的情况下建立了通用的频率模型平均框架，并通过最小化均方误差来获得自适应权重。Zhongqi Liang等人 [4] 提到通过极小化Mallows准则来确定模型平均的权重。并证明在宽松条件下，利用最小化Mallows准则所确定的权重具有渐近最优性。Shaobo Jin和Sebastian Ankargren [5] 将FMA方法应用于连续数据的结构方程模型，并提出了有效的置信区间和 ${\chi }^2$ 检验统计量。Carroll R.等人 [6] 提出简约模型平均法，将权重1分配给误差最小的准确模型，使其他准确模型的权重为0。Zhang X.等人 [7] 探讨MMA和JMA方法对参数模型平均估计的分布特性。

多任务高斯过程在统计学和概率论中占据重要地位。Zexun Chen等人 [8] 提供了对多任务高斯过程的精确定义，并介绍了其基本属性。Chen Z.等人 [9] 介绍多任务高斯过程的关键属性，如平稳性和独立性。多任务高斯过程近几年在各个领域应用逐渐广泛，并在多个方面不断改进，以提高预测精度。Xiaodan Hong等人 [10] 通过多任务高斯过程回归和依赖性多任务高斯过程回归，考虑了相关噪声的影响并提升了建模性能。R. Kontar等人 [11] 提出了一种可扩展和正则化的方法，以减小多任务高斯过程中的负迁移。J. Zapata等人 [12] 通过引入部分可分性概念和函数图模型，实现了对多任务高斯过程的建模。Rong Zhu等人 [13] 在本文中提及了如何在多个相关因变量的情况下，用模型平均与选择的方法对建立的模型提高预测精度。Zexun Chen等人 [14] 提及高斯过程中如何设置超参数。Pengfei Wei等人 [15] 在结构系统可靠性的研究中，基于多任务输出高斯过程的主动学习来提高系统可靠性。ChenZ等人 [16] 提出MV-TPR，并与MV-GPR作比较，证明在所考虑的数据集上的MV-TPR的有效性。MickaB等人 [17] 提出在高斯过程中用潜在变量解决异方差问题，并运用大量实例证明自己的方法最优。Hansen等人 [18] 提出JMA方法，通过删除一个样本的交叉验证法最小化误差向量。Tang Qinghu等人 [19] 设计了两个输出任务相关的试验。

本文的其余部分结构如下：在第2节中，我们介绍多任务高斯过程回归模型、核、超参数以及五种模型选择及模型平均方法；在第3节中，通过数值模拟构建MGPR及GPR，并分别进行模型平均及选择，最后将以上方法运用到第六次人口普查中部分地区数据及苄醇转化率和转频率的统计研究中；第4节是对文章的总结和讨论。

本文在以上文献的基础上，考虑现实情况中多个输出任务具有相关性，将模型平均引入多任务高斯过程回归中，这是本文的创新点。

2. 多任务高斯过程及相关介绍

2.1. 多任务高斯过程模型

如果f是X上的多任务高斯过程，具有向量型均值函数 $\mu :X\mapsto R^d$ ，协方差函数(也称为核) $k:X\times X\mapsto R$ 且半正定参数矩阵 $\Omega \in R^{d\times d}$ ，那么 $f~MGP\left(\mu ,k,\Omega \right)$ 。并且在此文中，考虑将噪声项包含在核函数内的模型。

给定n对样本观测值 ${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^n,x_i\in R^p,y_i\in R^{1\times d}$ ，我们假设模型为：

$f~MGP\left(\mu ,k\text{'},\Omega \right)$ (1)

$y_i=f\left(x_i\right),i=1,\cdots ,n$ (2)

其中 $k\text{'}=k\left(x_i,x_j\right)+{\delta }_{ij}{\sigma }_n^2$ ，且当 $i=j$ 时， ${\delta }_{ij}=1$ ，否则 ${\delta }_{ij}=0$ ，上式中的第二项代表的是随机噪声项。

在 $X_{\ast }={\left[x_{n+1},\cdots ,x_{n+m}\right]}^{\text{T}}$ 处的预测变量为 $f_{\ast }={\left[f_{\ast }{}_1,\cdots ,f_{\ast }{}_m\right]}^{\text{T}}$ ,训练集 $Y={\left[y_1^{\text{T}},\cdots ,y_n^{\text{T}}\right]}^{\text{T}}$ 和预测变量 $f_{\ast }$ 的联合分为

$\left[\begin{array}{c}Y\\ f_{\ast }\end{array}\right]~MGP\left(0,\left[\begin{array}{cc}K\text{'}\left(X,X\right)& K\text{'}\left(X_{\ast },X\right)\\ K\text{'}\left(X_{\ast },X\right)& K\text{'}\left(X_{\ast },X_{\ast }\right)\end{array}\right],\Omega \right)$ (3)

其中 $K\text{'}\left(X,X\right)$ 是一个 $n\times n$ 矩阵，第 $\left(i,j\right)$ 个元素是 $K\text{'}{\left(X,X\right)}_{ij}=k\text{'}\left(x_i,x_j\right)$ ，且 $K\text{'}\left(X_{\ast },X\right)$ 是 $m\times n$ 矩阵，第 $\left(i,j\right)$ 个元素是 $K\text{'}{\left(X_{\ast },X\right)}_{ij}=k\text{'}\left(x_{n+i},x_j\right)$ ， $K\text{'}\left(X_{\ast },X_{\ast }\right)$ 是 $m\times m$ 矩阵，第 $\left(i,j\right)$ 个元素是 $K\text{'}{\left(X_{\ast },X_{\ast }\right)}_{ij}=k\text{'}\left(x_{n+i},x_{n+j}\right)$ 。因此，根据多任务高斯过程的条件分布，其预测分布为 $p\left(f_{\ast }|X,Y,X_{\ast }\right)~MGP\left(\hat{\mu },\hat{\Sigma },\hat{\Omega }\right)$ ，其中

$\hat{\mu }=K\text{'}{\left(X_{\ast },X\right)}^{\text{T}}K\text{'}{\left(X,X\right)}^{-1}Y$ (4)

$\hat{\Sigma }=K\text{'}\left(X_{\ast },X_{\ast }\right)-K\text{'}{\left(X_{\ast },X\right)}^{\text{T}}K\text{'}{\left(X,X\right)}^{-1}K\text{'}\left(X_{\ast },X\right)$ (5)

$\hat{\Omega }=\Omega$ (6)

此外，期望和协方差如下所示：

$E\left|f_{\ast }\right|=\hat{\mu }=K\text{'}{\left(X_{\ast },X\right)}^{\text{T}}K\text{'}{\left(X,X\right)}^{-1}Y$ (7)

$cov\left(vec\left(f_{\ast }^{\text{T}}\right)\right)=\hat{\Sigma }\otimes \hat{\Omega }=\left[K\text{'}\left(X_{\ast },X_{\ast }\right)-K\text{'}{\left(X_{\ast },X\right)}^{\text{T}}K\text{'}{\left(X,X\right)}^{-1}\right]\otimes \Omega$ (8)

2.2. 核介绍

上述回归模型中有两个协方差阵，列协方差阵k与行协方差阵 $\Omega$ ，但只有列协方差取决于输入，并且被视为内核，因为它包含有关我们希望学习和定义数据点之间接近度和相似性的函数的假设，行协方差阵取决于输出，构建多任务高斯过程回归，既要考虑到列协方差在又要考虑到行协方差阵。与传统GPR一样，核的选择对多任务高斯过程回归的性能有深远的影响。现有方法中已经提出了许多有用的核，如线性核、有理二次核和matern。内核还可以通过ARD定义：

$k_{SEard}=s_f^2\exp \left(-\frac{{\left(x-x\text{'}\right)}^{\text{T}}{\Theta }^{-1}\left(x-x\text{'}\right)}{2}\right)$ (9)

其中 $s_f^2$ 指的是信号方差，也可以认为是输出尺度的振幅， $\Theta$ 是一个对角矩阵，具有组成元素为 ${\left\{{\mathcal{l}}_i^2\right\}}_{i=1}^p$ ，它表示每个相应输入维度的长度尺度，参数 ${\mathcal{l}}_i$ 是第i个输入的尺度。

2.3. 超参数介绍

参数估计需要从训练数据中估计核所涉及的超参数和行协方差矩阵。传统GP模型中使用的许多方法，例如最大似然估计(MLE)、最大后验(MAP)和马尔可夫链蒙特卡罗(MCMC)，可用于我们提出的模型。尽管蒙特卡罗方法可以执行GPR，而不需要估计超参数，但由于蒙特卡罗方法的计算成本较高，常用的方法是通过MLE来估计它们。因此，我们考虑使用MLE进行参数估计。与传统的GPR模型相比，行协方差阵 $\Omega$ 是一个额外的参数。因此，未知参数包括核中的参数，噪声方差 ${\sigma }_n^2$ 及行协方差参数矩阵 $\Omega$ 。

因为 $\Omega$ 是半正定矩阵，可以被分解为 $\Omega =\Phi {\Phi }^{\text{T}}$ ，其中

$\Phi =\left[\begin{array}{cccc}{\varphi }_{11}& 0& \cdots & 0\\ {\varphi }_{21}& {\varphi }_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\varphi }_{d1}& {\varphi }_{d2}& \cdots & {\varphi }_{dd}\end{array}\right]$ (10)

为了保证 $\Phi$ 的唯一性，对角元素被限制为正，且令 ${\phi }_{ii}=\ln \left({\varphi }_{ii}\right)$ ， $i=1,2,\cdots ,d$ 。

在MV-GPR模型中， $Y~MGP\left(0,K\text{'},\Omega \right)$ ，其中 $K\text{'}$ 是带有噪声项的列协方差矩阵，有 $K{\text{'}}_{ij}=k\text{'}\left(x_i,x_j\right)$ ，即为 $K\text{'}=K+{\sigma }_n^{2}I$ ,其中K是不带噪声项的列协方差阵，有 $K_{ij}=k\left(x_i,x_j\right)$ 。内核中的参数表示为 $K=K_{\theta }$ ，超参数集为 $\Theta =\left\{{\theta }_1,{\theta }_2,\cdots \right\}$ ，因此，

$\frac{\partial K\text{'}}{\partial {\sigma }_n^2}=I_n,\frac{\partial K\text{'}}{\partial {\theta }_i}=\frac{\partial K}{\partial {\theta }_i}$ (11)

其似然函数可表示为 $L=\frac{nd}{2}\ln \left(2\pi \right)+\frac{d}{2}\ln \det \left(K\text{'}\right)+\frac{n}{2}\ln \det \left(\Omega \right)+\frac{1}{2}tr\left({\left(K\text{'}\right)}^{-1}Y{\Omega }^{-1}{Y}^{\text{T}}\right)$ 。

利用似然函数对超参数进行估计有：

$\frac{\partial L}{\partial {\sigma }_n^2}=\frac{d}{2}tr\left({\left(K\text{'}\right)}^{-1}\right)-\frac{1}{2}tr\left({\alpha }_{K\text{'}}{\Omega }^{-1}{\alpha }_{K\text{'}}^{\text{T}}\right)$ (12)

$\frac{\partial L}{\partial {\theta }_i}=\frac{d}{2}tr\left({\left(K\text{'}\right)}^{-1}\frac{\partial K_{\theta }}{\partial {\theta }_i}\right)-\frac{1}{2}tr\left({\alpha }_{K\text{'}}{\Omega }^{-1}{\alpha }_{K\text{'}}^{\text{T}}\frac{\partial K_{\theta }}{\partial {\theta }_i}\right)$ (13)

$\frac{\partial L}{\partial {\varphi }_{ij}}=\frac{n}{2}tr\left[{\Omega }^{-1}\left(E_{ij}{\Phi }^{\text{T}}+\Phi E_{ji}\right)\right]-\frac{1}{2}tr\left({\alpha }_{\Omega }{\left(K\text{'}\right)}^{-1}{\alpha }_{\Omega }^{\text{T}}\left(E_{ij}{\Phi }^{\text{T}}+\Phi E_{ji}\right)\right)$ (14)

$\frac{\partial L}{\partial {\phi }_{ii}}=\frac{n}{2}tr\left[{\Omega }^{-1}\left(J_{ii}{\Phi }^{\text{T}}+\Phi J_{ii}\right)\right]-\frac{1}{2}tr\left({\alpha }_{\Omega }{\left(K\text{'}\right)}^{-1}{\alpha }_{\Omega }^{\text{T}}\left(J_{ii}{\Phi }^{\text{T}}+\Phi J_{ii}\right)\right)$ (15)

其中 ${\alpha }_{K\text{'}}={\left(K\text{'}\right)}^{-1}Y$ ， ${\alpha }_{\Omega }={\Omega }^{-1}{Y}^{\text{T}}$ ， $E_{ij}$ 是 $d\times d$ 单位阵， $J_{ii}$ 与 $E_{ij}$ 形式相同，但对角线上元素为 ${\text{e}}^{{\phi }_{ii}}$ 。因此，采用共轭梯度法用来最小化负对数似然函数，得到参数估计。且由于随机噪声被纳入核函数，噪声方差与其他超参数一起估计。

2.4. 模型选择及模型平均方法

实际情况中，输入变量之间及输出变量之间存在相关性。考虑输出任务之间不存在相关性时，构建GPR的模型平均；考虑输出任务之间的相关性时，构建MGPR的模型平均，我们想得到哪种模型的计算效率及预测准确率较高，故我们将这两种模型进行比较。

线性回归模型中应用较多的是以下几种传统的模型选择和模型平均的方法，本文为了研究分析各种模型平均及选择方法在MGPR及GPR下的预测结果，考虑了以下模型选择及模型平均方法：

AICc: $AI{C}_m=n\log \left({\sigma }_m^2\right)+2n\left(k_m+p+1\right)/\left(n-k_m-p-2\right)$ (16)

其中 ${\sigma }_m^2=\frac{1}{n}{\left(Y-F_m{\hat{\beta }}_m\right)}^{\prime }{\hat{C}}^{-1}\left(Y-F_m{\hat{\beta }}_m\right)$ ，P为与自相关函数相关的参数个数，C为单任务输出或者多任务输出时的协方差矩阵。具有最小AICc值挑选出来作为最优模型。

(2) $w_m=\exp \left(-AI{C}_m/2\right)/{\displaystyle \underset{i=1}{\overset{M}{\sum }}\exp \left(-AI{C}_i/2\right)}$ ，其中 $AI{C}_m=n\log \left({\hat{\sigma }}_m^2\right)+2\left(k_m+p+1\right)$ (17)

(3) $w_m=\exp \left(-BI{C}_m/2\right)/{\displaystyle \underset{i=1}{\overset{M}{\sum }}\exp \left(-BI{C}_i/2\right)}$ ，其中 $BI{C}_m=n\log \left({\hat{\sigma }}_m^2\right)+\log \left(n\right)\left(k_m+p+1\right)$ (18)

(4) MMA：采用同方差线性回归模型提出的模型平均方法。MMA方法通过极小化Mallows准则得到权重：

$w=\arg \min C_n\left(w\right)$ (19)

$C_n\left(w\right)={\Vert A\left(w\right)y\Vert }^2+2{\sigma }^{2}tr\left(P\left(w\right)\right)$ (20)

其中 $A_k=I_k-P_k$ ， $A\left(w\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}{w}_i{A}_i}$ ， $P_k=X_k{\left({X^{\prime }}_k{X}_k\right)}^{-1}{X}_k$ 为模型k的投影矩阵。

(5) MMMA：

$w=\arg \min C_n\left(w\right)$ (21)

$C_n\left(w\right)={\left\{Vec\left(Y\right)-Vec\left(\hat{\mu }\left(w\right)\right)\right\}}^{\prime }{\hat{C}}^{-1}\left\{Vec\left(Y\right)-Vec\left(\hat{\mu }\left(w\right)\right)\right\}$ (22)

其中 ${\hat{C}}^{-1}$ 指的是多个任务构建多任务高斯过程回归得到的列协方差矩阵。

(6) JMA：在异方差情形下提出的模型平均方法。

$w=\arg \min C{V}_n\left(w\right)$ (23)

(24)

记 $\tilde{u}\left(w\right)={\displaystyle \underset{m=1}{\overset{M}{\sum }}{w}^m{\tilde{u}}^m}$ ， ${\hat{e}}^m=y-{\hat{u}}^m$ ， $\tilde{e}=\left({\tilde{e}}^1,{\tilde{e}}^2,\cdots ,{\tilde{e}}^M\right)$ ，令 $P^m=x^m{\left({\left(x^m\right)}^{\prime }{x}^m\right)}^{-1}{\left(x^m\right)}^{\prime }$ ， $h_{ii}^m=x_i^m{\left({\left(x^m\right)}^{\prime }{x}^m\right)}^{-1}{\left(x_i^m\right)}^{\prime }$ 为 $P^m$ 中的第i个对角元素，令对角矩阵 $D_m={\left(1-h_{ii}^m\right)}^{-1}$ ，可得 ${\tilde{e}}^m=D_m{\hat{e}}^m$ 。其中， ${\tilde{u}}^m={\left({\tilde{u}}_1^m,{\tilde{u}}_2^m,\cdots ,{\tilde{u}}_n^m\right)}^{\prime }$ ， ${\tilde{u}}_i^m$ 是从样本中删除第i组观测值 $\left(x_i,y_i\right)$ 用最小二乘估计得到的 ${\hat{u}}^m$ 估计值，有 ${\hat{u}}^m=P^{m}y$ 。

(7) MT-JMA：将JMA扩展到多任务上的模型平均方法：

$w=\arg \min C_n\left(w\right)$ (25)

$C_n\left(w\right)={\left\{Vec\left(Y\right)-Vec\left(\tilde{u}\left(w\right)\right)\right\}}^{\prime }{\hat{C}}^{-1}\left\{Vec\left(Y\right)-Vec\left(\tilde{u}\left(w\right)\right)\right\}$ (26)

其中 ${\hat{C}}^{-1}$ 同MMMA方法指代的一样。 $\tilde{u}\left(w\right)$ 为JMA方法计算步骤相同所得到的估计值。

为了评估估计量，我们计算损失风险：

$Risk=\sqrt{\frac{1}{n}{\left(y-\hat{y}\right)}^2}$ (27)

其中y是真实值， $\hat{y}$ 是多任务高斯过程回归得到的候选模型进行模型平均及模型选择后的预测值。

3. 数值模拟

输出噪声同方差时，考虑输出间的相关性，即采用MGPR，通过数值模拟比较函数在aicc、MMMA、S-AIC、S-BIC这四种方法下的损失风险。在不考虑输出间的相关性时，分别构建GPR时，通过数值模拟比较函数在aicc、MMA、S-AIC、S-BIC这四种方法下的损失风险。并将MGPR与GPR的进行不同模型平均及选择方法后得到的损失风险进行比较。

3.1. 参数设置

输出噪声同方差时，为了衡量多任务与单任务在预测值精度这个方面的性能，我们考虑如下的一个已知函数模型。这个已知函数模型是定义在12维(p = 12)的输入空间[0,1]¹²上的，模型中的前六个变量 $x_1,x_2,\cdots ,x_6$ 对计算机试验输出结果的影响逐渐减弱，后六个变量的系数为零，即后六个变量为与输出没有相关性的变量。真模型表达成以下形式：

$y_i\left(x\right)=0.4x_1+0.3x_2+0.2x_3+0.1x_4+0.05x_5+0.01x_6+{\epsilon }_i\left(x\right)$ , $i=1,2$ (28)

其中 $\left[\epsilon \left(1\right),\epsilon \left(2\right)\right]~MGP\left(0,k_{SE},\Omega \right)$ 。 $k_{SE}$ 中的超参数为：

$\left[{\mathcal{l}}_i,s_f^2\right]=\left[1.001,1.001,1.001,1.001,1.001,1.001,1.001,1.001,1.001,1.001,1.001,1.001,0.1\right]$ .

且 $\Omega =\left(\begin{array}{cc}1& 0.25\\ 0.25& 1\end{array}\right)$ ，响应值y₁、y₂分别通过上式产生。在Matlab中，通过拉丁超立方抽样生成维数p = 12，样本量分别为N = 120，140，160，180的样本X，我们假设候选模型是嵌套的。即第m个模型包含前m个回归项。

使用损失风险来衡量多任务高斯过程模型平均后的平均模型对于随机产生的测试集G的预测精度，损失风险越小，预测精度越高。测试集G是使用计算机产生的1000个样本点。重复500次实验，计算平均的Risk。Risk越小时，我们认为MGPR进行模型平均比单个任务的GPR进行模型平均更有优势。

3.2. 试验结果

使用aicc、MMMA、MMA、S-AIC、S-BIC方法进行分析，仿真结果见表1~4：

Table 1. Risk of model averaging and selection (×10)

表1. 模型平均及选择的损失风险(×10)

Table 2. Risk of model averaging and selection (×10)

表2. 模型平均及选择的损失风险(×10)

Table 3. Risk of model averaging and selection (×10)

表3. 模型平均及选择的损失风险(×10)

Table 4. Risk of model averaging and selection (×10)

表4. 模型平均及选择的损失风险(×10)

(1) 将MGPR与GPR分别进行模型平均及模型选择，可以得到结论：不论样本量大小是多少，在两种情况中都能够进行比较的模型平均方法中，多个任务同时输出的表现要优于单个任务分别输出。因此可以得知，在已经考虑到输入变量之间的相关性的前提下，再考虑输出之间相关性时，可以提高预测结果的准确性。

(2) 由上述所有表可以得出，先构建MGPR，再进行模型选择及模型平均得到的Risk中，大部分情况下，S-AIC与S-BIC的损失风险十分接近。

(3) 在两个任务分别构建GPR后，再进行模型选择及模型平均得到的结果中，可以得出MMA方法的表现要差于其他模型平均方法，这是因为MMA方法未考虑输入变量之间的相关性，而其他模型平均方法考虑到输入变量间的相关性，这表明，考虑输入变量之间的相关性可以降低预测的损失风险。

3.3. 时间比较

本节在不同样本量下，分别将两个任务构建MGPR及GPR，将得到的候选模型采用上述提到的五种模型选择及模型平均方法，同时输出两个任务的预测风险，对以上过程重复500次，得到两个过程所耗费的时间。将以上不同情况的耗费时间进行比较，得到结果见表5：

Table 5. Time spent (minute)

表5. 耗费时间(分)

由表5可知，在不同样本量及不同模型平均及选择的算法下，构建MGPR所耗费的时间比构建GPR所耗费的时间要短。因此，结合上一节结论可知，在多个输出具有相关性的前提下，选择MGPR构建模型，并进行模型平均，可以降低损失风险，提高预测精度。

3.4. 实证分析

3.4.1. 输出噪声同方差

本节采用噪声同方差下提及的模型选择及模型平均方法，分析中国77个少数民族地区的数据，收集2010年中国第六次全国人口普查的数据。我们将城市化率、非农业人口占比以及受教育年限视为三个输出变量，我们考虑了6个与输出变量相关的输入变量:少数民族人口占比、城市和农村收入比的比例、固定资产人均投资、人均GDP、第二产业增加值占比和第三产业增加值占比。

我们考虑嵌套候选模型，总共6个候选模型。为了比较上述提到的模型选择及平均方法的性能，我们使用任何 T 个种族少数区域数据来进行模型建立，其中T = 50。然后，我们通过77-T个民族少数地区数据来预测风险。这个过程完成N = 500次，我们计算模型选择及平均方法的损失风险。结果见表6：

Table 6. Risk of model averaging and selection

表6. 模型平均及选择的损失风险

从表6可得出，构建MGPR后进行模型平均或选择得到的预测精度总体上要优于GPR，表明在多个任务同时输出的情况中，既考虑输入变量间的相关性又考虑输出变量间的相关性可以降低损失风险。

在该实例中，构建MGPR后，采用MMMA方法得到的损失风险总体程度上比S-AIC、S-BIC方法得到的损失风险要低。而构建GPR得到的结果中，总体上是MMA方法的损失风险要高于其他模型平均方法，表明输出为单个任务时，考虑输入变量之间的相关性可以降低损失风险。

3.4.2. 输出噪声异方差

采用上述提及的噪声异方差下的模型选择及模型平均方法，对本节使用的数据进行分析。本节数据来源于文献 [19] ，共38个样本。通过实验研究了5个工艺因素(反应温度、氧气分压、苯甲醇浓度(以10 ml甲苯稀释的mmol为单位)、Mn的百分比和K:Mn的比)在苄醇转化率和转频率(TOF)内稀释的影响。苄醇转化率和TOF被视为过程响应变量。需要注意的是，苄醇转化率和TOF高度相关(相关系数0.65)，即可以使用MGPR。

我们考虑嵌套候选模型，总共5个候选模型。为了比较上述提到的模型选择及平均方法的性能，我们使用27个样本数据来进行模型建立，使用剩余的11个数据进行预测。这个过程完成500次，我们计算模型选择及平均方法的损失风险。结果见表7：

Table 7. Risk of model averaging and selection

表7. 模型平均及选择的损失风险

由表7可得，大多数情况下，采用的模型平均方法，总体上是MGPR得到的预测风险要低于GPR且在单个任务构建GPR时，可以得出JMA方法优于MMA方法，与上述数值模拟结论一致。

4. 总结和讨论

在本文中，考虑到现实情况大多数为多个相关的输出变量，可以构建MGPR，并将不同的候选模型进行模型选择及平均，然后将不同任务分别构建GPR，同样将候选模型进行模型选择及模型平均，将二者得到的损失风险进行比较。数值结果表明，均值函数为线性时，噪声同方差时，构建多任务高斯过程回归模型后，再进行模型平均的损失风险要低于多个任务分别进行高斯过程回归后模型平均的损失风险；而对于模型选择方法，MGPR与GPR没有明显的优势。对于GPR建模来讲，MMA方法表现较差，表明在输出任务之间存在相关性时，采用MGPR建模后进行模型平均得到的损失风险更低，预测结果更优。

NOTES

^*通讯作者。