How to model multi-input gaussian process?

Questions v5
1

I want to use GP within my graphical model, but I’m not sure on how to apply it. Assume that I have 5 input features, and 1 observed output i.e. y_t = \epsilon + \sum_{i=1}^{5} \beta_{i,t} * x_{i, t}. I want to model \beta_{i, t} using gaussian process. But when I follow this example, it gave me this error

ValueError: Incompatible Elemwise input shapes [(100, 100), (1, 5)]

The code snippet

import pymc as pm
import numpy as np
import matplotlib.pyplot as plt

n_size = 100
n_feat = 5

# Generate data
np.random.seed(1337)
x = np.cumsum(np.random.normal(0, 1, size=(n_size, n_feat)), axis=0)
beta = np.tanh(np.cumsum(np.tanh(np.random.normal(0, 0.1, size=(n_size, n_feat))), axis=0))
noise = np.random.normal(0, 0.5, size=(n_size, ))
y = (beta * x).sum(axis=-1) + noise

# Plots
fig, axes = plt.subplots(3, 1, figsize=(8, 9), sharex=True)
axes[0].plot(x, label=[f"x{i}" for i in range(1, n_feat+1)])
axes[0].set_title("input x")
axes[0].legend(loc="best")
axes[1].plot(beta, label=[f"beta{i}" for i in range(1, n_feat+1)])
axes[1].set_title("beta")
axes[1].legend(loc="best")
axes[2].plot(y, label="y")
axes[2].set_title("output y")
plt.show()

# Modeling
with pm.Model() as model:
    ell = pm.Gamma("ell", alpha=2, beta=1, shape=(n_feat,))
    eta = pm.HalfNormal("eta", sigma=5, shape=(n_feat,))
    cov = eta**2 * pm.gp.cov.ExpQuad(n_feat, ell)
    gp = pm.gp.Latent(cov_func=cov)
    f = gp.prior("f", X=x)
    sigma = pm.HalfNormal("sigma", sigma=1)
    pm.Normal("likelihood", mu=f.sum(axis=1), sigma=sigma, observed=y)
    
    model_trace = pm.sample(
        nuts_sampler="numpyro",
        draws=1_000,
        tune=1_000,
        chains=4,
        idata_kwargs={"log_likelihood": True},
    )
    model_posterior_predictive = pm.sample_posterior_predictive(
        trace=model_trace
    )

image
image683×757 89.8 KB