通过贝叶斯非参数模型探索在物流风险评估中的大数据

Abstract

在货物物流中，一个关键的绩效指标是运输风险，定义为实际到达时间与计划到达时间的偏差。对于客户和货运代理来说，无论是早到还是晚到都是不可取的。本文使用了服务于20家航空公司、拥有1336条航线的一家领先的货运代理公司提供的半年空运数据，并对运输风险进行了评估和预测。有趣的是，我们的初步数据分析显示，运输风险具有很强的多模态（multimodal） 以两种及以上运输方式完成的货物运输;另外，intermodal，货物全程由一种且不变的单元或道路车辆装卸。 特征，这是由未观察到的事件驱动的，比如货物失踪航班。为了适应这一特点，我们引入了贝叶斯非参数模型—probit stick-breaking process mixture model—来灵活估计运输风险的条件密度函数。我们证明，使用其他方法可能导致误导的推论。
我们的模型为货运代理提供了定制价格和服务报价的功能。 **（根据我的理解，这边应该是为货运代理在加上了各种风险预估以后，可以在保证货运代理不折本的情况下，给客户一个合理的报价）它还可以生成基准航空公司（baseline airline）业绩，以提供公平的供应商评估。 （我的理解是货运代理在运送货物的时候需要选择不同的航空公司，所以我们的模型对航空公司的业绩也进行了一个评估，货运代理可以选择合适自己需要的航空公司）此外，该方法允许我们将周期性风险（recurrent risks）与破坏性风险（disruption risks）分离开来。这一点很重要，因为针对这两种风险的对冲（减低另一项投资风险的投资）**策略往往截然不同。

Introduction

近几十年来，全球贸易大幅增长。许多公司现在拥有海外设施和供应链合作伙伴。因此，国际货物物流管理在全球经济中发挥着越来越重要的作用。航空运输以具有竞争力的价格向世界各地的客户配送对时间敏感、价格昂贵、易腐烂或用于即时供应网络的货物。事实上，航空货运每年运送价值超过6.4万亿美元的货物。这大约占世界贸易总值的35%。Boeing预测，包括快递（express traffic）在内的这一行业在未来20年的年均增长率将达到4.7%，达到2013年创下的收益延吨公里(revenue tonne-kilometers，RTKs)的两倍多。然而，对这个行业的关注却少得出奇:航空货运一直是相对于更引人注目的乘客业务而言的穷亲戚。

忽视这一点的后果是，货物运输的服务水平已成为各公司的一大担忧。在货物物流中，一个关键的性能指标就是运输风险（运输的可靠性），定义为实际到达时间和计划到达时间的误差：
t r a n s p o r t r i s k = a c t u a l a r r i v a l t i m e − p l a n n e d a r r i v a l t i m e transport\ risk = actual\ arrival\ time−planned\ arrival\ time transport risk=actual arrival time−planned arrival time

既不早到也不晚到对于客户和货运代理是最理想的情况。延迟交货会导致生产和向所有下级客户交付产品/服务的延迟，提前交货会导致额外的库存和处理成本。极端风险，即超过48小时的延误或超过24小时的提前，被认为是(运输)破坏风险（disruption risks），因为它严重影响客户和货运代理的运营。为了区分一天内的破坏性风险（disruption risks）和常规偏差，我们将后者称为周期性风险（recurrent risks）。2011年PRTM的一项调查显示，69%的公司将提高交付配送水平作为其供应链管理的首要策略。在Infosys 2010年的一份报告中，物流行业将“运营商延迟交货和交货不佳”排名为主要的风险。此外，在国际航空运输协会(International Air Transport Association，IATA) 2014年对主要货运代理及其客户进行的一项调查中，低可靠性被认为是第二个最重要的因素(仅次于运输成本)。

本文以2012年至2013年半年的国际航空货代数据为基础，对由一个领先的货运代理提供的20家航空公司提供的1336条航线的数据进行了研究。使用贝叶斯非参数（Bayesian nonparametric, BNP）模型-the probit stick-breaking(PSBP) mixture model-我们就获得了运输风险分布和破坏性风险（disruption risks）可能性的准确评估。

Empirical Air Cargo Transport Risk Distribution

我们的工作似乎是供应链文献中对全球航空物流的第一个实证研究。从数据中可以观察到一个有趣的现象，即取决于预测因子(即，自变量包括航空公司、航线、发货时间、货物重量等)的运输风险分布是多模态（multimodal）的（意思就是由两种或以上的交通工具完成的货物运输），如图1所示。图1的左侧是数据中观察到的所有发货(近9万发货)的运输风险的经验分布，很明显，这是一个不对称的长尾分布，在正值处有几个凸起。为了更好地可视化凸起，我们只在图1的右侧绘制范围(0,150)内的数据。在这里，我们可以清楚地看到，大的凸起集中在的天数(24小时，48小时，72小时等等)和小凸起集中在这些天之间。**这些系统性的高峰（systematic peaks）在很大程度上是因为飞机没能按照原定计划飞行，后来又在同一航线上装载了一架飞机。**国际航班的飞机间隔通常在24小时左右，国内航班之间的航班间隔通常在几个小时左右，这在很大程度上取决于航线。因此，定期航班之间的时间间隔就会转化为在运输风险的条件分布中的不同高峰之间的间隔，从而形成多模态（multimodal）分布。

