SKlearn Gaussian Process с постоянной, задаваемой вручную корреляцией

Я хочу использовать приближение гауссовского процесса для простой одномерной тестовой функции, чтобы проиллюстрировать несколько вещей. Я хочу перебрать несколько разных значений для корреляционной матрицы (поскольку это 1D, это всего лишь одно значение) и показать, какое влияние разные значения оказывают на приближение. Насколько я понимаю, тета является параметром для этого. Поэтому я хочу установить значение тета вручную и не хочу его оптимизации / изменений. Я думал, что постоянное ядро ​​и функция clone_with_theta могут дать мне то, что я хочу, но я не заставил это работать. Вот что у меня есть на данный момент:

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


def f(x):
"""The function to predict."""
return x/2 + ((1/10 + x)  * np.sin(5*x - 1))/(1 + x**2 * (np.sin(x - (1/2))**2))

# ----------------------------------------------------------------------
#  Data Points
X = np.atleast_2d(np.delete(np.linspace(-1,1, 7),4)).T
y = f(X).ravel()

# Instantiate a Gaussian Process model
kernel = ConstantKernel(constant_value=1, constant_value_bounds='fixed')

theta = np.array([0.5,0.5])
kernel = kernel.clone_with_theta(theta)

gp = GaussianProcessRegressor(kernel=kernel, optimizer=None)

# Fit to data using Maximum Likelihood Estimation of the parameters
gp.fit(X, y)

# Make the prediction on the meshed x-axis (ask for MSE as well)
y_pred, sigma = gp.predict(x, return_std=True)

# Plot 
# ...

python scikit-learn gaussian-process
person Greeco-K    schedule 21.11.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Я сам сейчас запрограммировал простую реализацию, которая позволяет установить корреляцию (здесь 'b') вручную:

import numpy as np
from numpy.linalg import inv
    
def f(x):
    """The function to predict."""
    return x/2 + ((1/10 + x)  * np.sin(5*x - 1))/(1 + x**2 * (np.sin(x - (1/2))**2))

def kriging_approx(x,xt,yt,b,mu,R_inv):
    
    N = yt.size
    one = np.matrix(np.ones((yt.size))).T
    
    r = np.zeros((N))
    for i in range(0,N):
        r[i]= np.exp(-b * (xt[i]-x)**2)
    

    y = mu + np.matmul(np.matmul(r.T,R_inv),yt - mu*one)
    y = y[0,0]
    return y

def calc_R (x,b):
    N = x.size
    # setup R
    R = np.zeros((N,N))
    for i in range(0,N):
        for j in range(0,N):
            R[i][j] = np.exp(-b * (x[i]-x[j])**2)
            
    R_inv = inv(R)
    return R, R_inv

def calc_mu_sig (yt, R_inv):
    N = yt.size
    one = np.matrix(np.ones((N))).T
    mu = np.matmul(np.matmul(one.T,R_inv),yt) / np.matmul(np.matmul(one.T,R_inv),one)
    mu = mu[0,0]
    
    sig2 = (np.matmul(np.matmul((yt - mu*one).T,R_inv),yt - mu*one))/(N)
    sig2 = sig2[0,0]
    
    return mu, sig2

# ----------------------------------------------------------------------

#  Data Points
xt = np.linspace(-1,1, 7)
yt = np.matrix((f(xt))).T

# Calc R   
R, R_inv = calc_R(xt, b)

# Calc mu and sigma   
mu_dach, sig_dach2 = calc_mu_sig(yt, R_inv)

# Point to get approximation for
x = 1

y_approx = kriging_approx(x, xt, yt, b, mu_dach, R_inv)
person Greeco-K    schedule 23.11.2020