DOI: 10.12677/sa.2025.143071, PDF, HTML, XML,  被引量    科研立项经费支持
作者: 张 钰, 肖鸿民^*：西北师范大学数学与统计学院，甘肃 兰州
关键词: 死亡率预测；Heligman-Pollard模型；DNN；GMRF模型；Mortality Prediction； Heligman-Pollard Model； DNN； GMRF Model

文章引用：张钰, 肖鸿民. 贝叶斯推断与深度神经网络在死亡率预测中的应用与对比研究[J]. 统计学与应用, 2025, 14(3): 194-204. https://doi.org/10.12677/sa.2025.143071

1. 引言

死亡率数据分析和未来预测一直是精算师、人口学家和统计学家的研究重点。虽然静态问题较为简单，但动态问题更加复杂，至今只有部分解决方案。现有的死亡率预测模型各异，通常通过选择一个适合的模型，基于拟合优度或其他因素进行未来预测，不仅估计死亡率，还能给出预测区间以体现不确定性。图1展示了1960至2020年英国威尔士的死亡率数据。

国际上各个国家的对数死亡率的平均水平随着年龄的增长呈现出一种特定的趋势，这表明不同国家在随年龄增长的变化模式上却展现出显著的共性。进一步分析，这种共性揭示了一个重要现象：随着年龄的增长，不同国家的死亡率变化遵循着相似规律。在初始阶段，由于医疗保健、营养状况和生活环境的改善，死亡率可能呈现下降趋势。然而，随着年龄的进一步增长，人体机能逐渐衰退，慢性疾病和衰老相关疾病的发病率增加，导致死亡率重新上升。这一发现不仅深化了我们对全球人口死亡率变化规律的理解，也强调了死亡率预测在全球范围内的重要性和应用价值。因此，死亡率预测不仅是一个统计学问题，更是一项对社会福祉和全球健康具有重要意义的工作。

在1992年，Li和Lee提出了影响力深远的单因子Lee-Carter模型。这个模型简洁而经典，仅包含年龄项和时间项，被认为是最早的单因子动态死亡率模型之一，并考虑将死亡率的变化分解成随着时间而变化的时间因子和不随着时间而变化的年龄因子，进而提出了Lee-Carter模型。在文中运用美国死亡率数据对其进行拟合与预测，发现效果较好[1]。在此基础上持续研究，Li和Lee发现Lee-Carter死亡率建模和预测方法在长时间序列数据的情况下运行良好，对于数据缺失间隔不均匀的情况，采用带漂移项的随机游走模型来进行预测[2]。黄顺林和王晓军利用考虑出生年效应的Lee-Carter模型，对中国男性死亡率数据进行了精准拟合。相比其他模型，其残差分析显示拟合效果更佳。预测表明，中国男性预期寿命

Figure 1. Diagram of mortality rate data for the UK

图1. 英国的死亡率数据

将持续上升，但增幅逐步放缓[3]。李志生和刘恒甲在1992~2007年中国人口分年龄组死亡率数据基础上，对Lee-Carter模型参数拟合采用四种不同方法，即对奇异值分解法(SVD)、最小二乘法(OLS)、加权最小二乘法(WLS)和极大似然法(MLE)的拟合结果和预测能力进行了比较分析。

Dellaportas等在2001年探索了使用加权最小二乘法估计Heligman-Pollard模型的方法，并提出使用马尔科夫链蒙特卡洛(MCMC)进行贝叶斯推理[4]。随着MCMC技术的进步和应用范围的扩大，越来越多的研究开始采用贝叶斯方法进行死亡率建模。Angelos Alexopoulos等在2017年提出了一种基于Heligman-Pollard模型公式的动态模型[5]。Christian Evan Chandra等2022年结合使用截断的HP模型和Makeham模型并采用Metropolis-Hastings算法进行采样，利用贝叶斯方法估计死亡率模型的参数，这种方法允许模型在面对数据限制时整合先验信息并提供更稳健的估计[6]。