以往的实证研究主要集中在国内旅客的航班到达或离开的延误上。大多数文献资料认为延误服从单模态分布 （unimodal,使用单一的运输方式，与多模态的意思相反） ，并且是线性模型。然而，延迟分布呈现出明显的多模态（multimodal），使得线性模型不适用于航空货物运输风险评估和预测。先前对线性模型的关注可能是由于使用了美国交通部(DOT，Department of Transportation)在每次飞行时收集的数据。相反，我们的数据是在每次货运时收集的（at the level of each cargo trip），包括从开始到结束的运输信息(information on trip from the beginning to end)，通常包括几个中转航班。全程数据让我们可以探索以前没有考虑过的新的交通不确定性(explore new transport uncertainties not considered before)。具体地说，我们包括运输风险中因航班延误的信息。 由于全行程时延分布(delay distribution)具有明显的多模态性，本文提出了一种新的建模方法。 然而，我们所开发的建模方法并不局限于货物运输风险，还可以应用于其他运输风险(如旅客航空运输风险)。对于航空旅客来说，运输风险是由旅客到达目的地的时间决定的。旅客到达时间可能与计划的最后一次航班的到达时间不同，因为旅客可能会因为他或她之前在中转机场的航班延误而错过最后一次航班。从这个角度看，客运与我们研究的货物运输风险问题是相似的。

据我们所知，最接近我们工作的人是Tu et al。作者使用了一个包括季节性、日变化趋势和残差三部分的模型来研究航班的起降（departure and push-back）延误，建立了一个使用高斯核和固定混合权重的混合分布模型。研究人员用该模型拟合了美国联合航空公司(United Airlines)运营的丹佛国际机场(Denver International Airport)航班一年的数据。然而，他们模型的简单性使其不适用于我们的研究，因为我们的研究具有更不稳定的运输风险(如上所述)，而且数据集更大，包含运往95个国家20多家航空公司运营的200多个机场的货物。为了适应这种复杂性并更好地匹配我们的数据，我们使用了BNP模型。在与替代模型的广泛的模型比较中（in extensive model comparisons with alternative models），包括一个灵活的混合模型，我们的模型显示了优越的性能。

BNP模型和条件分布函数

我们的第二个贡献是方法论。为了适应经验运输风险分布中的多模态特征，我们引入了一种最先进的贝叶斯统计工具——BNP混合模型。据我们所知，之前没有任何工作在实证经营管理（Empirical Operation Management）中使用过相关技术，到目前为止，实证经营管理（Empirical Operation Management）主要应用频率统计，如普通最小二乘估计或最大似然估计。

贝叶斯统计在过去的二十年中经历了快速的发展，计算能力的不断提高使其得到了加速。在这些工具中，BNP混合模型在最近几年变得流行起来，在金融、计量经济学、遗传学和医学等领域都有应用。

在我们的应用程序中，我们采用了一个特定的BNP模型—PSBP混合模型。 该方法具有灵活性强、通用性强、计算方便等优点。此外，PSBP可以在弱正则条件下对任意条件下的条件密度进行一致性估计。Rodriguez等人(2009)利用该技术建立了一个非参数因子模型，研究预测DNA损伤和修复的遗传因子的预测。Chung和Dunson(2009)应用这项技术开发了一个非参数变量选择框架。我们的模型旨在在所有范围内捕捉运输风险的分布特征，包括周期性风险（recurrent risks）和破坏性风险（disruption risks）。

特别地， 我们关注于在PSBP框架中对运输风险的条件分布进行建模。 通过对条件分布的建模，我们可以研究运输风险与潜在预测因子（包括航空公司、航线、运输时间、货物重量等）之间的关系，并在此基础上进一步探索提高运输可靠性的方法。

为了证明PSBP的价值，我们将我们的运输风险估计与一个朴素线性模型得出的结果进行了比较。我们证明这两种方法的结果有显著的不同。例如，朴素线性模型没有捕捉到航空公司在 运输服务级别（transport service levels） 中起着关键的作用，更重要的是，低估了破坏性风险（disruption risks），这可能导致风险管理策略不足。把我们的模型和两个广义和先进的替代模型（generalized additive models (GAM) and flexible mixture models ）进行了比较。 我们的 PSBP 模型具有较强的样本内和样本外预测能力，但计算时间相对较长。

数据驱动的风险评估工具

我们的方法为供应链风险评估提供了一种强大而通用的工具，但它并没有得到应有的重视。 尤其是，尽管供应链风险管理越来越受到从业者和学术界的关注，但是最近麦肯锡公司（McKinsey&Co.）对企业高管的全球调查显示，“近四分之一的公司表示，他们的公司没有正式的风险评估。”相反，正如Van Mieghem(2011)所述，通过运行进行风险管理包括四个步骤，1)危害识别;(2)风险评估;(3)战术风险决策;(4)实施战略风险缓解或对冲。这四个步骤必须反复执行和更新。在这四个步骤中，第一步是更多的基于经验和上下文的，这通常需要从轶事记录或与特定业务流程相关的长期经验中获得信息。第四步是以行动为基础的，需要详细的组织设计和信息系统来执行步骤3中制定的对冲策略。这两个步骤可能不需要定量方法。相反，步骤2和步骤3需要严格的分析和量化，因此需要进行分析研究。虽然大多数供应链风险管理文献集中在第三步，这涉及到制定策略来降低负面事件发生的概率和（或）它们发生的后果， 本文集中在第2步风险评估。

