Gaussian Processes - combining constant features with covariance

That’s what I was thinking, but the Constant class doesn’t currently support that:

def __init__(self, c):
    super(Constant, self).__init__(1, None)
    self.c = c

https://github.com/pymc-devs/pymc3/blob/2cb6f79481f9130526e33cb679003d659868c4aa/pymc3/gp/cov.py#L155