近年来，神经网络及其深度学习扩展——深度神经网络(DNN)，因其在复杂模式识别和数据拟合方面的优势，逐渐成为死亡率预测的一个重要工具。Donatien Hainaut提出了一种预测和模拟人类死亡率的神经网络方法，与有和没有群体效应的Lee-Carter模型相比，神经方法具有极好的预测能力[7]。此外，Andrea Nigri等提出了一个深度神经网络（DNN）模型，通过利用类似于人口学间接估计技术的深度学习算法，从观测或预测的预期寿命中推导出特定年龄死亡率[8]。

2. 方法与模型

2.1. 经典死亡率模型

Lee和Carter (1992)提出了著名的Lee-Carter模型，这是一种简洁而经典的单因子模型，包含年龄项和时间项。作为最早的单因子动态死亡率模型之一，Lee-Carter模型采用了对数转换，形成了一个关于年龄和时间的双线性函数，即：

$\ln m_{x,t}={\alpha }_x+{\beta }_x+{\epsilon }_{x,t}$ (1)

其中， $m_{x,t}$ 为t时刻年龄为x岁的人的中心死亡率， ${\alpha }_x$ 为不同年龄段死亡率对数变动的均值； ${\beta }_x$ 为不同年龄段死亡率对数变动的趋势； $k_t$ 为时间因素变量； ${\epsilon }_{x,t}$ 是一个误差项(白噪声)。同样的还可以通附加下述约束条件以保证参数估计结果的唯一性： ${\displaystyle {\sum }_t{k}_t}=0$ ， ${\displaystyle {\sum }_x{\beta }_x}=1$ 。

该模型根据观察到的对数死亡率可通过奇异值分解法(SVD)、最小二乘法(OLS)、加权最小二乘法(WLS)、极大似然法(MLE)拟合参数。

2.2. 深度神经网络DNN模型

深度神经网络(DNN)是一种由多层神经元组成的人工神经网络，能够自动从数据中提取特征并学习复杂模式。其结构通常包括输入层、多个隐藏层和输出层。每层神经元接收前一层输出，通过加权求和及非线性激活将信息传递至下一层。隐藏层的数量决定了网络的“深度”，使其能够捕捉数据中的深层

特征。DNN的目标是对于给定的训练数据 $\left\{\left(x_1,y_1\right),\cdots ,\left(x_n,y_n\right)\right\}$ ，找到一个函数 $f(\cdot )$ ，使得在遵循未知分布 $P\left(x,y\right)$ 的测试数据上预期误差最小化：

${\displaystyle \int L\left(y,f\left(x\right)\right)P\left(x,y\right)\text{d}x\text{d}y}$ (2)

其中， $L\left(y,f\left(x\right)\right)$ 是用于衡量给定x相对于实际值y预测误差的损失函数。本文应用的DNN模型图如图2所示，指定了在一个通用时间 $t_i$ 出生的生命期望，一个特定年龄段的死亡率向量，其结构如图所示，其中输入是出生时生命期望随时间 $t\in T$ 的向量 $e_0=\left(e_{0,t_1},e_{0,t_2},\cdots ,e_{0,t_n}\right)$

Figure 2. Diagram of DNN model

图2. DNN模型图

DNN模型采用三层架构，包括输入层、隐藏层和输出层。输入是随时间 $e_0$ 的预期寿命的向量，输出为每个年龄对应的对数死亡率。隐藏层 $H^{\left(k\right)}$ 的计算公式为：

$H^{\left(k\right)}=f^{\left(k\right)}\left(W^{\left(k\right)}{H}^{\left(k-1\right)}+b^{\left(k\right)}\right)$ (3)

其中， $f^{\left(k\right)}$ 为激活函数(使用由 $f\left(z\right)=\max \left(z,0\right)$ 定义的整流线性单元(ReLU)函数)， $W^{\left(k\right)}$ 为权重矩阵， $H^{\left(k-1\right)}$ 为前一层的输出， $b^{\left(k\right)}$ 为偏置项。

输入层接收预期寿命向量 $e_0$ ，经过隐藏层逐层传递，最终输出对数死亡率 $M=\log \left(m_{a,t}\right)$ ，其中 $a\in A$ 表示年龄， $t\in T$ 表示年份。其中设 $M=\log \left(m_{a,t}\right)$ ( $a\in A$ , $t\in T$ )为一个死亡率矩阵，行代表年龄，列代表日历年。其中 $K=3$ (输入层、隐藏层和输出层)，对于一个包含输入层、隐藏层和输出层的三层标准架构，关系如下：

