# 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
#' Compute the Kullback-Leibler divergence between two models
#'
#' @param model0 first distribution
#' @param model1 second distribution
#' @param ... arguments to be passed to methods
#' @export
kl.divergence <- function(model0, model1, ...)
{
UseMethod("kl.divergence")
}
#' Compute the Kullback-Leibler divergence between two Gaussian
#' processes
#'
#' @param model0 first Gaussian process
#' @param model1 second Gaussian process
#' @param x positions at which the Gaussian processes are evaluated
#' @param ... unused
#' @method kl.divergence gp
#' @export
kl.divergence.gp <- function(model0, model1, x, ...)
{
gp0 <- model0
gp1 <- model1
list[mean0, variance0, covariance0] <- summarize(gp0, x, return.covariance=TRUE)
list[mean1, variance1, covariance1] <- summarize(gp1, x, return.covariance=TRUE)
stopifnot(!is.null(covariance0))
stopifnot(!is.null(covariance1))
n <- length(mean0)
L0inv <- solve(chol(covariance0))
L1inv <- solve(chol(covariance1))
tmp1 <- log(det((covariance1 %*% L0inv) %*% t(L0inv)))
tmp2 <- sum(diag(L1inv %*% (t(L1inv) %*% covariance0)))
tmp3 <- drop(t(mean0 - mean1) %*% (L1inv %*% (t(L1inv) %*% (mean0 - mean1))))
return (1/2*(tmp1 + tmp2 + tmp3 - n))
}
#' Compute the entropy of a model
#'
#' @param model for instance a Gaussian process
#' @param ... arguments to be passed to methods
#' @export
entropy <- function(model, ...)
{
UseMethod("entropy")
}
#' Compute the entropy of a Gaussian process
#'
#' @param model Gaussian process
#' @param x positions at which the Gaussian process is evaluated
#' @param ... unused
#' @method entropy gp
#' @export
entropy.gp <- function(model, x, ...)
{
gp <- model
list[mean, variance, covariance] <- summarize(gp, x, return.covariance=TRUE)
n <- length(mean)
return (1/2*(n*log(2*pi*exp(1)) + log(det(covariance))))
}
#' Compute the gain in information between two models
#'
#' @param model for instance a Gaussian process
#' @param ... arguments to be passed to methods
#' @export
information.gain <- function(model, ...)
{
UseMethod("information.gain")
}
#' Compute the gain in information between two Gaussian processes
#'
#' @param model first Gaussian process
#' @param gp2 second Gaussian process
#' @param x positions at which the Gaussian process is evaluated
#' @param ... unused
#' @method information.gain gp
#' @export
information.gain.gp <- function(model, gp2, x1, x2, ep=0.0, ...)
{
gp1 <- model
list[mean1, variance1, cov1] <- summarize(gp1, x1, return.covariance=TRUE)
list[mean2, variance2, cov2] <- summarize(gp2, x2, return.covariance=TRUE)
diag(cov1) <- diag(cov1) + ep
diag(cov2) <- diag(cov2) + ep
Linv <- solve(chol(cov2))
return (1/2*log(det(cov1 %*% Linv %*% t(Linv))))
}
# ------------------------------------------------------------------------------
if (FALSE) {
xp <- c(1, 2, 3)
yp <- c(0.7, 0.7, 0.7)
ep <- c(0.01, 0.01, 0.01)
gp0 <- new.gp(0.5, kernel.squared.exponential(1, 1))
gp1 <- new.gp(0.5, kernel.squared.exponential(1, 2))
gp2 <- posterior(gp1, xp, yp, ep)
kl.divergence(gp0, gp0, xp)
kl.divergence(gp0, gp1, xp)
kl.divergence(gp1, gp2, xp)
entropy(gp1, xp)
entropy(gp2, xp)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.