风险评估包括对两个组成部分的估计:(a) 风险可能性，即:“发生不良事件或危险的可能性”及 (b) 风险影响，即“不良事件的后果”。 长期的风险评估是这两部分的整合。Kleindorfer等人(2003)利用环境保护署收集的数据评估了灾难性化学事故的风险影响。Kleindorfer和Saad(2005) 提出了一个供应链面临中断时风险评估和风险缓解的概念框架。与这些研究不同的是，我们的工作侧重于 使用统计方法 来准确估计风险可能性，这需要更先进的科学计算和分析工具。正确识别危害和评估风险对替代管理政策的有效性具有重要意义。我们的研究显示，我们的研究表明，谨慎的风险评估对于为客户（即托运人）开发定制服务（develop tailored services）和选择服务供应商至关重要。此外，我们的方法能够将周期性风险(recurrent risks)与破坏性风险（disruption risks）分离开来。这一点很重要，因为针对这两种风险的对冲策略（处理办法）往往大相径庭。

本文研究的运输风险与许多作者研究的制造业随机产量/产能风险相似。 此外，运输中断风险是文献中考虑的随机供应中断风险的一个重要类型或组成部分 虽然这些作者关注的是假设特定风险分布的风险缓解策略，比如破坏风险（disruption risks）的伯努利分布， 但是本文介绍的贝叶斯PSBP混合模型可用于在数据可用时生成经验随机收益率分布（empirical random yield distributions）和破坏概率(disruption probabilities)。

论文的其余部分组织如下:第二部分，我们简要介绍了航空货运物流业的概况及其面临的挑战，我们在本研究中使用的数据以及我们提出的研究问题。在第三部分中，我们描述了模型选择的方法，介绍了PSBP混合模型和后验计算算法，并基于可靠度和预测性能与其他模型进行了比较。在第四部分中，我们解释了结果。第五部分，我们建议应用我们的模型来设计更有效的操作策略。第六部分对全文进行了总结，并对未来的研究方向进行了展望。附录A包含数据清理步骤，以及说明数据的表和图。附录B包含了一定的算法细节、模型实现步骤、模型检验和比较结果。附录C包含估计结果和选择的图。

模型

在本节中，我们将详细解释该模型。3.1节为使用PSBP混合模型估算运输风险和优势的条件分布提供了动力（motivation）。这个模型可以分解为两部分：混合权重（mixture weight）和混合核（mixture kernel），具体分别在3.2节和3.3节。在3.4节展示了贝叶斯后验采样算法来估计权值和核数的未知参数估计方法。在3.5节中讨论了模型的选择和比较。详细的补充材料见附录B部分。

条件风险分布

多模态特征不仅存在于总体的数据层次（at the aggregate data level）(见图1)，还存在于粒度层次(at the granular level)(如每条航线)。图3中的直方图显示了两家航空公司提供的两条样本航线上的经验分布。为了根据这些数据做出准确的预测和推断，第一步是选择一个足够灵活的模型来很好地拟合数据。多模态数据的通常的模型选择依赖于混合（mixtures）,例如正常内核的混合（mixtures of normal kernels）。

我们在研究基于需求和决策变量的运输风险分布（the distribution of transport risks conditional on demand and continuous predictors），包括分类预测因子和连续预测因子时，不能只依赖于简单的混合模型。这会导致条件分布估计的问题。一系列关于弹性条件分布估计（flexible conditional distribution）的文献都采用了频率方法（frequentist methods）。Fan等人(1996)提出了一种“双核局域线性研究方法（a double-kernel local linear approach）”，Hall等人(1999)和Hyndman和Yao 等人(2002)还提出了一种与频率相关的方法（related frequentist methods）。其他的一些主流的选择是BNP混合模型。Muller等人（1996）提出了一种非线性回归的贝叶斯方法，其中作者使用一种Dirichlet过程正态混合模型（Dirichlet process mixture of normal）来模拟因变量和自变量的联合分布。这种方法通过响应和预测器的联合模型得出响应的条件分布模型。虽然这种联合模型证明是灵活的，但在实践中，与直接针对条件响应分布的模型相比，它们有明显的缺点，而不需要对与预测器联合密度对应的高维损害参数进行建模（model the high-dimensional nuisance parameter corresponding to the joint density of the predictors）。这些缺点包括将自变量视为随机变量，而它们通常是设计变量(例如，将航线或航空公司视为随机变量似乎是不自然的)，以及在估计条件分布方面具有相对较差的实际性能。

相反，我们只关注给定预测因子 X = ( x 1 , … , x p ) ′ ∈ X X = (x_1,\dots,x_p)^{'}\in \mathcal{X} X=(x1,…,xp)′∈X（ X \mathcal{X} X 是预测因子 X X X的样本空间）的未知条件下运输风险 y y y分布的直接建模，而不指定 X X X的边缘模型。在上下文中，预测因子 X X X={航空公司（ a a a, airline）、航线（ r r r, route）、月份（ m m m, month）、中转次数（ l e g leg leg, number of leg）、初始偏差（ d e v s t a r t dev_{start} devstart, initial deviation）、计划持续时间（ d u r dur dur, planned duration）、货物重量（ w g t wgt wgt, cargo weight）、货物件数（ p c s pcs pcs, cargo number of pieces）（如表1所示）。特别是，我们认为运输风险 y y y是由卷积得到的(arises from a convolution):

y ∣ X ∼ ∫ k ( y ∣ ψ ) G ( d ψ ) ) , ( 1 ) y|X \sim \int k(y| \mathbf{\psi})\mathbf{G}(d\mathbf{\psi})),(1) y∣X∼∫k(y∣ψ)G(dψ)),(1)