$M=f^{\left(3\right)}\left(w^{\left(3\right)}{f}^{\left(2\right)}\left(W^{\left(2\right)}{f}^{\left(1\right)}\left(W^{\left(1\right)}{e}_0+b^{\left(1\right)}\right)+b^{\left(2\right)}\right)+b^{\left(3\right)}\right)$ (4)

DNN目标通过训练数据 $\left\{\left(x_1,y_1\right),\cdots ,\left(x_n,y_n\right)\right\}$ 学习一个函数使得预测输出 $\hat{M}$ 与实际输出 $M$ 之间的误差最小化。误差通过在多种数据集上表现优异的均方误差(MSE)损失函数衡量：

$\text{MSE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(M_i-{\hat{M}}_i\right)}^2}$ (5)

采用梯度下降法优化模型参数，通过反向传播算法计算损失函数关于权重和偏置的梯度，按照以下公式来更新参数：

${\omega }_{n,r}^{\left(k\right)}\leftarrow {\omega }_{n,r}^{\left(k\right)}-\eta \nabla \text{MSE}$ (6)

其中而 $\eta$ 是学习率，决定了参数更新的步幅。 $\nabla \text{MSE}$ 表示损失函数关于权重的梯度。通过这种方式，神经网络能够在训练过程中逐渐减小误差，提升模型的预测性能。

2.3. Heligman-Pollard动态模型

Heligman和Pollard在1980年提出的Heligman-Pollard模型，被广泛用于死亡率分析，包括参数化、分解和预测。该模型特别在描述不同年龄段的死亡率模式方面表现出独到之处，为精确反映特别是老年人群的死亡率变动，该模型通过引入灵活的参数设置，其复合公式有效捕捉了死亡率在不同年龄段的变化特征，公式如下：

$\frac{p_z}{1-p_z}=A^{{\left(Z+B\right)}^C}+D\exp \left[-E{\left\{\log \left(Z\right)-\log \left(F\right)\right\}}^2\right]+G{H}^z$ (7)

$p_z$ 是年龄Z的死亡率概率，公式由三部分组成，根据Alexopoulos的设置，每个参数的值的逻辑范围如下： $A\in \left(0,1\right)$ ， $B\in \left(0,1\right)$ ， $C\in \left(0,1\right)$ ， $D\in \left(0,1\right)$ ， $E\in \left(0,\infty \right)$ ， $G\in \left(0,1\right)$ ， $H\in \left(0,\infty \right)$ 。而在时间序列或空间建模中，潜在状态指的是数据生成过程中的隐藏变量，不易直接观察，而需要通过模型推断。在死亡率建模中，潜在状态可能表现在年龄或时间维度上参数的动态变化趋势，对Heligman-Pollard模型的八个参数引入动态分量，假设这些参数作为随机游走参数演化，从而放宽了死亡率曲线的平稳性假设。

为定义动态版本，设假设模型的八个参数作为随机游走参数演化，放宽死亡率曲线的平稳性假设。为定义动态版本，设，

$\phi \left(t\right)={\left({\tilde{A}}_t,{\tilde{B}}_t,{\tilde{C}}_t,{\tilde{D}}_t,{\tilde{E}}_t,{\tilde{F}}_t,{\tilde{G}}_t,{\tilde{H}}_t\right)}^{\prime }$ (8)

$\phi \left(t\right)$ 为模型参数在时间t的潜在状态，其中 $\phi \left(t\right)\in R^8$ ，元素通过变换从原始的参数变量中获得。其中t表示某一年，T表示有数据的过去年份数量。死亡的概率假设由Heligman-Pollard模型给出：

$\frac{p_z}{1-p_z}=A_t^{{\left(Z+B_t\right)}^{C_t}}+D_t\exp \left[-E_t{\left\{\log \left(Z\right)-\log \left(F_t\right)\right\}}^2\right]+G_t{H}_t^z$ (9)

其中， $Z=0,1,\cdots ,\omega$ ， $t=1,2,\cdots ,T$ ， $\omega$ 表示所采取数据集中最年长的年龄。将上述等式右边记作 $K\left(Z,{\phi }_t\right)$ ，则有

