Python空间地表联动贝叶斯地震风险计算模型

🎯要点

1. 使用贝叶斯推断模型兼顾路径和场地效应，量化传统地理统计曲线拟合技术。
2. 使用破裂和场地特征等地质信息以及事件间残差和事件内残差描述数学模型
3. 模型使用欧几里得距离度量、角距离度量和土壤差异性度量
4. 确定贝叶斯先验分布和后验分布参数，评估模型性能
   
   

🍇Python高斯过程

高斯过程是机器学习工具箱中的一个强大工具。它们使我们能够通过结合先验知识对数据进行预测。它们最明显的应用领域是用函数拟合数据。这称为回归，例如，用于机器人技术或时间序列预测。但高斯过程不仅限于回归——它们还可以扩展到分类和聚类任务。对于给定的一组训练点，可能有无数个函数可以拟合数据。高斯过程通过为每个函数分配一个概率为这个问题提供了一个优雅的解决方案。这个概率分布的平均值代表了数据最可能的特征。此外，使用概率方法我们可以将预测的置信度合并到回归结果中。

首先，我们将从连续视图转向函数的离散表示：我们感兴趣的不是寻找隐式函数，而是预测具体点处的函数值，我们将其称为测试点 $X$。那么，我们如何从迄今为止考虑过的多元正态分布中推导出这种函数视图呢？随机过程（例如高斯过程）本质上是一组随机变量。此外，这些随机变量中的每一个都有相应的索引 $i$。我们将使用此索引来指代我们的 $n$ 维多元分布的第 $i$ 维。

现在，高斯过程的目标是从训练数据中学习这种底层分布。相对于测试数据 $X$，我们将训练数据表示为 $Y$。高斯过程的关键思想是将 $X$ 的底层分布与 $Y$ 一起建模为多元正态分布。这意味着联合概率分布 $P_{X,Y}$ 跨越了我们想要预测的函数的可能函数值空间。请注意，测试和训练数据的联合分布具有 $|X|+|Y|$ 维度。

为了对训练数据进行回归分析，我们将这个问题视为贝叶斯推理。贝叶斯推理的基本思想是随着新信息的出现而更新当前假设。在高斯过程的情况下，此信息是训练数据。因此，我们感兴趣的是条件概率 $P_{X\mid Y}$。最后，我们回想一下，高斯分布在条件下是封闭的——因此 $P_{X\mid Y}$ 也是正态分布的。

现在我们已经了解了高斯过程的基本框架，只剩下一件事：如何建立这个分布并定义平均值 $\mu$ 和协方差矩阵 $\Sigma$ ？协方差矩阵 $\Sigma$ 由其协方差函数 $k$ 确定，后者通常也称为高斯过程的核。

在高斯过程中，我们将每个测试点视为一个随机变量。多元高斯分布的维数与随机变量的维数相同。由于我们想要预测 $|X|=N$ 个测试点处的函数值，因此相应的多元高斯分布也是 $N$ 维的。使用高斯过程进行预测最终归结为从该分布中抽取样本。然后，我们将结果向量的第 $i$ 个分量解释为与第 $i$ 个测试点相对应的函数值。

高斯模型

$$p(x\mid \pi ,\Sigma )=(2\pi {)}^{-k/2}|\Sigma {|}^{-1/2}\exp \left\{-\frac{1}{2}(x-\mu {)}^{\prime }{\Sigma }^{-1}(x-\mu )\right\}$$

