本文将介绍如何应用深度高斯过程建模量化信号中的不确定性
先进的机器学习 (ML) 技术能够从数据中得出的非常复杂的问题的解答。然而因为其“黑盒”的性质，很难评估这些答案的正确性。如果想在照片中找到特定的人或者物，例如在照片中找到猫的照片，这可能是很实用的。但在解决医疗数据时，因为可解释性的起因个别都不会被人们所承受，这导致 ML 模型在理论临床利用中的理论应用的概率很低。在这篇文章中，将介绍一种剖析生物数据的办法，它联合了古代 ML 的复杂性和经典统计办法的正当置信度评估。
本文的指标是提供
- 对生物数据应用高斯过程建模的合理性，特地是深度高斯过程 (DGP)，而不是一般的静止过程
- 哈密顿蒙特卡罗 (HMC) 办法概述及介绍其对 DGP 建模的帮忙
- 该办法在模仿数据示例中的利用阐明
剖析医疗数据的挑战
设想一下，医生正在监测一位承受过医治的患者。假如医治应该影响某种激素的程度。随着工夫的推移，咱们记录测试后果，并逐步失去这样的图表。
这显然是一个带有噪声的数据。咱们无奈立刻晓得产生了什么，但咱们很想晓得激素水平是否正在发生变化。个别状况下第一件事将数据拟合到线性回归模型，因为这是最简略的模型。然而在大多数事实世界的临床数据中，这简直不会失去给出任何信息。所以咱们要抉择一种更好的办法，比方将其建模为高斯过程（GP）为什么呢是通知过程呢？
首先，让咱们回顾一下什么是高斯过程（GP）。
将临床信号视为安稳高斯过程
当执行 GP 建模时，所有数据点都被认为是从多元高斯分布中提取的
这里有两点须要留神。1)这里的 K(X)是数据大小的方阵，蕴含非零的非对角元素; 这与将数据视为 n 次独立随机抽取不同。2) 数据点由 X 索引，例如在咱们的示例中，将程序定义为工夫 – 因而，过程就是咱们须要理解的对于 GP 的全部内容。
当提到拟合 GP 时，咱们的意思是找到矩阵 K(X) 的参数。为了使其更明确，让咱们将 GP 模型重写为
当初有了三个参数，在这里咱们提出一个论点：三个参数能够被解释为捕获尝试建模的过程的不同方面。
当咱们剖析随工夫变动的医学信号时，咱们像解决的问题是通过医学测试来评估某些生物在很多乐音的背景下进化的过程。这里能够将噪声源大抵分为两类。一个是测量噪声，随着现代医学的测试变得非常复杂，会在测量时产生不同水平的偏差，这个是无奈防止的。另一个起源是生物过程引起的变异，这与咱们感兴趣的过程无关并且更加简单，绝对的钻研也少得多。所以理论的信号（咱们试图检测的过程）可能是高度非线性的，这些都会造成数据的凌乱，所以个别不会对其进行准确建模，至多在可预感的将来是这样的。
然而咱们能够尝试在 GP 框架中对这三个参数进行半独立的建模。比如说取得一个最有可能的后果，两个噪声源的频率和振幅不同。
- g 参数是线性回归剖析中随机噪声的产物。从技术上讲它是 τ²g，这是“半独立”的起源，也很容易解释。将此参数与测量噪声分割起来是最天然的，只有咱们议论的是雷同类型的测量，这种类型的噪声不会因患者而异。它甚至能够通过反复测试并在进行每次训练时并保持一致。这里指定的生物噪声用 τ² 示意是因为该参数负责变动的幅度。这很可能取决于患者，但在同一患者的监测工夫内或多或少放弃不变。
- θ 负责形容数据点彼此之间的相干水平。如果咱们正在察看正在承受医治的患者，并且医治实际上正在发挥作用，那么它将在间断测量中反映为强相干。因而，θ 参数是为咱们试图检测的信号保留的。
GP 模型能够奇妙地满足咱们将一些数学实践放在医学数据上的需要。这是比线性回归的一个很好的提高。然而有一个问题，假如模型的所有参数放弃不变。这对于 g 和 τ² 参数可能很好，因为正如下面所探讨的，这两个参数归因于测量噪声和不受医治影响的身材生物变动。但 θ 的就不一样了，医治能够进行也能够扭转。这会影响信号的相关性，如果咱们想晓得这一点就必须思考 θ 的变动。
该解决方案带有一个深度高斯过程模型。就像深度神经网络一样，增加更多层。
两层 GP
有许多办法能够为 GP 增加层。这里咱们应用简略的办法：在工夫指标和测量信号之间插入另一个 GP
遵循援用 [2] 的办法，为什么这就足够了? 正如咱们下面提到的，咱们只须要适应 θ 的变动。通过为工夫变量引入额定的 GP，咱们以一种灵便的形式“扭曲”了测量工夫点之间的距离，从而产生了预期的成果。
但它也使拟合复杂化了! 咱们必须从这个后验散布中找到参数的最优值
这两种可能性被定义为
无关此表达式的推导，请参阅参考资料[2]。既然咱们写出了这个“奇怪”的公式，咱们就须要思考如何利用它做些有用的货色。最间接的办法是对这个散布进行数值抽样。最好的数值抽样办法是哈密顿蒙特卡洛法（HMC）!
哈密顿蒙特卡洛办法
咱们须要一种理智的办法来对的散布进行抽样，这有两个起因（实际上是同一件事件的两种观点）。咱们在多维空间中定义了概率分布。空间很大但概率是无限的。这意味着非零概率区域将被限度在空间的小体积中。但这个区域却是咱们想要失去的。所以咱们不能自觉地“扔飞镖”:1)浪费时间——咱们会在概率为零的地区破费大量工夫，这将是无用的;2)不精确的预计——在某些时候咱们将不得不进行采样，并且咱们有可能会错过非零区域。
本文的解决方案来自于将过程类比为物理学中家喻户晓粒子静止。设想一个带电粒子在相同电荷左近飞过空间，而你正在将一个球滑到一边。在所有这些状况下，静止的物体都被吸引并朝着锥体的底部或带相同电荷的粒子所在的中央挪动。但如果它一开始有足够的速度，例如你把球推到一边，它就会在这个吸引核心四周反弹。这张图显示了——两条可能的轨迹，它们具备不同的初始速度或动量，在物理学中，绿色的动量比红色的大，所以它在离核心更远的中央花更多的工夫。但他们俩都在朝着这个方向后退。
这种状况在实践力学中由哈密顿方程形容
这里
是动能和势能的总和。这个方程不仅简略而且看起来还十分的优雅。
应用 HMC
当初要做的就是概率分布是 U, 也就是势能。为了把它颠倒过去，须要在它后面加上一个减号，因为它曾经写在咱们的方程中，并通过求解具备不同动量值的哈密尔顿方程，摸索最大值左近的区域。HMC 的很多细节能够在参考文献 [3] 中找到。这是我为生成这样的轨迹而编写的示例代码。
def hamiltonian_monte_carlo(U, grad_U, current_position, steps=1, delta_t=0.5, change_step = False, masses = None):
q = np.copy(current_position)
p = np.random.normal(0.0,1.0,len(q))
current_momentum = np.copy(p)
if masses is None:
masses = np.ones(len(p))
q, p, path = leapfrog(q,p,grad_U,steps,delta_t,change_step, masses)
current_U = U(current_position)
current_K = sum(current_momentum*current_momentum/masses) / 2
proposed_U = U(q)
proposed_K = sum(p*p/masses) / 2
if (np.log(np.random.rand()) < (current_U-proposed_U+current_K-proposed_K)):
return (q,proposed_U,1,path) # accept
else:
return (current_position, current_U,0,None) # reject正如咱们看到这里有一个小问题——须要提供概率函数 grad_U 的导数。
不可能针对任何参数得出对数似然导数的解析表达式。然而能够尝试一种预计梯度的数值方所以咱们应用 tensorflow 库进行测试，上面是应用 tensorflow 函数实现对数似然的代码片段（代码看起来非常简单）。对于导数，它应用 GradientTape，这个函数为咱们发明了奇观，因为它是可用的，这缩小了咱们很多的工作。
def log_likelihood(self,params): # -log(p)
params = tf.convert_to_tensor(params, dtype=tf.float64)
noise_matrix = tf.math.multiply(tf.eye(num_rows =
self.n,dtype=tf.float64),params[0])
...
def log_likelihood_grad(self,params):
params = tf.convert_to_tensor(params, dtype=tf.float64)
with tf.GradientTape() as tape:
tape.watch(params)
y = self.log_likelihood(params)
return tape.gradient(y,params)它是如何工作的?
首先要看下咱们的数据，因为咱们演示的数据是生成的，模仿办法如下：
num_points = 20
target_mean = np.zeros(num_points)
target_mean[8] = 4.
target_mean[9] = 4.
target_mean[10] = 4.
target_cov = np.eye(num_points)
for i in range(num_points-1):
target_cov[i+1,i] = target_cov[i,i+1] = 0.1
data_sample_y = 2.*np.random.multivariate_normal(target_mean,target_cov,size=60).T
noise = np.random.normal(0.,1.,num_points)
data_sample_x = np.linspace(0,100,num_points)
data_to_model = data_sample_y[:,30] + noise它来自一个 20 维的多元高斯函数，除两头 3 个地位的均值为 4 外，其余地位的均值都为 0，第一个非对角元素的协方差矩阵都为 0.1。再把整个数据乘以 2（没有其余的非凡起因只是为了好玩）。这里生成了 60 个样本。随机挑出一个，而后加上方差为 1 的随机正态噪声。如下图所示，蓝色的线是筛选的样本，蓝色的点是增加了噪声的样本，灰色的虚线是来自雷同散布的残余样本，黑线是它们的平均值。
本文中用于概念验证的模仿数据中的所有样本。蓝点是用于拟合的数据，灰色虚线是雷同散布的类似样本，黑线是代表这些样本的平均值信号。
灰色线条是为了给咱们一个来自这个散布的数据的不确定性的视觉感官。它还能够作为咱们办法的一个额定性能，将试图依据给定的一个样原本预计其不确定性。当然次要的指标是预计黑线——信号。
咱们为一个带有噪声的样本 (上图中的蓝点) 增加了 HMC 的两层 GP 实现，并得出了以下后果。这里有 24 个参数(τ²，g，θy，θy，和 20 Ws)，所以不可能显示整个空间。下图只显示了其中两个参数。右边的两个 θs 是成心从远离正确值的中央开始的，轨迹向一个静止点挪动的速度很快，当达到了可能的概率最大值，就开始围绕它回旋。左边的图显示了总能量 vs 迭代次数。
应用参数的抽样后验值推断了 200 个样本，后果如下图所示，其中蓝线和黑线别离代表原始样本和实在数据的平均值。红线是推断样本的均值。灰色虚线是所有要与上图模仿数据进行比拟的单个推断样本。
最初，这就是咱们的最终后果。彩色——咱们试图复原的信号，红色——咱们对它的最佳猜想，灰色虚线——咱们猜想的 95% 置信区间。
这些置信区间可能比实在的要宽一点，是因为明确检测到了两头信号的变动，咱们不会被两边的摆动所蛊惑。思考到咱们的样本很少（一个样本有多种变异起源，没有反复），这是一个十分弱小的后果，对吧？
援用
[1] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, (2006), The MIT Press ISBN 0–262–18253-X.
[2] A. Sauer, R. B. Gramacy, and D. Higdon, (2021) Active Learning for Deep Gaussian Process Surrogates, arXiv:2012.08015v2
[3] M. Betancourt, A Conceptual Introduction to Hamiltonian Monte Carlo, (2018), arXiv:1701.02434
https://avoid.overfit.cn/post/be3e5d168e854c0daafc2845aad6d754
作者：Svetlana Rakhmanova Shchegrova