What is the difference between GPy.models.MultioutputGP and GPy.models.GPCoregioanlizedRegression?

Question

The GPy library given examples on multioutput GP given here. There also exists another example within the GPy.examples.Regression module here.

"""
compares GPCoregionalizedRegression with MultioutputGP
"""
import matplotlib.pyplot as plt
import numpy as np
import GPy


f = lambda x: np.sin(x) + 0.1 * (x - 2.0) ** 2 - 0.005 * x ** 3
fd = lambda x: np.cos(x) + 0.2 * (x - 2.0) - 0.015 * x ** 2
N = 10  # Number of observations
M = 10  # Number of derivative observations
Npred = 100  # Number of prediction points
sigma = 0.05  # Noise of observations
sigma_der = 0.05  # Noise of derivative observations
x = np.array([np.linspace(1, 10, N)]).T
y = f(x) + np.array(sigma * np.random.normal(0, 1, (N, 1)))

xd = np.array([np.linspace(2, 8, M)]).T
yd = fd(xd) + np.array(sigma_der * np.random.normal(0, 1, (M, 1)))

xpred = np.array([np.linspace(0, 11, Npred)]).T
ypred_true = f(xpred)
ydpred_true = fd(xpred)

# squared exponential kernel:
se = GPy.kern.RBF(input_dim=1, lengthscale=1.5, variance=0.2)
# We need to generate separate kernel for the derivative observations and give the created kernel as an input:
se_der = GPy.kern.DiffKern(se, 0)

# Then
gauss = GPy.likelihoods.Gaussian(variance=sigma ** 2)
gauss_der = GPy.likelihoods.Gaussian(variance=sigma_der ** 2)

# Then create the model, we give everything in lists, the order of the inputs indicates the order of the outputs
# Now we have the regular observations first and derivative observations second, meaning that the kernels and
# the likelihoods must follow the same order. Crosscovariances are automatically taken care of
m = GPy.models.MultioutputGP(
    X_list=[x, xd],
    Y_list=[y, yd],
    kernel_list=[se, se_der],
    likelihood_list=[gauss, gauss])


# Optimize the model
m.optimize(optimizer="lbfgs",messages=True,max_f_eval = 1000)

#%% CoregionalizedRegression on the same data, using the same kernel

K=GPy.kern.RBF(1, lengthscale=1.5, variance=0.2)
kicm = GPy.util.multioutput.ICM(input_dim=1,num_outputs=2,kernel=K)
mcoreg = GPy.models.GPCoregionalizedRegression(X_list = [x, xd], 
                                               Y_list = [y, yd], 
                                               kernel = kicm)
mcoreg.optimize(optimizer="lbfgs",messages=True,max_f_eval = 1000)

# # making predictions for the values:
# mu, var = m.predict_noiseless(Xnew=[xpred, np.empty((0, 1))])

# Two separate results!!!
m.param_array 
mcoreg.param_array

I wrote an MWE comparing the two using the module mentioned above. The model parameters obtained from the two models are different. I am unsure which one should be used, if a MultiOutput model is to be built. Or rather, when to use either?