bounty将在7小时后过期。回答此问题可获得+300的声望奖励。kilojoules希望吸引更多人关注此问题。
我有一组观测值f_i=f(x_i)，我想构造一个概率替代值f(x) ~ N[mu(x), sigma(x)]，其中N是正态分布。每个观测输出f_i与一个测量不确定度sigma_i相关联。我想将这些测量不确定度合并到我的替代值f_i中。从而mu(x)预测观测值f_i(x_i)，并且预测的标准偏差sigma(x_i)包含观测输出epsilon_i中的不确定性。
我能想到的唯一方法是通过蒙特卡罗抽样和高斯过程建模的结合。用单个高斯过程完成这一点是理想的，而不需要蒙特卡罗样本，但我不能使这一工作。
我展示了实现我的目标的三种尝试。前两种方法避免了蒙特卡罗抽样，但没有预测f(x_i)的平均值，不确定性带包围了epsilon(x_i)。第三种方法使用蒙特卡罗抽样，并实现了我想要做的事情。
有没有一种方法可以创建一个高斯过程，平均预测平均观测输出，不确定性将包络观测输出中的不确定性，而不使用这种蒙特卡罗方法？
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, Matern, ExpSineSquared, WhiteKernel
# given a set of inputs, x_i, and corresponding outputs, f_i, I want to make a surrogate f(x).
# each f_i is measured with a different instrument, that has a different uncertainty.
# measured inputs
xs = np.array([44, 77, 125])
# measured outputs
fs = [8.64, 10.73, 12.13]
# uncertainty in measured outputs
errs = np.array([0.1, 0.2, 0.3])
# inputs to predict
finex = np.linspace(20, 200, 200)
#############
### approach 1: uncertainty in kernel
# - the kernel is constant and cannot change as a function of the input
# - uncertainty in measurements can be incorporated using a whitenoisekernel
# - the white noise uncertainty can be specified as the average of the observation error
# RBF + whitenoise kernel
kernel = 1 * RBF(length_scale=9, length_scale_bounds=(10, 1e3)) + WhiteKernel(errs.mean(), noise_level_bounds=(errs.mean() - 1e-8, errs.mean() + 1e-8))
gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True)
gaussian_process.fit((np.atleast_2d(xs).T), (fs))
mu, std = gaussian_process.predict((np.atleast_2d(finex).T), return_std=True)
plt.scatter(xs, fs, zorder=3, s=30)
plt.fill_between(finex, (mu - std), (mu + std), facecolor='grey')
plt.plot(finex, mu, c='w')
plt.errorbar(xs, fs, yerr=errs, ls='none')
plt.xlabel('input')
plt.ylabel('output')
plt.title('White Noise Kernel - assumes uniform sensor error')
plt.savefig('gp_whitenoise')
plt.clf()
####################
### Aproach 2: incorporate measurement uncertainty in the likelihood function
# - the likelihood function can be altered throught the alpha parameter
# - this assumes gaussian uncertainty in the measured input
kernel = 1 * RBF(length_scale=9, length_scale_bounds=(10, 1e3))
gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True, alpha=errs)
gaussian_process.fit((np.atleast_2d(xs).T), (fs))
mu, std = gaussian_process.predict((np.atleast_2d(finex).T), return_std=True)
plt.scatter(xs, fs, zorder=3, s=30)
plt.fill_between(finex, (mu - std), (mu + std), facecolor='grey')
plt.plot(finex, mu, c='w')
plt.errorbar(xs, fs, yerr=errs, ls='none')
plt.xlabel('input')
plt.ylabel('output')
plt.title('uncertainty in likelihood - assumes measurements may be innacruate')
plt.savefig('gp_alpha')
plt.clf()
####################
### Aproach 3: Monte Carlo of measurement uncertainty + GP
# - The Gaussian process represents uncertainty in creating the surrogate f(x)
# - The uncertainty in observed inputs can be propogated using Monte Carlo
# - downside: less computationally efficient, no analytic solution for mean or uncertainty
kernel = 1 * RBF(length_scale=9, length_scale_bounds=(10, 1e3))
posterior_history = np.zeros((finex.size, 100 * 50))
for sample in range(100):
simulatedSamples = fs + np.random.normal(0, errs)
gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True)
gaussian_process.fit((np.atleast_2d(xs).T), (simulatedSamples))
posterior_sample = gaussian_process.sample_y((np.atleast_2d(finex).T), 50)
plt.plot(finex, posterior_sample, c='orange', alpha=0.005)
posterior_history[:, sample * 50 : (sample + 1) * 50] = posterior_sample
plt.plot(finex, posterior_history.mean(1), c='w')
plt.fill_between(finex, posterior_history.mean(1) - posterior_history.std(1), posterior_history.mean(1) + posterior_history.std(1), facecolor='grey', alpha=1, zorder=5)
plt.scatter(xs, fs, zorder=6, s=30)
plt.errorbar(xs, fs, yerr=errs, ls='none', zorder=6)
plt.xlabel('input')
plt.ylabel('output')
plt.title('Monte Carlo + RBF Gaussian Process. Accurate but expensive.')
plt.savefig('gp_monteCarlo')
plt.clf()
2条答案
按热度按时间zpgglvta1#
使用第二种方法，只稍微改变Alpha
d4so4syb2#
正如@amance所建议的，你需要对
errs求平方--参见this scikit example。然而，得到的平均值和标准差似乎与你的蒙特-卡罗结果并不一致，而且增加校准和路径的数量似乎也不能改善这种情况。我不知道如何解释这种差异。现在更详细地介绍一下，我又添加了两个Monte-Carlos方法，将它们与修正后的
errs**2方法#1（白色噪声内核计数为方法#0）进行比较：方法#2是在核中有误差的蒙特-卡罗方法（就像方法#1一样），方法#3是0误差的蒙特卡罗方法。原始蒙特卡罗方法是从底部开始的第二个图。在最后一个图中，我比较了所有方法（除了方法#0）的平均值和标准差-测量值是相对欧几里得距离。如您所见，方法1和方法2产生的平均值和标准差几乎相同，方法3的平均值方法#4产生非常不同的平均值和标准偏差，而且路径看起来非常不同--存在某种聚集现象 * 即使我注解掉错误 *。（只有5行，对我来说没有什么明显的）或与模块。
无论哪种方式，似乎方法#1是你需要的！