其中， k ( ⋅ ∣ ψ ) k(\cdot|\psi) k(⋅∣ψ) 是由 ψ \psi ψ限定的一个参数核（例如，正常核 k ( ⋅ ∣ ψ k(\cdot|\psi k(⋅∣ψ）是由 ψ = ( m e a n , s t a n d a r d d e v i a t i o n ) \psi=(mean, standard deviation) ψ=(mean,standarddeviation)),并且允许混合分布 G X G_X GX 随预测器 X ∈ X X\in\mathcal{X} X∈X 灵活变化。BNP文献中的典型形式是：
G X = ∑ l = 1 L ω l ( X ) δ ψ l ( X ) ′ G_X = \sum_{l=1}^{L}\omega_l(X)\delta_{\psi_l(X)'} GX=l=1∑Lωl(X)δψl(X)′
其中， ∑ l = 1 L ω l ( X ) = 1 \sum_{l=1}^L\omega_l(X) = 1 ∑l=1Lωl(X)=1并且 ω l ( X ) ≥ 0 , ( 2 ) \omega_l(X)\geq0,(2) ωl(X)≥0,(2)

这里，原子 { ψ l ( X ) : X ∈ X } l = 1 L \{\psi_l(X):X\in \mathcal{X}\}_{l=1}^L {ψl(X):X∈X}l=1L 是 X \mathcal{X} X 上随机过程的独立同分布（identitical distribution）样本路径， { ψ l ( X ) , X ∈ X } \{\psi_l(X),X\in \mathcal{X}\} {ψl(X),X∈X} 是对所有 X X X来说的一个和为1的相关预测因子的概率权重。（predictor-dependent probability weights that sum to one for all x）上述形式太笼统，而并不能实用性，所以有必要对实际应用作一些简化。

一个常见的可能性方法是仅在 G X G_X GX 原子和 ψ l ( X ) \psi_l(X) ψl(X) 中引入预测因子依赖性，同时保持权重 ω l ( X ) = ω l \omega_l(X) = \omega_l ωl(X)=ωl 固定。然而，这种方法在我们的实验中往往表现得相对较差，包括空运货物运输风险数据，正如与Flexmix模型的比较中可以看出。

在我们的例子中，因变量运输风险的峰值位置几乎是恒定的(即在美国，国际货运每天都有高峰，而国内货运除了每天的高峰外，每小时还有几个高峰)。然而，高峰的高度会随着 X X X 的变化而发生很大的变化(例如，航线、航空公司、需求变量)。每一个峰值的高度(大致)代表了观测结果落入以峰值为中心的内核的可能性。例如，如果以某一特定的 X 1 X_1 X1 为条件，24小时左右的峰值相对较高，那么以 X 1 X_1 X1 为条件的发货有很大的可能被延迟24小时。相反，如果以某个 X 2 X_2 X2 为条件，在0小时左右只有一个显著的峰值，那么以 X 2 X_2 X2 为条件的发货可能接近计划的到达时间。因此，在我们的上下文中，找出每一峰值的高度如何依赖于 X X X 是一个重要的问题。

在权重中诱导依赖结构是困难的（Inducing dependence structure in the weights can be difficult），导致计算算法复杂和效率低下，限制了模型的应用。 为了克服这些困难，我们采用了PSBP混合模型，它具有在弱正则条件下(weak regularity conditions)计算能力强、一致性好等优点。

Bayesian PSBP

回顾一下公式(2)中混合测度的一般形式，破棒权重（stick-breaking weights）被定义为 ω l = u l ∑ p < l ( 1 − u p ) \omega_l = u_l\sum_{p<l}(1-u_p) ωl=ul∑p<l(1−up) ，其中，破棒率(stick-breaking ratios)是独立分布的 u l ∼ H l ( l < L ) u_l \sim H_l(l<L) ul∼Hl(l<L) 并且 u l = 1 ( L 是有限个的实例（ c a s e ） ) u_l=1(L是有限个的实例（case）) ul=1(L是有限个的实例（case）)。在没有预测因子的基线情况下（in the baseline case ），概率破棒 ( p r o b i t s t i c k − b r e a k i n g ) (probit \ stick-breaking) (probit stick−breaking) 权重构造为
u l = Φ ( γ l ) , γ l ∼ N ( μ , ϕ ) , u_l = \Phi(\gamma_l), \gamma_l \sim \textbf{N}(\mu, \phi), ul=Φ(γl),γl∼N(μ,ϕ),
其中， Φ ( ⋅ ) \Phi(\cdot) Φ(⋅) 表示累计分布函数（CDF，Cumulative Distribution Function）， μ \mu μ是平均数（mean）， ϕ \phi ϕ 是正态分布的精度（方差的倒数）,例如 x ∼ N ( μ , ϕ ) x\sim \mathbf{N}(\mu, \phi) x∼N(μ,ϕ) ,概率密度函数(PDF，probability Density Function)是 f ( x ) = ϕ 2 π e − ϕ 2 ( x − μ ) 2 f(x) = \sqrt\frac{\phi}{2\pi}e^{-\frac{\phi}{2}(x-\mu)^2} f(x)=2πϕ e−2ϕ(x−μ)2 .对于一个有限的 L L L 要保证 ∑ l = 1 L ω l = 1 \sum_{l=1}^{L}\omega_l = 1 ∑l=1Lωl=1。当 L = ∞ L = \infty L=∞的时候， ∑ l = 1 ∞ ω l = 1 \sum_{l=1}^{\infty}\omega_l = 1 ∑l=1∞ωl=1是无疑的。

利用概率变换定义权值，建立实数 γ l ∈ ( − ∞ , + ∞ ) \gamma_l\in(-\infty,+\infty) γl∈(−∞,+∞) 到 u l ∈ ( 0 , 1 ) u_l\in(0,1) ul∈(0,1) 的映射。因此，转换（transformation）允许研究人员使用正态分布的潜在变量（latent variables） γ l \gamma_l γl 重新表述模型，通过第3.4节中介绍的数据增强Gibbs采样算法促进计算。这种转换还使得模型扩展能够包含额外的结构(例如，预测器)。此外，概率单位变换简化了先前第3.2节末尾所述的启发式（prior elicitation）。

为了让权重预测器 ω l （ x ) \omega_l（x) ωl（x) 是独立的，我们进一步解释了隐藏变量 γ l \gamma_l γl 是 X X X 的线性函数 { γ l ( X ) , X ∈ X } \{\gamma_l(X),X\in \mathcal{X}\} {γl(X),X∈X} 。(在本文中，我们使用上标作为指标而不是参数的指数):
ω l ( X ) = Φ ( γ l ( X ) ) ∏ p < l ( 1 − Φ ( γ p ( X ) ) ) , ( 3 ) \omega_l(X) = \Phi(\gamma_l(X))\prod_{p<l}(1-\Phi(\gamma_p(X))),(3) ωl(X)=Φ(γl(X))p<l∏(1−Φ(γp(X))),(3)
γ l ( X ) = θ l 1 + θ a 2 + θ r 3 + θ ( a , r ) 4 + θ m 5 + θ l e g 6 + θ ( a , l e g ) 7 + f 1 ( d e v s t a r t ∣ θ 8 ) + f 2 ( d u r ∣ θ 9 ) + f 3 ( l o g ( w g t ) ∣ θ 10 ) + f 4 ( l o g ( p c s ) ∣ θ 11 ) , ( 4 ) \gamma_l(X) = \theta^1_{l}+\theta_a^2+\theta_r^3+\theta_{(a,r)}^4+\theta^5_m+\theta^6_{leg}+\theta_{(a,leg)}^7\\+f_1(dev_{start}|\mathbf{\theta}^8)+f_2(dur|\mathbf{\theta}^9)\\+f_3(log(wgt)|\mathbf{\theta}^{10})+f_4(log(pcs)|\mathbf{\theta}^{11}),(4) γl(X)=θl1+θa2+θr3+θ(a,r)4+θm5+θleg6+θ(a,leg)7+f1(devstart∣θ8)+f2(dur∣θ9)+f3(log(wgt)∣θ10)+f4(log(pcs)∣θ11),(4)

