# Need help in modeling 5D input, multi-output GP

**URL:** <https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190>\
**Category:** version agnostic\
**Tags:** modeling\
**Created:** [November 27, 2024, 10:32pm UTC](https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190 "2024-11-27T22:32:10Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![chiiingky](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/chiiingky/32/6323_2.png) [@chiiingky](https://discourse.pymc.io/u/chiiingky)\
**Post date:** [November 27, 2024, 10:32pm UTC](https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190/1 "2024-11-27T22:32:10Z")

</div>

Hi,

Let me start by saying I am very new to GP modeling.

I have adapted the codes from the [Multi-output Gaussian Processes Notebook](https://www.pymc.io/projects/examples/en/latest/gaussian_processes/MOGP-Coregion-Hadamard.html) to my own system. The ICM model in the notebook worked on my data for 1 batch of the data set. However, I need help in modifying the model when I want to train it on 6 batches of data. I’m guessing the modification has to do with the shapes of the parameters.

Each batch has an output of y\_species(t), for 11 time points and for 6 species.

Here is what my input array looks like for all 6 Batches:

```py
# First column
time = np.tile([1, 15, 30, 45, 60, 120, 180, 300, 420, 720, 1080], 36)

# Second, third, and fourth columns. Features related to each batch.
X1 = np.repeat([0.427, 0.427, 0.639, 0.249, 0.860, 0.427], 66)
X2 = np.repeat([0.389, 0.389, 0.389, 0.227, 0.782, 0.389], 66)
X3 = np.repeat([0.467, 0.542, 0.467, 0.273, 0.935, 0.427], 66)

# Fifth column
speciesIDX = np.tile(np.repeat([0, 1, 2, 3, 4, 5], 11), 6)

# Combine into a single array
input_array = np.column_stack((time, X1, X2, X3, speciesIDX))

```

The following model works when it is trained on the Batch 1 data (first 66 elements):

```py
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 = pm.gp.cov.Coregion(input_dim=input_dim, W=W, kappa=kappa, B=B, active_dims=active_dims)
    icm_cov = kernel * coreg # Use Hadamard Product for separate inputs
    return icm_cov

X = input_array[:66,:]
Y = output[:66]

with pm.Model() as model:
    # Priors
    ell = pm.Gamma("ell", alpha=5, beta=2)
    eta = pm.Gamma("eta", alpha=0.5, beta=0.1)
    kernel = eta**2 * pm.gp.cov.ExpQuad(input_dim=5, ls=ell, active_dims=[0])
    sigma = pm.HalfNormal("sigma", sigma=0.06)

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

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

```

How do I modify this model when I want to train it on 6 batches of data?

Many thanks,

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [November 30, 2024, 8:06pm UTC](https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190/2 "2024-11-30T20:06:10Z")

</div>

What do you mean by batch exactly? Is each batch independent, or is there some relation between them?

---

<div class="post-metadata">

**Author:** ![chiiingky](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/chiiingky/32/6323_2.png) [@chiiingky](https://discourse.pymc.io/u/chiiingky)\
**Post date:** [November 30, 2024, 8:28pm UTC](https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190/3 "2024-11-30T20:28:26Z")

</div>

Hi,

I have 6 batches of data. Each data is the concentration at a particular time: Ci(t) for i = 1, 2, …, 6 species, with 11 time points. I have successfully made a multi-output GP for 1 batch of data just by using the following input:

```py
time = np.tile([1, 15, 30, 45, 60, 120, 180, 300, 420, 720, 1080], 6) # for 6 species
speciesIDX = np.repeat([0, 1, 2, 3, 4, 5], 11) #species index
input_array = np.column_stack((time, speciesIDX))

```

Now, I have X1, X2, and X3 which has unique values per batch. You can think of them as additional features that describe the batch data. How do I combine this multi-ouput GP with a multi dimension input? I have been reading about the Kronecker product. Do you think I can use that here?

My ultimate goal is to have some prediction of the Ci(t) for i = 1, 2, …, 6 species, given X1, X2, X3.

---

<div class="post-metadata">

**Author:** ![chiiingky](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/chiiingky/32/6323_2.png) [@chiiingky](https://discourse.pymc.io/u/chiiingky)\
**Post date:** [December 4, 2024, 7:55pm UTC](https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190/4 "2024-12-04T19:55:53Z")

</div>

Update: I created a covariance per dimension and set the overall covariance as the overall product.

```auto
with pm.Model() as model:
    # Priors
    ell_t = pm.Gamma("ell", alpha=5, beta=2)
    ls1 = pm.Gamma("ls1", alpha=0.2, beta=0.2)
    ls2 = pm.Gamma("ls2", alpha=0.2, beta=0.2)
    ls3 = pm.Gamma("ls3", alpha=0.2, beta=0.2)
    
    eta = pm.Gamma("eta", alpha=0.5, beta=0.1)
    kernel = eta**2 * pm.gp.cov.ExpQuad(input_dim=5, ls=ell_t, active_dims=[0])
    
    # Get the ICM kernel
    W = pm.Normal("W", mu=0, sigma=1, shape=(6, 2), initval=np.random.randn(6, 2))
    kappa = pm.Gamma("kappa", alpha=3.0, beta=2.0, shape=6)
    B = pm.Deterministic("B", pt.dot(W, W.T) + pt.diag(kappa))
    coreg = pm.gp.cov.Coregion(input_dim=5, B=B, active_dims=[4])

    cov_x1 = pm.gp.cov.Matern52(5, ls=ls1, active_dims=[1]) 
    cov_x2 = pm.gp.cov.Matern52(5, ls=ls2, active_dims=[2])
    cov_x3 = pm.gp.cov.Matern52(5, ls=ls3, active_dims=[3])
    cov_product = kernel*cov_x1*cov_x2*cov_x3*coreg
    
    sigma = pm.HalfNormal("sigma", sigma=0.06)

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

```

I will report the results as soon as I get them. The sampling seems to take a while. Any suggestions on making this faster?

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [December 5, 2024, 9:21pm UTC](https://discourse.pymc.io/t/need-help-in-modeling-5d-input-multi-output-gp/16190/5 "2024-12-05T21:21:09Z")

</div>

> [@chiiingky](#):
>
> `pm.Gamma("ls1", alpha=0.2, beta=0.2)`

Sorry it took a while to get back to you. This is a tricky question. I think your model is going to be a little slow inherently, but I would definitely recommend moving away from those Gamma(0.2, 0.2) priors. I’d recommend using InverseGamma priors as per [this post](https://betanalpha.github.io/assets/case_studies/gaussian_processes.html). You can set the upper and lower bounds using `pm.find_constrained_prior` or the `maxent` funcion in Preliz so you don’t have to choose `alpha` and `beta` directly. The lengthscale priors can have a very big effect on sampling speed.