The purpose of the mean function is to make both prediction and gp.Marginal work. The mean function basically acts as a super simple passthrough for the non-GP parts of the model. I think for your case it’d be easiest to either use the Linear mean function or write your own subclass in the same way as (3). The idea is that the predictors you’d change with pm.set_data, so X, are always the only input to __call__, and then you can stash anything else you need when you define the mean function in __init__, like alpha and beta. I think the Linear mean function might be what you need, or a subclass that’s very similar.