其中， { θ l 1 } \{\theta_l^1\} {θl1} 控制着潜在类 l （ l = 1 , 2 , ⋯   , L ) l（l=1,2,\cdots,L) l（l=1,2,⋯,L) 的基线概率(baseline probability of latent class)； { θ a 2 } \{ \theta^2_a\} {θa2} 控制着航空公司 a ( a = 1 , 2 , ⋯   , 20 ) a(a = 1,2,\cdots,20) a(a=1,2,⋯,20) 的基线异质性（baseline heterogeneity）(我的理解是航空公司由很多因素所影响着)； { θ r 3 } \{\theta_r^3\} {θr3} 控制着航线 r ( r = 1 , 2 , ⋯   , 336 ) r(r=1,2,\cdots,336) r(r=1,2,⋯,336) 的异质性（heterogeneity）； θ ( a , r ) 4 \theta^4_{(a,r)} θ(a,r)4 表示权重（weight）对航空公司和航线之间可能的相互作用的依赖性； θ m 4 , θ l e g 5 , θ ( m , l e g ) 6 \theta_m^4,\theta_{leg}^5,\theta_{(m,leg)}^6 θm4,θleg5,θ(m,leg)6 的意义是相似的。此外， f 1 、 f 2 、 f 3 f_1、f_2、f_3 f1、f2、f3 和 f 4 f_4 f4 是样条函数(样条是一种特殊的函数，由多项式分段定义)被表示为4阶B样条的线性组合（参见附录第B.4节，了解所使用的样条的细节）,其中 d e v s t a r t dev_{start} devstart 的节点（knots）为 [ − 3 , − 2 , − 1 , 0 , 1 , 2 , 3 ] [−3,−2,−1,0,1,2,3] [−3,−2,−1,0,1,2,3] ， d u r dur dur 的节点(knots)为 [ 1 , 2 , 4 , 6 , 8 , 10 ] [1,2,4,6,8,10] [1,2,4,6,8,10] ， l o g （ w e i g h t ） log（weight） log（weight）的节点（knots）为 [ 2 , 4 , 6 , 8 ] [2,4,6,8] [2,4,6,8] ， l o g （ p c s ） log（pcs） log（pcs）的节点（knots）为 [ 1 , 3 , 5 ] [1,3,5] [1,3,5] 。这里，我们使用货物重量(wgt)和件数(pcs)的对数形式作为预测因子，因为原来的分布是高度倾斜(highly skewed)的。为了确保参数的识别（identification of parameters），我们让对于所有的 a , r a,r a,r 来说 θ 1 2 = θ 1 3 = θ ( a , 1 ) 4 = θ ( 1 , r ) 5 = θ 1 5 = θ ( 1 , l e g ) 6 = θ ( a , 1 ) 7 = 0 \theta_1^2=\theta_1^3=\theta_{(a,1)}^4=\theta_{(1,r)}^5=\theta_1^5=\theta_{(1,leg)}^6=\theta_{(a,1)}^7=0 θ12=θ13=θ(a,1)4=θ(1,r)5=θ15=θ(1,leg)6=θ(a,1)7=0。

