Mean functions: better defaults and a formula interface

[bwengals]

Most GP libraries treat the mean function as an afterthought. The default is usually zero, and most examples never change it. But for practitioners, the mean function matters:

• A zero mean function means the GP reverts to zero away from data, which is rarely what you want in practice
• A constant mean is a better default (the GP reverts to the data mean), but still leaves trend on the table
• A linear mean function captures obvious trends and lets the kernel focus on residual structure, which is the semi-parametric GP pattern that works well in practice

ptgp currently has a Zero mean function. Some things to consider:

Better default

Should the default mean function be a constant (estimated from data) rather than zero? This is a small change that would make out-of-the-box behavior much more reasonable for practitioners who don't think about mean functions.

Formula interface

A random idea that might be interesting: a Wilkinson-style formula language (like R's lm) for specifying mean functions:

gp = pg.VFE(kernel=k, mean=pg.mean.Formula("1 + x1 + C(x2)"), ...)

Open questions

• What should the default mean function be? Constant or linear?
• Do mean functions need their own class, or can they just be a Python callable?
• What do people actually use for mean functions in practice?

[jessegrabowski]

We might be able to use formulae and rip off some of what Bambi does. Tagging @tomicapretto for input since he's well versed in the subject, and might know sharp edges in proceeding this way.

[tomicapretto]

What kind of things do you want to put in the mean function? I'm not sure what are the cases you are considering here. But, if it just a bunch of additive literals that represent terms, maybe you can use your own thing?

If you need to handle more complex things, such as function calls, perhaps automatic resolution of redundancies and different kind of terms, you could go with formulae or formulaic.

If you care about stuff such as having all the resulting objects being sparse and/or being able to turn on/off the resolution of redundancies, formulaic is what you need IMO. If you need the pipe operator, as we do in Bambi, then go with formulae. I wish we had the ability to turn off automatic resolution of redundancies, but we don't have it yet.

Turning off the resolution of redundancies is important if you want to handle them in a more bayesian way (i.e. through a zero sum distribution rather than fixing a coefficient to zero or constraints of the same family).