这样做似乎没有任何好处，因为正态分布本身并不是特别灵活的分布。然而，采用一组高斯分布（多元正态向量）有许多好处。首先，多元正态分布中任何元素子集的边际分布也是正态的：
$$\begin{array}{c}p(x,y)=N\left(\left[\begin{array}{l}{\mu }_x\\ {\mu }_y\end{array}\right],\left[\begin{array}{cc}{\Sigma }_x& {\Sigma }_{xy}\\ {\Sigma }_{xy}^T& {\Sigma }_y\end{array}\right]\right)\\ p(x)=\int p(x,y)dy=N\left({\mu }_x,{\Sigma }_x\right)\end{array}$$
此外，多元正态分布元素子集的条件分布（取决于剩余元素）也是正态的：
$$p(x\mid y)=N\left({\mu }_x+{\Sigma }_{xy}{\Sigma }_y^{-1}\left(y-{\mu }_y\right),{\Sigma }_x-\Sigma xy{\Sigma }_y^{-1}\Sigma x{y}^T\right)$$
高斯过程将多元正态推广到无限维度。它被定义为随机变量的无限集合，任何边缘子集都具有高斯分布。因此，边缘化属性在其定义中是明确的。另一种思考无限向量的方式是将其视为函数。当我们编写一个以连续值作为输入的函数时，我们本质上是在暗示一个无限向量，它只在调用该函数时返回值（由输入索引）。同样，这种以函数表示的无限维高斯概念使我们能够以计算方式处理它们：我们永远不需要存储高斯过程的所有元素，只需根据需要计算它们即可。

从高斯过程采样

import numpy as np
def exponential_cov(x, y, params):
    return params[0] * np.exp( -0.5 * params[1] * np.subtract.outer(x, y)**2)

我们将利用多变量高斯分布的条件性质，逐点顺序生成实现。条件如下：
$$p(x\mid y)=N\left({\mu }_x+{\Sigma }_{xy}{\Sigma }_y^{-1}\left(y-{\mu }_y\right),{\Sigma }_x-\Sigma xy{\Sigma }_y^{-1}\Sigma x{y}^T\right)$$
这是实现它的函数：

def conditional(x_new, x, y, params):
    B = exponential_cov(x_new, x, params)
    C = exponential_cov(x, x, params)
    A = exponential_cov(x_new, x_new, params)
    mu = np.linalg.inv(C).dot(B.T).T.dot(y)
    sigma = A - B.dot(np.linalg.inv(C).dot(B.T))
    return(mu.squeeze(), sigma.squeeze())

我们将从高斯过程开始，其超参数为 ${\sigma }_{-}0=1,{\sigma }_{-}1=10$。我们还将假设零函数为平均值，因此我们可以绘制一个表示平均值一个标准差的带。

import matplotlib.pylab as plt
θ = [1, 10]
σ_0 = exponential_cov(0, 0, θ)
xpts = np.arange(-3, 3, step=0.01)
plt.errorbar(xpts, np.zeros(len(xpts)), yerr=σ_0, capsize=0)

让我们选择一个任意的起点进行采样，比如 x=1。由于没有先前的点，我们可以从无条件高斯中采样：

x = [1.]
y = [np.random.normal(scale=σ_0)]
print(y)

[0.4967141530112327]

现在，我们可以更新我们的置信带，给出我们刚刚采样的点，使用协方差函数生成新的逐点区间，条件是值 [$x_0,y_0$]。

def predict(x, data, kernel, params, sigma, t):
    k = [kernel(x, y, params) for y in data]
    Sinv = np.linalg.inv(sigma)
    y_pred = np.dot(k, Sinv).dot(t)
    sigma_new = kernel(x, x, params) - np.dot(k, Sinv).dot(k)
 return y_pred, sigma_new

x_pred = np.linspace(-3, 3, 1000)
predictions = [predict(i, x, exponential_cov, θ, σ_1, y) for i in x_pred]

y_pred, sigmas = np.transpose(predictions)
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
plt.plot(x, y, "ro")

因此，基于这一点以及我们指定的协方差结构，我们基本上限制了其他点的可能位置。现在让我们再取样一下：

m, s = conditional([-0.7], x, y, θ)
y2 = np.random.normal(m, s)
print(y2)

-0.1382640378102619

该点被添加到实现中，并可用于进一步更新下一个点的位置。

x.append(-0.7)
y.append(y2)
σ_2 = exponential_cov(x, x, θ)
predictions = [predict(i, x, exponential_cov, θ, σ_2, y) for i in x_pred]

y_pred, sigmas = np.transpose(predictions)
plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
plt.plot(x, y, "ro")

👉参阅、更新：计算思维 | 亚图跨际