$\frac{P_{zt}}{1-P_{zt}}=K\left(Z,{\phi }_t\right)$ (10)

死亡率概率 $p_{zt}$ 及其补集 $1-p_{zt}$ ，于是有

$P_{zt}=\frac{K\left(Z,{\phi }_t\right)}{1+K\left(Z,{\phi }_t\right)}$ (11)

似然函数是针对所有时间点t和所有年龄组z的数据构建的。对于每个t和z，是观测到的个体数，而 $d_{zt}$ 是该组潜在的总人数。每个z和t的贡献是：

$\left(\begin{array}{l}{n}_{zt}\\ d_{zt}\end{array}\right){\left[K\left(Z,{\phi }_t\right)\right]}^{d_{zt}}{\left[1+K\left(Z,{\phi }_t\right)\right]}^{-n_{zt}}$ (12)

这里的二项式系数 $\left(\begin{array}{c}{n}_{zt}\\ d_{zt}\end{array}\right)$ 表示从 $d_{zt}$ 个个体中选择 $n_{zt}$ 个存活个体的不同方式的数目。 $K{\left(Z,{\phi }_t\right)}^{d_{zt}}$ 表示所有选择的 $d_{zt}$ 个存活个体的概率乘积。 ${\left[1-K\left(Z,{\phi }_t\right)\right]}^{-n_{zt}}$ 是所有可能的存活与死亡组合的总概率进行归一化。将所有年龄组和时间点的贡献乘在一起，形成整体的似然函数：

$\gamma \left(d|\phi \right)={\displaystyle \underset{t=1}{\overset{T}{\prod }}{\displaystyle \underset{z=0}{\overset{\omega }{\prod }}\left(\begin{array}{l}{n}_{zt}\\ d_{zt}\end{array}\right){\left[K\left(Z,{\phi }_t\right)\right]}^{d_{zt}}{\left[1+K\left(Z,{\phi }_t\right)\right]}^{-n_{zt}}}}$ (13)

其中，d是一个向量，其元素为 $d_{zt}$ ，且 $Z=0,1,\cdots ,\omega$ 和 $t=1,2,\cdots ,T$ 。为了对潜在状态进行动态建模，假设潜在状态 ${\phi }_t$ 遵循随机游走结构，

$\gamma \left({\phi }_t|{\phi }_{t-1},\mu ,\sum ,\eta ,\xi \right)={\varphi }_8\left({\phi }_t;{\phi }_{t-1}+\mu ,\sum ,\eta ,\xi \right)$ (14)

其中 $t=1,2,\cdots ,T$ ，这个公式表示的是一个概率分布 $\gamma \left({\phi }_{it}\right)$ ，其中 ${\phi }_{it}$ 是向量 ${\phi }_t$ 的第i个元素，当 ${\phi }_{i1}$ 属于区间 $\left[{\eta }_i,{\xi }_i\right]$ 时， $\gamma \left({\phi }_{it}\right)$ 与1成正比；当 ${\phi }_{i1}$ 不属于区间 $\left[{\eta }_i,{\xi }_i\right]$ 时， $\gamma \left({\phi }_{it}\right)$ 为0。随机游走过程通过(14)定义，对Heligman-Pollard模型的参数引入了大量的先验结构，同时放宽了这些参数在时间上的平稳性假设。本文固定了向量 $\eta$ 和 $\xi$ 的元素，向量 $\eta$ 具体为：

$\eta ={\left(-10.61,-10.61,-5.99,-11.29,-25.33,-\infty ,-17.5,-1.39\right)}^{\prime }$ (15)

向量 $\xi$ 具体为：

$\xi ={\left(-2.75,-0.2,2.2,-3.48,4.09,2.64,-3.48,0.18\right)}^{\prime }$ (16)

对于随机游走过程的漂移 $\mu$ ，假设 $\gamma \left(\mu \right)={\varphi }_8\left(\mu ,0,M^{-1}\right)$ ，其中M是一个8 × 8的对角矩阵，其对角线元素为0.001，对于协方差矩阵 $\sum$ ，采取逆Wishart分布： $\sum |\alpha ~IW\left(\nu ,\psi \right)$ ， $\nu$ 是自由度， $\psi$ 是一个 $p\times p$ 的正定矩阵，其密度函数为：

