I want to perform coregionalized regression in GPy, however I am using a Bernoulli likelihood and then to estimate that as a Gaussian, I use Laplace inference. The code below shows how I would usually run a single-output GP with this set up (with my custom PjkRbf kernel):

likelihood = GPy.likelihoods.Bernoulli()
laplace_inf = GPy.inference.latent_function_inference.Laplace()
kernel = GPy.kern.PjkRbf(X.shape[1])
m = GPy.core.GP(X, Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf)

Now I am trying to run the same set up but as multi-output. This is something I have not been able to do.

I have tried using the GPCoregionalizedRegression class with an ICM kernel, as shown in the code below:

likelihood1 = GPy.likelihoods.Bernoulli()
likelihood2 = GPy.likelihoods.Bernoulli()
laplace_inf = GPy.inference.latent_function_inference.Laplace()
K = GPy.kern.PjkRbf(X.shape[1])
icm = GPy.util.multioutput.ICM(input_dim=2,num_outputs=2,kernel=K)
m = GPy.models.GPCoregionalizedRegression([X,X],[Y1,Y2],
kernel=icm,
likelihoods_list=[likelihood1, likelihood2])

Running this code throws an AssertionError, with a long stack-trace, but the last section shows the following. The likelihood cannot be asserted as a Gaussian.

~\Anaconda3\lib\site-packages\GPy\likelihoods\mixed_noise.py in gaussian_variance(self, Y_metadata)
22
23 def gaussian_variance(self, Y_metadata):
---> 24 assert all([isinstance(l, Gaussian) for l in self.likelihoods_list])
25 ind = Y_metadata['output_index'].flatten()
26 variance = np.zeros(ind.size)
AssertionError:

This is because I am **not able to pass the Laplace inference** to the GPCoregionalizedRegression model.

Can anyone offer advice on how I can resolve this, or if there is a different model I can use to perform multioutput regression with a Bernoulli likelihood and Laplace inference method?