# (C) Petteri Pajunen 2006 # # # gphier.R # # gaussian processes in regression with priors on covariance parameters samples <- 50 noisestd <- 0.1 alphaa <- 2.0 betaa <- 5.0 alphab <- 0.5 betab <- 1.2 alphac <- 5.0 betac <- 1.0 # x <- runif(samples,min=-1.0,max=1.0) x <- matrix(runif(samples/2,min=-1.0,max=0.0),samples/2,1) x2 <- matrix(runif(samples/2,min=0.4,max=1.0),samples/2,1) x <- rbind(x,x2) y <- x^3+rnorm(samples,mean=0.0,sd=noisestd) X11(width=5,height=5) plot(x,y) # gradient ascent on the MAP estimate of covariance parameters stepscale <- 0.2 stepa <- 0.08 stepb <- 0.02 stepc <- 0.0005 a <- 1.0 b <- 1.0 c <- 1.0 Ca <- Cb <- Cc <- C <- matrix(0,samples,samples) # compute the covariance matrix C3 <- C2 <- C1 <- matrix(0,samples,samples) for (i in 1:samples) { for (j in 1:samples) { C1[i,j] <- a*exp(-(x[i]-x[j])^2/(2*(b)^2))+c*(i==j) } } iC1 <- solve(C1) t <- seq(-1.0,1.0,by=0.01) k <- matrix(0,samples,1) gp <- matrix(0,1,length(t)) gpvar <- gp for (it in 1:length(t)){ k <- as.vector(a)*exp(-(x-t[it])^2/(2*(as.vector(b))^2)) gamma <- a+c gp[it] <- t(k) %*% iC1 %*% y gpvar[it] <- gamma-t(k) %*% iC1 %*% k } lines(t,gp,lwd=1) # lines(t,gp+sqrt(gpvar),col="red") # lines(t,gp-sqrt(gpvar),col="red") cat("a,b, c: ") cat(format(c(a,b,c),digits=1)) cat("\n") for (iter in 1:4) { for (iter2 in 1:100) { # compute C dm <- x %*% matrix(1,1,length(x))-matrix(1,length(x),1) %*% t(x) dm <- dm^2 scaling <- as.double((2*b^2)^(-1)) C <- a*exp(-dm*scaling)+c*diag(length(x)) iC <- solve(C) # compute C' for a,b,c Ca <- exp(-dm*scaling) Cb <- C*b^(-3)*dm Cc <- diag(length(x)) # derivatives of the likelihood lla <- -0.5*sum(diag(iC %*% Ca))+0.5*t(y)%*% iC %*% Ca %*% iC %*% y llb <- -0.5*sum(diag(iC %*% Cb))+0.5*t(y)%*% iC %*% Cb %*% iC %*% y llc <- -0.5*sum(diag(iC %*% Cc))+0.5*t(y)%*% iC %*% Cc %*% iC %*% y # prior derivatives lpa <- -(alphaa+1)*a^(-1)+betaa*a^(-2) lpc <- -(alphac+1)*c^(-1)+betac*c^(-2) lpb <- (alphab-1)*b^(-1)-betab a <- as.double(a+stepscale*stepa*(lla+lpa)) b <- as.double(b+stepscale*stepb*(llb+lpb)) c <- as.double(c+stepscale*stepc*(llc+lpc)) } #cat("Gradient: ") #cat(lla+lpa,llb+lpb,llc+lpc) #cat("\n") cat("a,b, c: ") cat(format(c(a,b,c),digits=1)) cat("\n") # compute the covariance matrix C3 <- C2 <- C1 <- matrix(0,samples,samples) for (i in 1:samples) { for (j in 1:samples) { C1[i,j] <- a*exp(-(x[i]-x[j])^2/(2*(b)^2))+c*(i==j) } } iC1 <- solve(C1) t <- seq(-1.0,1.0,by=0.01) k <- matrix(0,samples,1) gp <- matrix(0,1,length(t)) gpvar <- gp for (it in 1:length(t)){ k <- as.vector(a)*exp(-(x-t[it])^2/(2*(as.vector(b))^2)) gamma <- a+c gp[it] <- t(k) %*% iC1 %*% y gpvar[it] <- gamma-t(k) %*% iC1 %*% k } lines(t,gp) } lines(t,gp+sqrt(gpvar),col="red") lines(t,gp-sqrt(gpvar),col="red") scan(quiet=TRUE) dev.off()