$p(\sum )=\frac{{\left|\psi \right|}^{\left(\nu -p-1\right)/2}\exp \left(-\frac{1}{2}tr\left(\psi {\sum }^{-1}\right)\right)}{2^{\left(\nu p\right)/2}{\Gamma }_p\left(\frac{\nu }{2}\right){|\sum |}^{\left(\nu +p+1\right)/2}}$ (17)

假设 $\Psi$ 是一个对角矩阵，其对角元素为 $1/{\alpha }_1,\cdots ,1/{\alpha }_8$ ，其中 $\alpha =\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_8\right)$ 。对于每一个参数 ${\alpha }_i$ 都独立同分布(IID)地从逆伽马分布抽取，

${\alpha }_1~IG\left(\frac{1}{2},\frac{1}{l^2}\right)$ (18)

设置 $l={10}^5$ ，这种逆伽马先验分布结构意味着 $\sum$ 对角线上地标准差 ${\sigma }_i$ 服从半 $t\left(\nu ,l\right)$ 分布，选 $\nu =2$ ，用 $\theta =\left(\Sigma ,\mu ,\alpha \right)$ 表示模型参数，模型后验分布为：

$\gamma \left(\phi ,\theta |d,\eta ,\xi \right)\propto \gamma \left(\theta \right)\gamma \left(d|\phi \right){\displaystyle \underset{t=2}{\overset{T}{\prod }}{\varphi }_8\left({\phi }_t;{\phi }_{t-1}+\mu ,\sum ,\eta ,\xi \right)}$ (19)

提出的模型是一个基于潜在高斯模型与 $\phi$ 和超参数 $\theta$ 的分布的框架。为了从后验分布中获得样本，利用Metropolis-within-Gibbs采样器，交替从 $\gamma \left(\theta |d\right)$ 和 $\gamma \left(\theta \right)$ 中采样，条件分布所需样本从 $\gamma \left(\theta |d\right)$ 采样。从 $\gamma \left(\theta |d\right)$ 采样通过7个Metropolis-Hastings步骤更新每个 $\phi$ 来执行。

2.4. GMRF动态模型

马尔可夫随机场描述一组变量 $\left(x_1,\cdots ,x_n\right)$ 的联合分布，该分布是通过其全条件分布 $\pi \left(x_i|x_{-i}\right)$ 来确定的，其中 $x_{-i}$ 表示除 $x_i$ 之外的所有变量集合。当这些条件分布是高斯分布时，这种马尔可夫随机场被称为高斯马尔可夫随机场(GMRF)，可通过内在高斯马尔可夫随机场模型模拟死亡率。基于似然函数公式(13)相关模型模拟死亡率，将每年各年龄组的死亡概率 $p_{zt}$ 转换成变量 $x_{zt}=\log \left[\frac{p_{zt}}{1-p_{zt}}\right]$ ，对所有 $z=0,\cdots ,\omega$ 和 $t=1,\cdots ,T$ 。其中，

$x_t={\left(x_{0t},\cdots ,x_{\omega t}\right)}^{\prime }$ (20)

并设

$x=\left({x^{\prime }}_1,\cdots ,{x^{\prime }}_T\right)$ (21)

这为一个 $\left(\omega +1\right)T$ 维向量，包含了所有时间点和所有年龄组的数据，想象成一个有 $\left(\omega +1\right)\times T$ 个节点的网格，其中z表示行(年龄组)，t表示列(年份)。对于向量x，假设遵循一个 $\left(\omega +1\right)T$ 维高斯分布时，均值为 $\mu \left(b_1{1}_{\omega +1},2b_1{1}_{\omega +1},\cdots ,T{b}_1{1}_{\omega +1}\right)$ ，其中 $1_{\omega +1}$ 是全1的 $\left(\omega +1\right)$ 维向量，精度矩阵为：

$Q=\tau \left({\rho }_{age,R_{\omega +1}}\otimes I_T+{\rho }_{year,1_{\omega +1}}\otimes R_T\right)$ (22)

