This causes an exception when I try to create a GP with 6-D inputs and 1-D output, e.g., X.shape=[100,6] and Y.shape=[100,1] In pymc3/gp/mean.py, the return values of Zero and Constant have the same shape as X: class Zero(Mean): def __call__(self, X): return tt.zeros(X.shape, dtype='float32') class Constant(Mean): def __init__(self, c=0): Mean.__init__(self) self.c = c def __call__(self, X): return tt.ones(X.shape) * self.c