Splitting range of independent variable to maximize prediction within the subranges I have a dataset with two independent variables $X,Z \in \mathbb{R}$ and a dependent variable $Y \in \mathbb{R}$.
This dataset has the following characteristics:


*

*given some number $z$ and a "small" $\varepsilon > 0$, if we consider only the observations where $Z \in [z, z+\varepsilon]$ we find out that the relation $Y=f(X)$ in that subset is estimable quite easily with a polynomial function (it is a physical phenomenon).

*if we consider an adjacent interval in some sense (i.e. $[z,z+2\varepsilon]$ or $[z+\varepsilon,z+2\varepsilon]$) the previous estimate may either still be a good estimate or not with respect to some criteria.
I'm asked to "find a clusterization in $N$ adjacent intervals with respect to $Z$ where each interval is the largest possible interval where there exists a good fit $Y=f(X)$".
How should approach this problem? Is there any literature where problems like this have been studied? What would you recommend? I think this problem may be a "clustering of functionals" or something like that.
I hope I made myself clear, it is a bit difficult to explain, feel free to ask more details.
 A: Here's a go at the (very interesting!) problem.
If we can assume some kind of a distribution on $Y|X, Z, \theta$, for example a polynomial regression as you suggested, for example, $Y|X, Z, \theta \sim N(f_{\theta_Z}(X), \sigma_n^2)$, we obtain the conditional distribution $p(Y|X, Z, \theta)$.
To construct a full posterior, one needs the joint distribution $p(Y, X, Z, \theta)$. The remaining piece, considering the conditional distribution above is $p(X, Z, \theta)$.
Answering the other part of your question, fitting $Y=f(X)$ for different intervals of $Z$, is equivalent to finding the set of parameters $\theta(Z)$ for the range of $Z$ - which a GP could model beautifully.
You enforce a different $\theta$ for a prespecified number of intervals by treating $\theta$s and the interval bounds $[Z_{(i)}]$ as parameters. On the other hand, you could enforce a continuity in the function $\theta(Z)$ by using, for example, an RBF kernel.
This way, we'd have the distribution $p(\theta|Z)$ (a multivariate normal), and the "learning" could be done by using MCMC techniques (using Stan for example) or by simply finding the set of parameters that maximises the joint likelihood described here (using Tensorflow or any other optimiser).
After obtaining a fit, I guess that it'd be possible to use other parametric models to replace the GP once you're more confident about what the function $\theta(Z)$ looks like.
Side note: I've used GPs in vaguely similar ways before (to model "latent mappings" I guess I'd call them - which are essentially functionals), and I've had some great results. They're quite powerful when framed well.