这里， $R_{\omega +1}$ 是一个 $\left(\omega +1\right)\times \left(\omega +1\right)$ 矩阵，其中 $R_{ij}$ 满足：若 $i=1$ 和 $j=1$ 时 $R_{ij}=1$ ；若 $i=\omega +1$ 和 $j=\omega +1$ 时 $R_{ij}=1$ ； $i=j$ 且 $i,j\ne 1$ 且 $i,j\ne \omega +1$ 时 $R_{ij}=2$ ；若 $\left|i-j\right|=1$ 时， $R_{ij}=-1$ ；其余情况 $R_{ij}=0$ 。x本质是一个高斯马尔可夫随机场，矩阵Q是奇异的。对于 $Z=0,1,\cdots ,\omega -1$ 和 $t=1,2,\cdots ,T-1$ ，随机变量 $x_{zt}$ 的完全条件分布是正态分布，其均值为：

$\frac{1}{4}\left\{{\rho }_{age}\left(x_{z-1,t}+x_{z+1,t}\right)+{\rho }_{year}\left(x_{z,t-1}+x_{z,t+1}\right)\right\}$

方差为 $\frac{1}{4\tau }$ 。参数 ${\rho }_{age}$ 和 ${\rho }_{year}$ 分别控制跨年龄组和跨年份的死亡概率的关联性，捕捉了年龄纬度和日历时间纬度的相关性，被期望为不同值，为保证能够唯一确定模型的参数，假设 ${\rho }_{age}+{\rho }_{year}=2$ 。

2.5. 贝叶斯推断

此模型的似然函数由等式(13)右侧各项的乘积给出：

$\gamma \left(d|x\right)={\displaystyle \underset{t=1}{\overset{T}{\prod }}{\displaystyle \underset{z=0}{\overset{\omega }{\prod }}\left(\begin{array}{l}{n}_{zt}\\ d_{zt}\end{array}\right)p_{zt}^{d_{zt}}{\left[1-p_{zt}\right]}^{n_{zt}-d_{zt}}}}$ (23)

其中d是一个 $\left(\omega +1\right)T$ 维向量，元素为 $d_{zt}$ 。模型参数表示为 $\theta =\left(b,{\rho }_{age},\tau \right)$ ，构建MCMC算法用于从联合后验分布中采样，其密度为：

$\gamma \left(\theta ,x|d\right)\propto \gamma \left(\theta \right)\gamma \left(x|\theta \right)\gamma \left(d|x\right)$ (24)

其中， $\gamma \left(\theta \right)$ 是参数的先验分布的密度， $\gamma \left(x|\theta \right)$ 是均值为 $\mu$ 精度矩阵为Q的 $\left(\omega +1\right)T$ 变量高斯分布的密度。模型参数b、 ${\rho }_{age}$ 和 $\tau$ 的完全条件分布中采样。对于 ${\rho }_{age}$ ，采用随机游走的Metropolis-Hastings步骤。 $\gamma \left(x|\theta ,d\right)$ 中采样，使用基于梯度的辅助MCMC采样器首先构建辅助变量 $u_k\in {ℝ}^{\left(\omega +1\right)T}$ 是从高斯分布中抽取的：

$N\left[x+\left(\frac{\delta }{2}\right)\nabla \log \left\{\gamma \left(d|x,\theta \right)\right\},\frac{\delta }{2}{I}_{\left(\omega +1\right)T}\right]$ (25)

计算状态x和 $\theta$ 下对数似然的梯度 $\nabla \log \gamma \left(d|\theta ,x\right)$ ，有

$q\left(x_{prop}|x\right)\propto {\left(\frac{\delta }{2}\right)}^{\left(\omega +1\right)T}{\text{e}}^{-\frac{1}{\delta }{\Vert x_{prop}-u_x\Vert }^2}\gamma \left(x_{prop}|\theta \right)$ (26)

从(26)提取新的值 $x_{prop}$ 中，并且选择 $x_{prop}$ 以Metropolis-Hastings接受概率 $\min \left(1,\alpha \right)$ 被接受，其中

$\alpha =\frac{\gamma \left(d|x_{prop},\theta \right)}{\gamma \left(d|x,\theta \right)}\exp \left\{f\left(u_x,x_{prop}\right)-f\left(u_x,x\right)\right\}$ (27)

