I’m pretty new to Pymc3, and really, bayesian modelling so apologies ahead of time.
I’m interested in fitting a model which includes both a two-dimensional (spatial) gaussian process and linear coefficient for an indicator variable, with a large dataset.
I know I can do this using gp.Latent, but fitting takes a long time with large datasets because it is an O(N^3) operation.
e.g. something like
with pm.Model() as model: ℓ = pm.Gamma("ℓ", alpha=2, beta=1) η = pm.HalfCauchy("η", beta=5) cov = η ** 2 * pm.gp.cov.Matern52(1, ℓ) gp = pm.gp.Latent(cov_func=cov) f = gp.prior('f', X = X) beta_intercept = pm.Normal('intercept', mu = 0, sd = 10) beta_indicator = pm.Normal('beta', mu = 0, sd = 5) σ = pm.HalfCauchy("σ", beta=5) y_ = pm.Normal("y", mu = beta_intercept + f + beta_indicator * indicator, y=y, noise=σ)
I’m intrigued by the gp.MarginalSparse construct, and wondering if it would be possible to achieve something similar with it. My first thought was to design a custom mean function for gp.MarginalSparse which would depend on my indicator variable (separate mean for state = 1 vs 0). I don’t know if this is even possible or how I would code it. Would anybody be able to steer me in the right direction or provide an example?
Thanks in advance