Блог

Я хочу использовать приближение гауссовского процесса для простой одномерной тестовой функции, чтобы проиллюстрировать несколько вещей. Я хочу перебрать несколько разных значений для корреляционной матрицы (поскольку это 1D, это всего лишь одно значение) и показать, какое влияние оказывают различные значения на аппроксимацию. Насколько я понимаю, «тета» является параметром для этого. Поэтому я хочу установить значение theta вручную и не хочу никакой оптимизации / изменений в нем. Я думал, что постоянное ядро и функция 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 
# ...
 

Сейчас я сам запрограммировал простую реализацию, которая позволяет устанавливать корреляцию (здесь «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)