且

$f\left(u_x,x\right)={\left(u_x-x-\left(\frac{\delta }{4}\right)\vee \nabla \log \gamma \left(d|\theta ,x\right)\right)}^{\text{T}}\vee \nabla \log \gamma \left(d|\theta ,x\right)$ (28)

调整参数 $\delta$ 以实现50%到60%的接受率。对于每个z和k，目标为预测在未来时间点 $t=T+1,T+2,\cdots ,T+k$ 的死亡概率 $p_{zt}$ ，通过向量 $x^{\star }=\left(x_{T+1},\cdots ,x_{T+k}\right)$ 表达。所需的预测密度为：

$\gamma \left(x^{\ast }|d\right)={\displaystyle \int \gamma \left(x^{\ast }|\theta ,x\right)\gamma \left(\theta |x,d\right)\gamma \left(x|d\right)\text{d}x\text{d}\theta }$ (29)

3. 实证分析

3.1. 数据分析及预处理

在本研究中，从人类死亡数据库中选取了智利和日本的死亡率数据作为研究对象。这些国家代表了不同的经济和社会发展水平，能够提供具有广泛适用性的分析结果。从公开的死亡率数据库中提取了1995年至2020年的死亡率数据。该数据涵盖了0到89岁年龄段的死亡率，并按年份和年龄段进行分类，以确保每个年龄组的死亡率数据充分且具有代表性。

为了后续模型训练和评估提供可靠的基础以确保数据的质量和一致性，首先，数据从人类死亡率数据库中读取数据。接着，对“年龄”进行处理，使用正则表达式移除可能存在的“+”符号，并将该列转换为整数类型，以确保年龄数据的格式统一。随后，对数据进行筛选，仅保留“年份”列在1995年至2020年之间且“年龄”列小于等于89的数据。这一步骤的目的是排除异常值或不符合要求的数据，确保数据的合理性和有效性。接下来，处理“死亡率”列中的异常值。由于后续步骤需要对死亡率取对数，而死亡率为0时会导致对数计算无意义，因此将死亡率为0的值替换为一个非常小的正数。然后，对“死亡率”进行对数变换，生成新的“对数死亡率”列。这一步骤的目的是将偏态分布的数据转换为接近正态分布的形式，使其更适合线性模型的假设。为了消除不同特征之间的量纲差异，数据被进一步缩放。通过最小-最大缩放方法，将“年份”和“年龄”列的值缩放到[0, 1]的范围内。最后，数据被分为训练集和测试集。这一步骤的目的是确保模型能够在未见过的数据上进行评估，从而验证其泛化能力。

图3、图4分别是Heligman-Pollard模型和GMRF在预测2018年日本死亡率中的表现。

Figure 3. Diagram of HP model result

图3. HP模型结果图

Figure 4. Diagram of GMRF model result

图4. GMRF模型结果图

为了比较四种模型的效果，最终选择使用智利的预测数据进行对比分析。在经典的LC模型中，我们分别对全年龄组(0~89岁)和老年组(40~89岁)进行了预测。而在DNN模型中，数据按照年份划分，2001年之前的数据用于训练集，2001年及之后的数据则作为测试集。

3.2. 模型评估

本文采用的评价指标有均方根误差(RMSE)、平均绝对误差(MAE)和平均绝对百分比误差(MAPE)，表达式为：

$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(y_i-{\hat{y}}_i\right)}}$ (26)

$\text{MAE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left|y_i-{\hat{y}}_i\right|}$ (27)

$\text{MAPE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left|\frac{y_i-{\hat{y}}_i}{y_i}\right|}\times 100\%$ (28)

其中 $y_i$ 为真实值， ${\hat{y}}_i$ 为预测值，n为样本总数。

3.3. 模型对比分析

利用测试集对组合模型进行测试，将建立之后不断调优的动态模型和神经网络模型得出最终的预测结果，对比结果如表1所示。

Table 1. Model parameter settings

表1. 模型对比结果

	RMSE	MAE	MAPE (%)
LC	0.1687	0.1306	13.0862
LC(局部)	0.0054	0.0023	10.6682
DNN	0.0150	0.0088	17.224
HP	0.0014	0.0006	8.4098
GMRF	0.0006	0.0002	3.8572

