# (C) Petteri Pajunen 2006 # # # gpreg.R # # gaussian processes in regression samples <- 50 noisestd <- 0.1 cat("Input vertical scaling: (1.0)") a <- scan(n=1) # 2.0 cat("Input length scaling: (0.5)") b <- scan(n=1) #0.5 cat("Input noise std: (0.05)") c <- scan(n=1) #0.1 x <- runif(samples,min=-1.0,max=1.0) y <- x^3+rnorm(samples,mean=0.0,sd=noisestd) X11(width=5,height=5) plot(x,y) # 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 <- a*exp(-(x-t[it])^2/(2*(b)^2)) gamma <- a+c gp[it] <- t(k) %*% iC1 %*% y gpvar[it] <- gamma-t(k) %*% iC1 %*% k } lines(t,gp,lwd=2) lines(t,gp+sqrt(gpvar),col="red") lines(t,gp-sqrt(gpvar),col="red") scan(quiet=TRUE) dev.off()