Calculating the expression for the derivative of a Gaussian process I know based on the answers to this question Derivative of a Gaussian Process that the derivative of a Gaussian process is another Gaussian process, but I was wondering if someone could tell (or show) me how to explicitly calculate the expression for the derivative of a Gaussian process. For example, assuming that I have the following
$$f\sim GP(X\beta,\sigma^2R)$$
where $R$ is the Gaussian correlation function
$$R=\exp\left\{-\sum_{i=1}^n\frac{|x_{ij}-x_{ik}|^2}{\phi_i}\right\}$$
What would the distribution of $f'$ be (i.e., the derivative of $f$)?
 A: Since your domain for $f$ is $n$-dimensional you will actually have $n$ derivative processes $\frac{\partial f}{\partial x_j}$ with $j=1,\ldots n$. You need to calculate the mean and correlation function of $\frac{\partial f}{\partial x_j}$. From the linked answer you know that those are the derivative of $f$'s mean function and the derivative with respect to both arguments of the correlation $R$. So there is nothing more to do than to calculate those derivatives:
I assume that $\beta$ are constant, hence the mean function is $\frac{\partial X\beta}{\partial x_j}=\beta_j$. 
To keep confusion to a minimum let's write $R$ as $R(x,y)=exp\{-\sum_i\frac{(x_i - y_i)^2}{\phi_i}\}$. Then the derivative with respect to the first argument (assuming $\sigma$ is constant) is $$ \frac{\partial R}{\partial x_j}= \left(- 2 \frac{x_j-y_j}{\phi_j}\right) R(x,y)$$ and with respect to both $$ \frac{\partial }{\partial y_j}\frac{\partial R}{\partial x_j}=\frac{\partial }{\partial y_j}\left( R(x,y)\left(- 2 \frac{x_j-y_j}{\phi_j}\right) \right)=\frac{2}{\phi_j}R(x,y)\left( 1 - \frac{2}{\phi_j}(x_j - y_j)^2\right).$$ Which means in your notation that $$ \frac{\partial f}{\partial x_j} \sim \text{GP}\left( \beta_j, \sigma^2 \frac{2}{\phi_j} R(x,y)\left( 1 - \frac{2}{\phi_j}(x_j - y_j)^2\right)\right). $$
A: For simplicity I assume x has one dimension.
$\frac{\partial f}{\partial x}$ is normally distributed with expectation:
$E[\frac{\partial f}{\partial x}] = \frac{\partial}{\partial x}E[f]$
And Covariance:
$\text{Cov}(\frac{\partial f_1}{\partial x_1},\frac{\partial f_2}{\partial x_2})$ = $\frac{\partial^2 }{\partial x_2\partial x_1}\text{Cov}(f_1,f_2))$
In case you are using a gaussian correlation function: $\text{Cov}(f_1,f_2)) = \sigma^2\exp(-\frac{1}{2}\frac{(x_1-x_2)^2}{a^2})$, then:
$\text{Cov}(\frac{\partial f_1}{\partial x_1},\frac{\partial f_2}{\partial x_2})$ = $\frac{\sigma^2}{a^2}(1.0-\frac{(x_1-x_2)^2}{a^2})\exp(-\frac{1}{2}\frac{(x_1-x_2)^2}{a^2})$.
If x has more than one dimensions, each of the partial derivatives are normally distributed, and the covariance and expectation for each dimension can be calculated in the same way as above.
A: You can use sympy in Python, it will calculate any derivatives including integral defined one.
diffn(ff,x0,kk) :     
 dffk= Derivative(ff(x),x,kk)
 dffk1= simplify( dffk.doit())
 dffx0=  simplify(Subs(dffk1, (x), (x0)).doit())
 return dffx0

