Search code examples
gaussiannormal-distribution

Conditioning Gaussian process regression


I was wondering how I could condition a Gaussian process (e.g., by adjusting a covariance function?) to produce only those trajectories that meet some pre-set constraints. As an example, the following code written in R produces smooth curves (Figure 1) and I'd like to make curves start from the origin, i.e, f(0)=0. Further, would it be possible to force additional characteristics, e.g., make all curves concave down and just one inflection point, etc.

require(MASS) ## mvrnorm function is needed 
## to generate covariance matrix using a squared exponential function
calc.sigma <- function( X1, X2, sigma_sq, phi_sq, tau_sq ) {
    Sigma <- matrix( rep( 0, length(X1)*length(X2) ), nrow=length(X1) )
    for( i in 1:nrow(Sigma) ) {
        for( j in 1:ncol(Sigma) ) {
            Sigma[i,j] <- sigma_sq*exp(-phi_sq*(X1[i] - X2[j])^2) + ifelse(i==j, tau_sq, 0.0);
        }
    }
    return(Sigma)
}

x <- seq(0,1,len=100) # input value
## parameters for covariance function 
sigma_sq <- 1 
phi_sq <- 10
tau_sq <- 1e-6

## parameters for mean function
a <- 1
b <- 5
mu <- a*( 1-exp(-b*x) ) # mean
# plot(mu)

sigma <- calc.sigma( x, x , sigma_sq, phi_sq, tau_sq )

N <- 10 # number of samples
samples <- matrix( rep( 0, length(x)*N ), ncol=N )
for ( i in 1:N ) {
    samples[,i] <- mvrnorm( 1, mu, sigma )
}

plot( x, samples[,1], type="l", xlim=c(0,1), ylim=c(-5,5) )
for( i in 2:ncol(samples) ){
    lines( x, samples[,i] ) 
}

Figure 1.


Solution

  • The first property that you are requesting is very easy to achieve:

    • the samplepaths of your GP will satisfy the property f(0)=0 almost surely if, and only if, the mean function and the variance function vanish at x=0.
    • Note that this is only possible for a non-stationary GP (or for the trivial GP with null mean and variance functions, of course).
    • A typical example of this the standard Brownian motion over [0; +\infty), which has a zero mean function and covariance function k(x,y) = min(x,y).
    • Starting from a generic covariance function k0, you can obtain this property using the conditioning formula: k(x,y) = k0(x,y) - k0(x,0) k0(0,y) / k0(0,0).

    All the above can be generalized to any property which is an equality constraint on some values of the function, or linear combination of such.

    Conditioning with respect to affine constraints is possible but harder:

    • Monotonicity of the samplepaths is equivalent to an infinite number of positivity constraints on the derivative (one at each point).
    • Concavity of the samplepaths is equivalent to an infinite number of positivity constraints on the second derivative.
    • Approximate versions of these properties can be enforced by conditioning a GP on a finite number of such constraints. The resulting process, however, is not a GP anymore. See, e.g., [1] for more details.
    • I don't know if there are any ready-to-use packages out there that propose this.

    [1] S. Da Veiga & A. Marrel (2012), Gaussian process modeling with inequality constraints, Annales de la faculté des sciences de Toulouse Mathématiques, 21(3), 529-555.