Hierarchical Gaussian Process

Questions
1

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

2

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?