Multi-dimensional gaussian process regression

I am using a multi dimensional gaussian process. The code is written in python. Note that ‘X_data’ is a list of 10 independent variables and ‘Y_data’ a dependent variable with a non-linear dependency on all of independent variables.
Y_data = fun_non_lin(X_data)

My questions are: 1. is the implementation correct?
2. Are there better implementations?
3. Where can i find necessary documentation for multidimensional gaussian process?
4. Here it is assumed that mean_func is a constant. In reality it’s a complex non-linear function. how can i introduce it?

def Gauss_pymc3(self, method = None):

    with pm.Model() as GP_Model:

        #Start here for Gaussian Process

        GP_const_func   =

        Sd_gp_len_scale = list(np.max(self.X_data, axis = 1)-
                               np.min(self.X_data, axis = 1))

        GP_len_scale    = []

        for i in range(len(Sd_gp_len_scale)):

                                            sigma = Sd_gp_len_scale[i]))

        GP_cov_func     =,ls= GP_len_scale)
        GP              = = GP_const_func, cov_func = GP_cov_func)

        noise_l             = pm.Gamma("noise_l", alpha=2, beta=2)
        GP_cov_func_noise   =, noise_l) \

        Y_data_gp       = GP.marginal_likelihood('Y_data_gp',
                                                 X = np.array(self.X_data).T,
                                                 y = self.Y_data,

        #GP_map          = pm.find_MAP()
        GP_trace        = pm.sample()