# Hierarchical Gaussian Process

Hi,
I’ve seen this post but I was wondering if there is any update or example of hierarchical gp.

My code attempt:

``````with pm.Model() as model:
# hyper priors
i = 0
η_per_mu = pm.Uniform("η_per_mu", 0, 1)
η_med_mu = pm.Uniform("η_med_mu", 0, 1)

# yearly periodic component x long term trend
η_per = pm.HalfCauchy("η_per", beta=η_per_mu, testval=1.0, shape=2)

# small/medium term irregularities
η_med = pm.HalfCauchy("η_med", beta=η_med_mu, testval=0.1, shape=2)

for ch in [df_4, df_5]:
# yearly periodic component x long term trend
ℓ_pdecay = pm.Gamma("ℓ_pdecay_"+str(i), alpha=10, beta=0.075)
period = pm.Normal("period_"+str(i), mu=1, sigma=0.05)
ℓ_psmooth = pm.Gamma("ℓ_psmooth_"+str(i), alpha=4, beta=3)
cov_seasonal = (
η_per[i] ** 2 * pm.gp.cov.Periodic(1, period, ℓ_psmooth) * pm.gp.cov.Matern52(1, ℓ_pdecay)
)
gp_seasonal = pm.gp.Marginal(cov_func=cov_seasonal)

# small/medium term irregularities
ℓ_med = pm.Gamma("ℓ_med_"+str(i), alpha=2, beta=0.75)
α = pm.Gamma("α_"+str(i), alpha=5, beta=2)
cov_medium = η_med[i] ** 2 * pm.gp.cov.RatQuad(1, ℓ_med, α)
gp_medium = pm.gp.Marginal(cov_func=cov_medium)

# long term trend
η_trend = pm.HalfCauchy("η_trend_"+str(i), beta=2, testval=2.0)
ℓ_trend = pm.Gamma("ℓ_trend_"+str(i), alpha=4, beta=0.1)
cov_trend = η_trend ** 2 * pm.gp.cov.ExpQuad(1, ℓ_trend)
gp_trend = pm.gp.Marginal(cov_func=cov_trend)

# noise model
η_noise = pm.HalfNormal("η_noise_"+str(i), sigma=0.5, testval=0.05)
ℓ_noise = pm.Gamma("ℓ_noise_"+str(i), alpha=2, beta=4)
σ = pm.HalfNormal("σ_"+str(i), sigma=0.25, testval=0.05)
cov_noise = η_noise ** 2 * pm.gp.cov.Matern32(1, ℓ_noise) + pm.gp.cov.WhiteNoise(σ)

# The Gaussian process is a sum of these three components
gp = gp_seasonal + gp_medium + gp_trend

# Since the normal noise model and the GP are conjugates, we use `Marginal` with the `.marginal_likelihood` method
y_ = gp.marginal_likelihood("y_"+str(i), X=ch['t'].values[:, None], y=ch['y'].values, noise=cov_noise)

i += 1
# this line calls an optimizer to find the MAP
mp = pm.find_MAP(include_transformed=True)
``````

Any help would be highly appreciated

No, no specific support for multiple outputs or a hierarchical GP has been added (though it would be great if it was there!). Does your code here run fine?