# Copyright (C) 2013 Philipp Benner
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
library("ggplot2")
library("gridExtra")
library("scales")
plot.gp.1d <- function(gp, x, main="",
xlabel=NULL, ylabel=NULL, alpha=0.5,
col.mean="red", plot.mean = TRUE, plot.variance=TRUE,
plot.scatter=TRUE, covariates=NULL, ...)
{
result <- summarize(gp, x, ...)
p <- ggplot(data.frame(x = x[,covariates], mean = result$mean), aes(x=x))
# add labels and title
p <- p + xlab(xlabel)
p <- p + ylab(ylabel)
p <- p + ggtitle(main)
# first add the scatter plot
if (plot.scatter && !is.null(gp$xp) && dim(gp$yp)[2] == 1) {
p <- p + geom_point(data=data.frame(x=gp$xp[,covariates], y=gp$yp),
aes(x = x, y = y))
}
# and the variance
if (plot.variance && !is.null(result$variance)) {
p <- p + geom_ribbon(data=data.frame(x=x[,covariates],
ymin=result$mean-2*sqrt(result$variance),
ymax=result$mean+2*sqrt(result$variance)),
aes(ymin=ymin, ymax=ymax),
col="gray", fill="black",
alpha=alpha)
}
# on top the regression line
if (plot.mean) {
p <- p + geom_line(aes(y=mean), col=col.mean)
}
return (p)
}
plot.gp.2d <- function(gp, x, plot.variance=TRUE, plot.scatter=FALSE,
low=muted("green"), mid="white", high=muted("red"),
midpoint=NULL,
covariates=NULL, ...)
{
# initialize all plot objects
p1 <- NULL
p2 <- NULL
# compute expectation and variance
result <- summarize(gp, x, ...)
# maximum and minimum for plotting
limits <- c(min(result$mean), max(result$mean))
# midpoint for the color gradient
midpoint <- sum(limits)/2.0
# first, plot the expectation
p1 <- ggplot(data = data.frame(x = x[,covariates[1]], y = x[,covariates[2]], z = result$mean),
aes_string(x = "x", y = "y", z = "z")) +
geom_tile(aes_string(fill="z"), limits=limits) +
stat_contour() +
scale_fill_gradient2(limits=limits, low=low, mid=mid, high=high, midpoint=midpoint) +
ggtitle("Expected value")
if (plot.scatter && !is.null(gp$xp) && dim(gp$yp)[2] == 1) {
p1 <- p1 + geom_point(data=data.frame(x=gp$xp[,covariates[1]], y=gp$xp[,covariates[2]], z = gp$yp),
aes(x = x, y = y, colour = z))
}
# plot varience only if plot.variance is TRUE
if (plot.variance) {
p2 <- ggplot(data = data.frame(x = x[,covariates[1]], y = x[,covariates[2]], z = result$variance),
aes_string(x = "x", y = "y", z = "z")) +
geom_tile(aes_string(fill="z")) +
stat_contour() +
scale_fill_gradient(low=mid, high=high) +
ggtitle("Variance")
}
# use different grids depending on what plots are available
if (!is.null(p2)) {
grid.arrange(p1, p2, ncol=2)
}
else {
grid.arrange(p1, ncol=1)
}
# return a list of all plot objects
return (invisible(list(p1=p1, p2=p2)))
}
#' Plot a Gaussian process
#'
#' @param x Gaussian process
#' @param y positions where to evaluate the Gaussian process
#' @param covariates plot only the given covariates
#' @param ... arguments to be passed to methods
#' @method plot gp
#' @export
plot.gp <- function(x, y, covariates=NULL, ...)
{
# rename variables
gp <- x
x <- as.matrix(y)
if (is.null(covariates)) {
covariates <- 1:dim(gp)
}
if (length(covariates) == 1) {
plot.gp.1d(gp, x, covariates=covariates, ...)
}
else if (length(covariates) == 2) {
plot.gp.2d(gp, x, covariates=covariates, ...)
}
else {
stop("Gaussian process has invalid dimension.")
}
}
#' Plot a heteroscedastic Gaussian process
#'
#' @param x heteroscedastic Gaussian process
#' @param y positions where to evaluate the Gaussian process
#' @param alpha transparency of confidence intervals
#' @param ... arguments to be passed to methods
#' @method plot gp.heteroscedastic
#' @export
plot.gp.heteroscedastic <- function(x, y, alpha=0.3, ...)
{
# rename variables
model <- x
x <- as.matrix(y)
result <- summarize(model, x)
p <- plot(model$gp, x, ...)
p <- p + geom_ribbon(data=data.frame(x=x,
ymin=result$mean-2*sqrt(result$variance),
ymax=result$mean+2*sqrt(result$variance)),
aes(ymin=ymin, ymax=ymax),
alpha=alpha, fill="red")
return (p)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.