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

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

🍪语言内容分比

🍇Python高斯过程

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

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

高斯模型

这样做似乎没有任何好处，因为正态分布本身并不是特别灵活的分布。然而，采用一组高斯分布（多元正态向量）有许多好处。首先，多元正态分布中任何元素子集的边际分布也是正态的：

此外，多元正态分布元素子集的条件分布（取决于剩余元素）也是正态的：

高斯过程将多元正态推广到无限维度。它被定义为随机变量的无限集合，任何边缘子集都具有高斯分布。因此，边缘化属性在其定义中是明确的。另一种思考无限向量的方式是将其视为函数。当我们编写一个以连续值作为输入的函数时，我们本质上是在暗示一个无限向量，它只在调用该函数时返回值（由输入索引）。同样，这种以函数表示的无限维高斯概念使我们能够以计算方式处理它们：我们永远不需要存储高斯过程的所有元素，只需根据需要计算它们即可。

从高斯过程采样

 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)

我们将利用多变量高斯分布的条件性质，逐点顺序生成实现。条件如下：

这是实现它的函数：

 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())

 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]

 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")

PreviousPython或R时偏移算法实现 NextPython气象辐射光谱能量平衡模型

Last updated 7 months ago

Was this helpful?

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())

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]