python - 如何在 Scikit 中构建线性加性模型？

Question

我正在尝试建立一个模型，该模型采用简化的中间模型 f_I 的谓词并将其乘以某个系数 c_p 并添加一些一般的结果模型f_g:

然后优化系数 C 以适应数据。

我选择的通用模型是RBF模型，目前我已经对数据进行了RBF拟合:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import  RBF

# The true formula
beta = np.pi
omega = 10*np.pi
alpha = 1
phi = .5*np.pi

def f(t):
    return alpha*np.tanh(beta*t)*np.sin(omega*t + phi)

# Generate some test data
t = np.linspace(0, 1, 1000)

X_train = np.random.choice(t.reshape(-1), size=12).reshape(-1, 1)
y_train = f(X_train)

kernel = RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0))

gpr = GaussianProcessRegressor(kernel=kernel, random_state=1)
gpr.fit(X_train, y_train)

fig, ax = plt.subplots(dpi=200)
ax.plot(t, f(t), 'grey', label="True")
ax.plot(X_train, y_train, 'ko', label="Training")
ax.plot(t, gpr.predict(t.reshape(-1, 1)), 'b--', label="GLM")
fig.legend()

但是，我不知道如何实现这个线性加法模型。我试图为中间模型创建一个函数:

def IM(t):
    return t*np.sin(omega*t)

然后将 GaussianProcessRegressor 中的内核更改为如下内容:

kernel = IM()*ConstantKernel() + RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0))

但是，这不起作用，因为 IM() 需要是某种可与 Scikit 库一起使用的类。但是，我在网上找不到太多关于如何执行此操作的信息。

是否为中间模型制作自定义内核，f_I 是正确的方法吗？如果是这样，我如何构建自定义内核以使用 Scikit？

我正在关注以实现此目的的论文是 available here, in section 2.2

Answer 1

我认为拟合函数可以改善或纠正，但我认为它是一个很好的基础。

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel
from sklearn.base import RegressorMixin


class MyModel(RegressorMixin):
    def __init__(self, model_fi, model_f_g):
        self.model_fi = model_fi
        self.coef_c_p = ConstantKernel()
        self.model_f_g = model_f_g

    def fit(self, X=None, y=None):
        self.model_fi = self.model_fi.fit(X, y)
        min_erreur = np.inf
        for i in np.arange(0.2, 1.0, 0.1):
            self.coef_c_p = ConstantKernel(constant_value=i)
            res = (self.model_fi.predict(X) * self.coef_c_p.constant_value) + self.model_f_g(X)
            erreur = self.score(y, res)
            if erreur < min_erreur:
                min_erreur = erreur
                best_value = i
        self.coef_c_p = ConstantKernel(constant_value=best_value)
        print(best_value)

    def predict(self, X=None):
        return (self.model_fi.predict(X) * self.coef_c_p.constant_value) + self.model_f_g(X)

# The true formula
beta = np.pi
omega = 10 * np.pi
alpha = 1
phi = .5 * np.pi


def f(t):
    return alpha * np.tanh(beta * t) * np.sin(omega * t + phi)


# Generate some test data
t = np.linspace(0, 1, 1000)

X_train = np.random.choice(t.reshape(-1), size=12).reshape(-1, 1)
y_train = f(X_train)

kernel = RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0))
gpr = GaussianProcessRegressor(kernel=kernel, random_state=1)


model = MyModel(gpr, RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0)))
model.fit(X_train, y_train)
print(model.predict(X_train))