后验计算（posterior computation）

在贝叶斯统计中，后验分布通常不具有解析性，涉及一个难以处理的归一化常数。因此后验计算通常依赖于在我们的运输风险应用中具有可以的准确性的大样本近似或马尔可夫链蒙特卡罗（MCMC）抽样。MCMC采样的基本思想是构造一个与模型参数的联合后验分布相对应的平稳分布的马尔可夫链，该方法避免了求解顽固性常数（intractable constant）的麻烦。为了使马尔可夫链具有适当的行为，需要谨慎地选择马尔可夫转移核（transition kernel），通常选择与 Metropolis-Hastings (MH) or Gibbs抽样相对应的样本。在参数较多的模型中，MH可以进行大量的调整，而Gibbs则通过从给定其他参数当前值的参数子集的条件后验分布中连续采样来避免调整。Gibbs抽样依赖于一种称为条件共轭的性质。以模型参数的一个子集为中心，以其他参数为条件，如果条件后验分布与前验分布形式相同，则先验概率分布是条件共轭的。我们的模型形式和先验分布的具体选择（将在下面描述）是由保留条件共轭性驱动的。

为了获得参数块的条件共轭，我们采用了一种常见的策略，即数据扩充（data augmentation）。数据扩充的基本思想是，通过谨慎地引入潜在变量，可以获得条件共轭性，从而得到一个简单的Gibbs抽样算法。然后对潜在变量和模型参数运行MCMC算法;尽管这增加了未知样本的数量，但通过允许模型参数直接从完全条件后验分布中选取块，可以提高效率。在频率统计中经常使用类似的扩充策略;例如，用EM算法拟合混合模型。这里，我们使用了Rodriguez等人的扩充策略。首先，我们关注的是 L < ∞ L<\infty L<∞的情况。对于每个观察值 y j ∣ X y_j|X yj∣X ,(如果有 n ( X ) n(X) n(X)个复制的话，在 X X X 的条件下对应复制 j ( j = 1 , 2 , ⋯   , n ( X ) ) j(j=1,2,\cdots,n(X)) j(j=1,2,⋯,n(X)) ,否则的话如果没有复制的话（比如 n ( X ) = 1 n(X)=1 n(X)=1 ）就去掉 j j j),我们引入一个潜在的指示剂变量（latent indicator variable） s j ( X ) s_j(X) sj(X) 例如当且仅当观察值（observation） y j ∣ X y_j|X yj∣X 从混合组成 l ( l = 1 , 2 , ⋯   , L ) l(l=1,2,\cdots,L) l(l=1,2,⋯,L) 取样时 s j ( X ) = l s_j(X) = l sj(X)=l 。这些潜在变量的使用是混合模型的标准。在潜在指标 s j ( X ) s_j(X) sj(X) 的帮助下，Gibbs对2000多个模型参数进行了抽样，并将其分类为四类。

分布的选择(Distributional Choices)

为了完成模型的指定，我们需要指定内核混合中的内核以及每个模型参数的先验概率分布。

正常核（Normal Kernel）

众所周知具有少量到中等分量的高斯混合模型（Gaussian mixture models with a small to moderate number of components）基本上可以近似任何单变量密度；此外，受正态分布（即条件共轭）计算可拓性(computational tractability)的影响，我们将PSBP混合模型的参数核 k ( ⋅ ∣ ψ ) k(\cdot|\psi) k(⋅∣ψ) 指定为正态分布 N ( μ , ϕ ) N(\mu,\phi) N(μ,ϕ) ,其中 ψ = ( μ , ϕ ) \mathbf{\psi} = (\mu,\phi) ψ=(μ,ϕ) 。回顾我们的混合模型采用式(1)中的形式，我们将上面方程中的核替换为normal Kernel，并使用 G X G_X GX 指定的先验PSBP。然后， y y y 的条件分布可以用简单的形式表示：
y ∣ X = ∑ l = 1 L ω l ( X ) N ( y ∣ μ l , ϕ l ) , ( 5 ) y|X = \sum_{l=1}^L\omega_l(X)N(y|\mu_l,\phi_l),(5) y∣X=l=1∑Lωl(X)N(y∣μl,ϕl),(5)
原子 { ( μ l , ϕ l ) , l = 1 , 2 , ⋯   , L } \{(\mu_l,\phi_l),l=1,2,\cdots,L\} {(μl,ϕl),l=1,2,⋯,L} 的先验是 N G ( ζ μ , ξ μ , a ϕ , b ϕ ) NG(\zeta_\mu,\xi_\mu,a_\phi,b_\phi) NG(ζμ,ξμ,aϕ,bϕ),一个条件共轭的正态gamma分布，例如
μ l ∼ N ( ζ μ , ξ μ , ϕ l ) , ϕ l ∼ G ( a ϕ , b ϕ ) \mu_l \sim \mathbf{N}(\zeta_\mu,\xi_\mu,\phi_l), \phi_l\sim \mathbf{G}(a_\phi,b_\phi) μl∼N(ζμ,ξμ,ϕl),ϕl∼G(aϕ,bϕ)
其中， l = 1 , 2 , ⋯   , L l=1,2,\cdots,L l=1,2,⋯,L。

