如何将个体测量不确定性纳入高斯过程？

Question

我有一组观察结果，

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

Answer 1

正如 @amance 所建议的，您需要平方

errs

- 请参阅 这个 scikit 示例。然而，所得的均值和标准差似乎与蒙特卡罗结果无关，并且增加校准和路径的数量似乎并不能改善这种情况。我不知道如何解释这种差异。

现在更详细。我又添加了两个蒙特卡洛，将它们与修正后的

errs**2

方法 #1 进行比较（白噪声内核算作方法 #0）：方法 #2 是蒙特卡洛，内核中存在错误（就像 #1 一样） ，方法 #3 是蒙特卡洛方法，错误率为 0。你原来的蒙特卡洛是倒数第二个图。在最后一张图中，我比较了所有方法（#0 除外）的平均值和标准差 - 度量是相对欧几里德距离。

如您所见，方法 #1 和方法 #2 产生几乎相同的均值和标准差，并且方法 #3 的均值（0 个错误）仍然非常接近。方法 #4 产生非常不同的方法和标准开发，并且路径看起来非常不同 - 即使我注释掉错误，也会发生某种聚类。看起来要么是代码的实现有问题（只有 5 行，对我来说没有什么明显的问题），要么是模块有问题。 无论哪种方式，似乎方法#1就是你所需要的！

import matplotlib.pyplot as plt
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, 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. 
xs = np.array([44, 77, 125]) # measured inputs
fs = [8.64, 10.73, 12.13] # measured outputs
errs = np.array([0.1, 0.2, 0.3]) # uncertainty in measured outputs
finex = np.linspace(20, 200, 200) # inputs to predict
finex_T = [[f] for f in finex]
xs_T = [[x] for x in xs]
calibrations, paths = 100, 100
kernel = RBF(length_scale=9, length_scale_bounds=(10, 1e3))

#######################################
fig, ax = plt.subplots(6)
means, std_devs, titles, posterior_history = [None] * 5, [None] * 5, [None] * 5, [None] * 5
#######################################
### 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
gaussian_process = GaussianProcessRegressor(
  kernel=kernel + WhiteKernel(errs.mean(), noise_level_bounds=(errs.mean() - 1e-8, errs.mean() + 1e-8)), 
  n_restarts_optimizer=9, 
  normalize_y=True)
gaussian_process.fit(xs_T, fs)
means[0], std_devs[0] = gaussian_process.predict(finex_T, return_std=True)
titles[0] = 'White Noise Kernel - assumes uniform sensor error'

####################
### Approach 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

gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True, alpha=errs**2)
gaussian_process.fit(xs_T, fs)
means[1], std_devs[1] = gaussian_process.predict(finex_T, return_std=True)
titles[1] = 'uncertainty in likelihood - assumes measurements may be innacruate'

####################
### Test Approach with Errors:  
gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True, alpha=errs**2)
gaussian_process.fit(xs_T, fs)
posterior_history[2] = gaussian_process.sample_y(finex_T, calibrations * paths)
titles[2] = 'Monte Carlo + RBF Gaussian Process (only one calibration, with errors)'

####################
### Test Approach No Errors:  
gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True)
gaussian_process.fit(xs_T, fs)
posterior_history[3] = gaussian_process.sample_y(finex_T, calibrations * paths)
titles[3] = 'Monte Carlo + RBF Gaussian Process (only one calibration, no errors)'

####################
### Approach 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

posterior_history[4] = np.zeros((finex.size, calibrations * paths))
for sample in range(calibrations):
   gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, normalize_y=True)
   gaussian_process.fit(xs_T, fs)# + np.random.normal(0, errs**2*0))
   posterior_history[4][:, sample * paths : (sample + 1) * paths] = gaussian_process.sample_y(finex_T, paths)
titles[4] = 'Monte Carlo + RBF Gaussian Process'

for i in range(2, 5):
  means[i], std_devs[i] = posterior_history[i].mean(1), posterior_history[i].std(1)

####################
i_j = [[i, j] for i in range(2, 5) for j in range(1, i)]
means_err = [np.linalg.norm(means[i] - means[j]) / np.linalg.norm(means[i] + means[j]) for i, j in i_j]
str_dev_err = [np.linalg.norm(std_devs[i] - std_devs[j]) / np.linalg.norm(std_devs[i] + std_devs[j]) for i, j in i_j]

width = 0.35  # the width of the bars
x = np.arange(len(i_j)) 
rects1 = ax[-1].bar(x - width/2, means_err, width, label='means_err')
rects2 = ax[-1].bar(x + width/2, str_dev_err, width, label='str_dev_err')
ax[-1].set_ylabel('Discrepacies')
ax[-1].set_xticks(x, i_j)

####################
for i, ax_ in enumerate(ax[:-1]):
  if posterior_history[i] is not None:
    ax_.plot(finex, posterior_history[i], c='orange', alpha=0.005, zorder=-8)
  ax_.fill_between(finex, (means[i] - 1.96 * std_devs[i]), (means[i] + 1.96 * std_devs[i]), facecolor='grey')
  ax_.plot(finex, means[i], c='w')
  ax_.set_title(titles[i])
  ax_.errorbar(xs, fs, errs, linestyle="None", color="tab:blue", marker=".", markersize=8)
  ax_.set_ylim([6, 17])
  if i == len(ax) - 2:
    ax_.set_xlabel('input')
  else:
    ax_.set_xticks([])
  ax_.set_ylabel('output')

plt.show()

Answer 2

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**2)
gaussian_process.fit((np.atleast_2d(xs).T), (fs))
mu, std = gaussian_process.predict((np.atleast_2d(finex).T), return_std=True)