Predict radio signal strength (RSS) using Gaussian Process Regression (GPR)

I want to use GPR to predict RSS from a deployed access point (AP). Since GPR gives mean RSS and its variance too, GPR could be very useful in positioning and navigation system. I read the GPR related published journals and got the theoretical insight of it. Now, I want to implement it with real data (RSS). In my system, the input and corresponding outputs (observations) are:

X: 2D cartesian coordinates points

y: an array of RSS (-dBm) at the corresponding coordinates

After searching online, I found that I can use sklearn software (using python). I installed sklearn and successfully tested the sample codes. The sample python scripts are for 1D GPR. Since my input sets are 2D coordinates, I wanted to modify the sample code. I found that other people have also tried to do the same, for example : How to correctly use scikit-learn's Gaussian Process for a 2D-inputs, 1D-output regression?, How to make a 2D Gaussian Process Using GPML (Matlab) for regression?, and Is kringing suitable for high dimensional regression problems?.

I am very new to Python and learned some syntax of python from this link. I modified the sample code as follows (I presume that there are some funny mistakes in my code):

print(__doc__)

import numpy as np
from matplotlib import pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

np.random.seed(1)

def f(x):
"""The function to predict."""
return x * y

# ----------------------------------------------------------------------
X = np.array([[0,0],[2,0],[4,0],[6,0],[8,0],[10,0],[12,0],[14,0],[16,0],[0,2],[2,2],[4,2],[6,2],[8,2],[10,2]
,[12,2],[14,2],[16,2]])

# Observations
y = np.array([-54,-60,-62,-64,-66,-68,-70,-72,-74,-60,-62,-64,-66,-68,-70,-72,-74,-76])

# Mesh the input space for evaluations of the real function, the prediction and
# its MSE
x1min=0;x1max=16;x2min=0;x2max=16
x1= np.linspace(x1min,x1max)
x2= np.linspace(x2min,x2max)
x=(np.array([x1,x2])).T

# Instanciate a Gaussian Process model
kernel = C(1.0, (1e-3, 1e3)) * RBF([5,5], (1e-2, 1e2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=15)

# 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, MSE = gp.predict(x, return_std=True)

print('Prediction:',y_pred,MSE)

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

While running the above lines of code, I got the following output:

None
Prediction: [-5.40000000e+01 -5.88802074e+01 -6.23801149e+01 -6.44250615e+01
-6.50311653e+01 -6.42850517e+01 -6.23227063e+01 -5.93118521e+01
-5.54395350e+01 -5.09041209e+01 -4.59094915e+01 -4.06591494e+01
-3.53489007e+01 -3.01581606e+01 -2.52410595e+01 -2.07190270e+01
-1.66763199e+01 -1.31592605e+01 -1.01790935e+01 -7.71766852e+00
-5.73476984e+00 -4.17587979e+00 -2.97937495e+00 -2.08249945e+00
-1.42580944e+00 -9.56070612e-01 -6.27784064e-01 -4.03618378e-01
-2.54057149e-01 -1.56553641e-01 -9.44384324e-02 -5.57673049e-02
-3.22368654e-02 -1.82417636e-02 -1.01045778e-02 -5.47896004e-03
-2.90797702e-03 -1.51068416e-03 -7.68099265e-04 -3.82199473e-04
-1.86106261e-04 -8.86757353e-05 -4.13434356e-05 -1.88610374e-05
-8.41972963e-06 -3.67822755e-06 -1.57267205e-06 -6.58205844e-07
-2.69700814e-07 -1.08211513e-07] [1.58637065e-04 2.39179570e+00 3.93253097e+00 4.50302334e+00
4.05844130e+00 2.63382155e+00 3.39817699e-01 2.65063741e+00
6.12264448e+00 9.84349987e+00 1.35871419e+01 1.71555182e+01
2.03942572e+01 2.32011479e+01 2.55271841e+01 2.73710856e+01
2.87690815e+01 2.97822154e+01 3.04834770e+01 3.09467005e+01
3.12385015e+01 3.14137041e+01 3.15139405e+01 3.15685809e+01
3.15969642e+01 3.16110178e+01 3.16176528e+01 3.16206408e+01
3.16219249e+01 3.16224516e+01 3.16226579e+01 3.16227351e+01
3.16227627e+01 3.16227722e+01 3.16227752e+01 3.16227762e+01
3.16227765e+01 3.16227766e+01 3.16227766e+01 3.16227766e+01
3.16227766e+01 3.16227766e+01 3.16227766e+01 3.16227766e+01
3.16227766e+01 3.16227766e+01 3.16227766e+01 3.16227766e+01
3.16227766e+01 3.16227766e+01]
[Finished in 2.538s]

The expected (predicted) values should be similar to y. The value I got is very different. The size of the testbed where I want to predict the RSS is 16*16 sq.meters. I want to predict RSS at every meter apart. I assume that the Gaussian Process predictor is trained with the Gaussian Decent algorithm in the sample code. I want to optimize the hyperparameter (theta: trained by using y and X) with Firefly algorithm.

In order to use my own data (2D input), in which line of code am I suppose to edit? Similarly, how can I implement Firefly algorithm (I've installed firefly algorithm using pip)?