权重上的参数的先验（prior for Parameters in weight）

我们选择正态先验作为参数 Θ = { { θ l 1 } , { θ a 2 } , { θ r 3 } , { θ ( a , r ) 4 } , { θ m 5 } , { θ l e g 6 } , { θ ( a , l e g ) 7 } , θ 8 , θ 9 , θ 10 , θ 11 } \Theta =\{\{\theta_l^1\},\{\theta_a^2\},\{\theta_r^3\},\{\theta_{(a,r)}^4\},\{\theta_m^5\},\{\theta_{leg}^6\},\{\theta_{(a,leg)}^7\},\mathbf{\theta}^8,\mathbf{\theta}^9,\mathbf{\theta}^{10},\mathbf{\theta}^{11}\} Θ={{θl1},{θa2},{θr3},{θ(a,r)4},{θm5},{θleg6},{θ(a,leg)7},θ8,θ9,θ10,θ11} , θ j i ∼ N ( v i , ϵ i ) \theta_j^i\sim N(v^i,\epsilon^i) θji∼N(vi,ϵi) 其中 i = 8 , ⋯   , 11 i=8,\cdots,11 i=8,⋯,11 并且 j = 1 , ⋯   , n ( i ) j=1,\cdots,n(i) j=1,⋯,n(i)

其中 n ( i ) n(i) n(i) 是用于预测的B样条基(B-spline bases)的数目。对于7个分类预测因子 θ 1 , ⋯   , θ 7 \mathbf{\theta}^1,\cdots,\theta^7 θ1,⋯,θ7 (比如 θ 1 = θ l 1 \theta^1={\theta^1_l} θ1=θl1 ) 的效率，我们构建了一个层次结构，允许在一个类别中的参数之间进行信息借用,
θ l 1 ∼ N ( Φ − 1 ( 1 L − l + 1 ) , ϵ 1 ) , ( 6 ) \theta_l^1 \sim N(\Phi^{-1}( \frac{1}{L-l+1}),\epsilon^1),(6) θl1∼N(Φ−1(L−l+11),ϵ1),(6)
θ a 2 ∼ N ( 0 , ϵ 2 ) , ⋯   , θ ( a , l e g ) 7 ∼ N ( 0 , ϵ 7 ) \theta_a^2 \sim N(0,\epsilon^2),\cdots, \theta_{(a,leg)}^7\sim N(0,\epsilon^7) θa2∼N(0,ϵ2),⋯,θ(a,leg)7∼N(0,ϵ7)
其中， ϵ i ∼ G ( c i , d i ) ( i = 1 , 2 , ⋯   , 7 ) \epsilon^i \sim G(c_i,d_i)(i=1,2,\cdots,7) ϵi∼G(ci,di)(i=1,2,⋯,7).在这里 G ( a , b ) G(a,b) G(a,b) 是一个Gamma分布，例如 x ∼ G ( a , b ) x\sim G(a,b) x∼G(a,b) ,这个的pdf是 f ( x ) = b a Γ ( a ) x a − 1 e − b x f(x) = \frac{b^a}{\Gamma(a)}x^{a-1}e^{-bx} f(x)=Γ(a)baxa−1e−bx。我们使用专门设计的 θ l 1 \theta_l^1 θl1 的先验概率来实施每个簇 l = 1 , 2 , ⋯   , L l= 1,2,\cdots,L l=1,2,⋯,L 的相同先验基线概率。

模型拟合评估

为了为我们的PSBP混合模型选择一组特定的预测因子，我们使用交叉验证比较了不同的可能性。我们选择了一个样本外（out-of-sample）预测性能最好的模型，如式所示：
γ l ( X ) = θ l 1 + θ a 2 + θ r 3 + θ ( a , r ) 4 + θ m 5 + θ l e g 6 + f 1 ( d e v s t a r t ∣ θ 8 ) + f 2 ( d u r ∣ θ 9 ) + f 3 ( l o g ( w g t ) ∣ θ 10 ) , ( 7 ) \gamma_l(X) = \theta_l^1+\theta_a^2+\theta_r^3+\theta_{(a,r)}^4+\theta_m^5+\theta_{leg}^6+f_1(dev_{start}|\theta^8)\\+f2(dur|\theta^9)+f_3(log(wgt)|\theta^{10}),(7) γl(X)=θl1+θa2+θr3+θ(a,r)4+θm5+θleg6+f1(devstart∣θ8)+f2(dur∣θ9)+f3(log(wgt)∣θ10),(7)
并将其用于模型比较和应用说明。 总的来说，贝叶斯方法避免了过度拟合，因为它对模型复杂度的隐式惩罚出现在后验分布中，而不是在最大似然估计中。 在混合模型中，由于使用了太多的混合组件，可能会导致过拟合。然而，PSBP 混合模型和相关的BNP模型自动地倾向于将概率权重的几乎可以忽略的部分放在一些组件上。这与我们观察到的相对于灵活频率回归模型的卓越的样本外预测性能一致。

