Good machine learning algo for partial derivatives? Does anyone know of good robust algos to estimate partial derivatives of a regression model? I am talking about a general regression model like this:
$\mathbb{E}(y|x_1, x_2, ... x_n) = f(x_1, x_2, ... x_n)$, where I want estimates of $\frac{\partial f}{\partial x_k}$.
These are of great importance in medicine, economics, social sciences, etc. Now a linear model would give an approximation of these partial derivatives (its parameters) around the mean point. 
This approximation can be made valid further around the mean by the inclusion of polynomial terms, but the exponential increase in numbers of parameters as the polynomial degree rises makes this infeasible with a linear model (except perhaps with Lasso?).
For better local estimates over the whole input space, natural candidates would be kernel regression and feedforward neural networks, which have smooth forms (unlike tree-based methods) and from which partial derivatives are easy to recover. But I don't know about the robustness of these estimates, as the main aim of these models is to reduce prediction error, not find the right derivatives...
I know gradient boosting (and really all regression models) give a "partial dependance", but they are not partial derivatives: they are a function $f_k(x_k)$ that gives the mean value of the predicted regression function $\widehat{f}$ for a given value of $x_k$, averaging over the realisations of the other variables $(x_1, x_2 ... x_n)$ in the training dataset.
Any help would be greatly appreciated!
 A: I am not sure if this is what you are searching for, but I'll give it a try.
I think you can use Gaussian process regression and also get a nice robustness measure. 
Let's say you have $m$ observations $x_1, ..., x_m\in \mathbb R^n$ and a target variable $y_i \in \mathbb R$ for each $x_i$. For a Gaussian process with covariance function $k$, the joint density of a function value $f$ at an new point $x^*$ (the idea is that $y=f+\varepsilon$ where $\varepsilon$ is Gaussian noise) is given by
$$\left(\begin{array}{c}
\mathbf{f}\\
f^{*}
\end{array}\right)=\mathcal{N}\left(\mathbf{0},\left(\begin{array}{cc}
K+\sigma^{2}I & \mathbf{k}^{*}\\
\mathbf{k}^{*} & k\left(x^{*},x^{*}\right)
\end{array}\right)\right)$$
where $\mathcal N$ denotes the Gaussian distribution, $\mathbf K$ is the matrix of pairwise covariance function values between the $x_i$ and $\mathbf k^*$ is the vector of covariance function values between $x^*$ and all $x_i$ from the training set. 
From that it is easy to compute the posterior mean (i.e. the prediction at $x^*$) and covariance for $f^*$ (i.e. the measure of certainty which I think you mean with robustness). For all the formulae see the book by Rasmussen and Williams.
The interesting fact is that if $k(x_i,x_j)=\mbox{cov}(f_i,f_j)$, then $$\frac{\partial}{\partial x_{i\ell}}k(x_i,x_j)=\mbox{cov}\left(\frac{\partial}{\partial x_{i\ell}} f_i,f_j\right)$$ and $$\frac{\partial^2}{\partial x_{i\ell}\partial x_{jt}}k(x_i,x_j)=\mbox{cov}\left(\frac{\partial}{\partial x_{i\ell}} f_i,\frac{\partial}{\partial x_{jt}}f_j\right)$$ (again, check chapter 9.4 of Rasmussen and Williams). This means, you can just treat the partial derivatives as "unobserved function value $f^*$" and adapt the covariance functions in the above equation to the respective partial derivatives of the original covariance function. Then you simply predict the derivative like you would predict $f^*$. As a measure of certainty, you can look at the posterior covariance which is easily computes as well. 
There is a nice  Gaussian process toolbox for matlab written by Hannes Nickisch. It can probably help you play around with GPs.
