How can I build a linear additive model in Scikit?

user_14720380 0
Viewed 267

I'm trying to build a model which takes the predication of a simplified intermediate-model f_I and multiplies it by some coefficient c_p and adds the result of some general model f_g:

Formula

The coefficients C are then optimized to fit the data.

The general model I am choosing is an RBF model, and so far I have fitted the RBF to the data:

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

RBF

However, I do not know how to implement this linear additive model. I have attempted to make a function for the intermediate model:

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

IM

And then change the kernel in the GaussianProcessRegressor to something like:

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

However, this doesn't work because the IM() needs to be a class of some sort that works with the Scikit library. However, I can’t find much information online on how to do this.

Is making a custom kernel for the intermediate model, f_I the right approach? If so, How do I build a custom kernel to work with Scikit?

The paper that I'm following to implement this is available here, in section 2.2

scikit-learn machine-learning python
2 Answers

I think the fit function can be ameliorated or corrected, but I think it's a good base.

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))
user_129862940
  • Utilize the inbuilt Sum and Product functions in gaussian_process.kernels
from sklearn.gaussian_process.kernels import RBF, Sum, Product, ConstantKernel
def IM(t):
    return t*np.sin(omega*t)
kernel_custom= Product(IM(X_train) ,ConstantKernel())
kernel = Sum(kernel_custom, RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0)))
  • Other option is to create Custom Kernels in RBF can be created as per the official documentation, sample code as shown here in line 1423 where you need to implement __call__ as StationaryKernelMixin, NormalizedKernelMixin already implements diag and is_stationary. Code added for reference
class RBF(StationaryKernelMixin, NormalizedKernelMixin, Kernel):

    def __call__(self, X, Y=None, eval_gradient=False):
        """Return the kernel k(X, Y) and optionally its gradient.
        Parameters
        ----------
        X : ndarray of shape (n_samples_X, n_features)
            Left argument of the returned kernel k(X, Y)
        Y : ndarray of shape (n_samples_Y, n_features), default=None
            Right argument of the returned kernel k(X, Y). If None, k(X, X)
            if evaluated instead.
        eval_gradient : bool, default=False
            Determines whether the gradient with respect to the log of
            the kernel hyperparameter is computed.
            Only supported when Y is None.
        Returns
        -------
        K : ndarray of shape (n_samples_X, n_samples_Y)
            Kernel k(X, Y)
        K_gradient : ndarray of shape (n_samples_X, n_samples_X, n_dims), \
                optional
            The gradient of the kernel k(X, X) with respect to the log of the
            hyperparameter of the kernel. Only returned when `eval_gradient`
            is True.
        """
user_10410850
Related