我们使用了一个常见的贝叶斯策略，通过使用后验或预测图来进行有效检验。其基本思想是利用PSBP混合模型下的后验预测分布生成新的数据，如果模型拟合良好，生成的数据与经验数据相似。此外，我们还使用样本内和样本外预测残差将我们的模型与频率GAMs和灵活的混合模型进行了比较。总的来说，我们的模型表现出了优越的性能。

结果

附录C.1中的表C.1显示了（选定）模型参数的后验平均值和95%概率区间。我们从表中得出以下结论：

1、50个核的意思是 μ 1 , μ 2 , ⋯   , μ 50 \mu_1,\mu_2,\cdots,\mu_{50} μ1,μ2,⋯,μ50 范围在（-70.0-77.5）（小时）之间，标明模型预测的误差在（-3,3）天之间的连续的数据。50个核的标准差 1 ϕ 1 , 1 ϕ 2 , ⋯   , 1 ϕ 50 \frac{1}{\sqrt{\phi_1}},\frac{1}{\sqrt{\phi_2}},\cdots,\frac{1}{\sqrt{\phi_{50}}} ϕ1 1,ϕ2 1,⋯,ϕ50 1 的范围在（0.64,84.4）之间，这意味着普通的内核可能非常窄或平坦，允许进行灵活的估计。

2、水平参数 θ 1 1 , θ 2 1 , ⋯   , θ 49 1 \theta_1^1,\theta_2^1,\cdots,\theta_{49}^1 θ11,θ21,⋯,θ491 在−10.9到6.74之间变化，宽范围表明风险变化很大。例如，如果一条航线有 γ l − θ l 1 \gamma_l-\theta_l^1 γl−θl1 接近于零,那么对于比-5小的 θ l 1 \theta_l^1 θl1 中特定的 l l l ，权重 ∝ Φ ( γ l ( X ) ) ≈ Φ ( − 5 ) ≈ 0 \propto \Phi(\gamma_l(X)) \approx \Phi(-5)\approx 0 ∝Φ(γl(X))≈Φ(−5)≈0 ，从而消除了包含组件。通过类似的论证， θ l 1 \theta_l^1 θl1 也可以帮助确定 γ l ( X ) − θ l 1 \gamma_l(X)-\theta_l^1 γl(X)−θl1 中哪个组成 l l l 在其中起主要作用。

3、所有系数的后验分布比其先验分布更集中，这表明数据提供了更新先验的大量信息；此外，95%的概率区间很窄。

4、航空公司系数的后验估计 θ a 2 \theta_a^2 θa2 显示出很大的异质性(heterogeneity);这与大标准差 1 e 2 \frac{1}{\sqrt{e^2}} e2 1 一致。进一步的检查发现，除了A1的有效性被确定为零，其余19家航空公司95%的概率间隔中的18家不包括0。此外，其中许多都离零很远，这意味着对运输风险的影响很大。然而，根据线性模型，19家航空公司中只有2家与0有显著差异（使用0.05的p值阈值）。这种巨大差异归因于原始线性模型对分布平均值独立变量影响的估计，这在航空公司中相对稳定。PSBP检测到航空公司在运输风险分布的其他方面存在较大差异，例如不同模式的相对权重。

5、由于航线数量及其与航空公司的相互作用很大（分别为1336和587条），因此我们不在表C.1中包括它们的后汇总。然而，标准差 1 ϵ 3 \frac{1}{\sqrt{\epsilon^3}} ϵ3 1 的后验总结显示，各路径之间存在很大的异质性。此外，较大的标准差 1 ϵ 4 \frac{1}{\sqrt{\epsilon^4}} ϵ4 1 表明不同航空公司在同一航线上运输风险的分布存在很大差异。这表明，仔细选择航空公司可能会产生显著不同的航运体验。

总结和未来方向

利用国际航空货运物流的数据，我们研究了评估和预测运输风险的方法，并将其定义为实际到达时间和计划到达时间之间的偏差。为了适应数据的特殊多模态特性，我们引入了BNP模型PSBP混合模型，以灵活估计运输风险的条件密度函数。具体来说，我们在核心权重中构建一个线性结构，包括需求变量和决策变量，以便概率权重随预测器的变化而变化。如果数据允许，模型结构可以很容易地扩展到考虑其他因素，例如长期影响，允许系数随时间动态变化。 PSBP的优点包括通用性、灵活性、相对简单的采样算法和理论支持。

我们的结果表明，我们的方法比替代模型（原始线性模型、GAM和灵活混合模型）更准确地预测。此外，线性模型还可能导致误导性的推论。我们还演示了准确估计运输风险CPDF如何帮助托运人从多种可用服务中进行选择，以及帮助货代设定目标价格等。此外，我们还演示了如何使用该模型来估计预测因子（如航空公司）的基线性能。我们将研究结果与国际航空运输协会发布的业绩报告进行了比较，指出了国际航空运输协会简单的航空公司排名方法的不足之处。我们注意到，我们的方法的使用范围比这里所示的例子要广泛得多。实际上，任何涉及分布函数的决策都需要估计的CPDF。

我们研究中一个有趣的发现是，航空公司对交通运输分布形态的影响至关重要，而不是线性模型所关注的主题。对于未来的研究，我们希望能够提供关于为什么和航空公司在美国表现如此不同的信息。通过了解航空公司服务性能的根本驱动因素，我们可以为如何改进每个航空公司的服务质量提供指导，而不是简单地选择性能最好的航空公司。