Noise-level estimation with scikit-learn GPR package for multidimensional data

Question

I am trying to estimate the noise level for Gaussian Process. Scikit-learn gives an example on this website https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_noisy.html. The last chunk of code seems to do exactly what I need, but it is only for the 1-d case

from matplotlib.colors import LogNorm

length_scale = np.logspace(-2, 4, num=50)
noise_level = np.logspace(-2, 1, num=50)
length_scale_grid, noise_level_grid = np.meshgrid(length_scale, noise_level)

log_marginal_likelihood = [
    gpr.log_marginal_likelihood(theta=np.log([0.36, scale, noise]))
    for scale, noise in zip(length_scale_grid.ravel(), noise_level_grid.ravel())
]
log_marginal_likelihood = np.reshape(
    log_marginal_likelihood, newshape=noise_level_grid.shape
)
vmin, vmax = (-log_marginal_likelihood).min(), 50
level = np.around(np.logspace(np.log10(vmin), np.log10(vmax), num=50), decimals=1)
plt.contour(
    length_scale_grid,
    noise_level_grid,
    -log_marginal_likelihood,
    levels=level,
    norm=LogNorm(vmin=vmin, vmax=vmax),
)
plt.colorbar()
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Length-scale")
plt.ylabel("Noise-level")
plt.title("Log-marginal-likelihood")
plt.show()

However, I have inputs X which is two-dimensional, so I have two length scales to optimize. I modified the code as follows:

from matplotlib.colors import LogNorm
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
kernel = 1.0 * RBF(length_scale=np.array([1e-1,1e-1])) + WhiteKernel( 
    determination
    noise_level=1e-2, noise_level_bounds=(1e-10, 1e1)
)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.0)
gpr.fit(T, y)


res=gpr.get_params

length_scale1 = np.logspace(-2, 4, num=50)
length_scale2 = np.logspace(-2, 4, num=50)
noise_level = np.logspace(-2, 1, num=50)
length_scale_grid1,length_scale_grid2, noise_level_grid = np.meshgrid(length_scale1, length_scale2,noise_level)

log_marginal_likelihood = [
    gpr.log_marginal_likelihood(theta=np.log([0.36, scale1, scale2, noise]))
    for scale1, scale2, noise in zip(length_scale_grid1.ravel(),length_scale_grid2.ravel(), noise_level_grid.ravel())
]

log_marginal_likelihood = np.reshape(
    log_marginal_likelihood, newshape=noise_level_grid.shape
)

vmin, vmax = (-log_marginal_likelihood).min(), 50
level = np.around(np.logspace(np.log10(vmin), np.log10(vmax), num=50), decimals=1) 

The result seems wrong since it always gives the largest length scale and noise estimation. I looked up the source code of scikit-learn and it says theta : array-like of shape (n_kernel_params,)

So my questions are: 1) for 1d case, theta should be of length 2 (one for length scale and one for noise), but why is there a third valur 0.36 in np.log([0.36, scale, noise]))? 2) is my modification correct? if not, how should I change my code for the 2d case?

Thank you very much!