4. 结论

Lee-Carter模型单因子假设限制了其对复杂死亡率模式的捕捉能力。实证结果显示其预测误差较高(RMSE: 0.1687, MAE: 0.1306, MAPE: 13.0862)，局部数据表现有所改善(RMSE: 0.0054, MAE: 0.0023, MAPE: 10.6682)但仍未达到最优。DNN模型通过多层神经元结构捕捉死亡率数据中的非线性关系，采用梯度下降法和反向传播算法优化参数。尽管其理论拟合能力强，但受数据量和超参数调优的影响，实证预测误差(RMSE: 0.0150, MAE: 0.0088, MAPE: 17.224)高于Heligman-Pollard和GMRF模型。Heligman-Pollard模型通过引入动态分量和随机游走结构，放宽了参数平稳性假设以更好地捕捉死亡率的动态变化。采用贝叶斯推断和Metropolis-within-Gibbs采样器，其预测精度较高(RMSE: 0.0014, MAE: 0.0006, MAPE: 8.4098)，在描述不同年龄段死亡率模式方面表现优异。GMRF模型基于梯度的辅助MCMC采样器提高了采样效率，预测误差最低(RMSE: 0.0006, MAE: 0.0002, MAPE: 3.8572)，能够有效捕捉时空依赖性，为死亡率预测提供了可靠框架。

总得来说，GMRF模型预测精度最优，适用于高精度预测和复杂时空依赖性场景；Heligman-Pollard模型适合描述死亡率动态变化；DNN模型适用于数据充足且需捕捉非线性关系的场景，但需注意超参数调优；Lee-Carter模型简洁且解释性强，但在复杂数据上表现有限。未来改进方向可结合GMRF的动态建模能力和DNN的非线性拟合能力开发混合模型，或引入社会经济因素等外部变量，可能进一步提升预测精度。

基金项目

国家自然科学基金项目(NO.12061066)；甘肃省自然科学基金项目(NO. 20JR5RA528)。

NOTES

^*通讯作者。

参考文献

[1]	Lee, R.D. and Carter, L.R. (1992) Modeling and Forecasting U. S. Mortality. Journal of the American Statistical Association, 87, 659-671. https://doi.org/10.2307/2290201
[2]	Li, N., Lee, R. and Tuljapurkar, S. (2004) Using the Lee-Carter Method to Forecast Mortality for Populations with Limited Data. International Statistical Review, 72, 19-36. https://doi.org/10.1111/j.1751-5823.2004.tb00221.x
[3]	黄顺林, 王晓军. 加入出生年效应的死亡率预测及其在年金系数估计中的应用[J]. 统计与信息论坛, 2010(5): 81-86.
[4]	Dellaportas, P., Smith, A.F.M. and Stavropoulos, P. (2001) Bayesian Analysis of Mortality Data. Journal of the Royal Statistical Society Series A: Statistics in Society, 164, 275-291. https://doi.org/10.1111/1467-985x.00202
[5]	Alexopoulos, A., Dellaportas, P. and Forster, J.J. (2018) Bayesian Forecasting of Mortality Rates by Using Latent Gaussian Models. Journal of the Royal Statistical Society Series A: Statistics in Society, 182, 689-711. https://doi.org/10.1111/rssa.12422
[6]	Chandra, C.E., Abdullah, S. and Devila, S. (2022) Estimating Indonesian Complete Life Table and Fair Annual Pure Premium Range from Abridged Life Table with Bayesian Method and Bootstrapping. International Journal on Advanced Science, Engineering and Information Technology, 12, 2226-2236. https://doi.org/10.18517/ijaseit.12.6.14947
[7]	Hainaut, D. (2018) A Neural-Network Analyzer for Mortality forecast. ASTIN Bulletin, 48, 481-508. https://doi.org/10.1017/asb.2017.45
[8]	Nigri, A., Levantesi, S. and Aburto, J.M. (2022) Leveraging Deep Neural Networks to Estimate Age-Specific Mortality from Life Expectancy at Birth. Demographic Research, 47, 199-232. https://doi.org/10.4054/demres.2022.47.8