MultiOutput Gassian Processes

Using this tutorial as a guide: Multi-output Gaussian Processes: Coregionalization models using Hamadard product — PyMC example gallery - I have implemented this method to analyse common effects reported across longitudinal studies. Following this and plotting as per the original authors suggestion I get separate estimates of the trajectory of the param for each included study - Im wondering; how can I combine these outputs to estimate the common effect between studies?

Code is as follows:

def get_icm(input_dim, kernel, W=None, kappa=None, B=None, active_dims=None):
This function generates an ICM kernel from an input kernel and a Coregion kernel.
coreg =, W=W, kappa=kappa, B=B, active_dims=active_dims)
icm_cov = kernel * coreg # Use Hadamard Product for separate inputs
return icm_cov

with pm.Model() as model_contin:
# Priors
ell = pm.Gamma(“ell”, alpha=2, beta=0.5)
eta = pm.Gamma(“eta”, alpha=3, beta=1)
kernel = eta**2 *, ls=ell, active_dims=[0])
sigma = pm.HalfNormal(“sigma”, sigma=3)

# Get the ICM kernel
W = pm.Normal("W", mu=0, sigma=3, shape=(n_outputs, 2), initval=np.random.randn(n_outputs, 2))
kappa = pm.Gamma("kappa", alpha=1.5, beta=1, shape=n_outputs)
B = pm.Deterministic("B",, W.T) + pt.diag(kappa))
cov_icm = get_icm(input_dim=2, kernel=kernel, B=B, active_dims=[1])

# Define a Multi-output GP
mogp =
y_ = mogp.marginal_likelihood("f", X, Y, sigma=sigma)

gp.conditional was then used to estimate new x values however this returns each constitutent study individually.

Apologies if the ability to combine these preferably with shrinking is obvious, however its not obvious to me at the time of writing.

Are you thinking about a more hierarchical model? Like, there’s a common effect between studies that is some GP, and then the param for each study is drawn where that effect is the mean?

If so, you might want to try a hierarchical GP. In the ICM, the outputs have varying amounts of correlation, but there isn’t some global average effect. You’re learning a structure that’s more like, the outputs of study 1 and 2 are closely correlated, but not correlated with study 3. That’s what the matrix B is representing.


Thanks so much for your help. I am definitely thinking of a hierarchical structure. Would this be achieved using another GP to model the mean of the one below? Or would it use it for the covariance function? Apologies I am a bit unclear on this and can’t find many examples online.

I am definitely thinking of a hierarchical structure. Would this be achieved using another GP to model the mean of the one below?

That’s right. And it’s true that there isn’t a ton online. It’s still in progress, but example 3 in this PR has some working code you could try: Add HSGP reference and example nb by bwengals · Pull Request #647 · pymc-devs/pymc-examples · GitHub