require(MASS) require(ggplot2) require(kernlab) require(nonparametricbayes) opts_knit$set(upload.fun = socialR::notebook.url) opts_chunk$set(dev.args=list(bg="transparent"), comment=NA, tidy=FALSE, message=FALSE) theme_update(plot.background = element_rect(fill = "transparent",colour = NA))
X <- seq(-5,5,len=50) obs <- data.frame(x = c(-4, -3, -1, 0, 2, 3), y = c(-2, 0, 1, 2, -1, -1)) l <- 1 sigma_n <- 0.1
#knit("beverton_holt_data.Rmd") #true <- data.frame(x = X, y = sapply(X, f, 0, p))
SE <- function(Xi,Xj, l=1) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2) cov <- function(X, Y) outer(X, Y, SE, l) n <- length(obs$x) K <- cov(obs$x, obs$x) I <- diag(1, n) L <- chol(K + sigma_n^2 * I) alpha <- forwardsolve(t(L), backsolve(L, obs$y)) # see http://stat.ethz.ch/R-manual/R-patched/library/base/html/chol2inv.html k_star <- cov(obs$x, X) Y <- t(k_star) %*% alpha v <- backsolve(L, k_star) Var <- cov(X,X) - t(v) %*% v loglik <- -.5 * t(obs$y) %*% alpha - sum(log(diag(L))) - n * log(2 * pi) / 2
d <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "direct") Ef <- d$Ef Cf <- d$Cf
s <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "sequential") ef <- s$Ef cf <- s$Cf
ggplot(data.frame(x=X, Ef=Ef, ef=ef))+ geom_point(aes(x,Ef), col='red') + geom_line(aes(x,ef)) ```` ### kernlab method ``` {r} k <- gp_fit(obs, X, c(sigma_n=sigma_n, l=l, tau=1), "kernlab") y_p <- k$Ef
Compare these results:
``` {r} require(reshape2) df <- data.frame(x = X, direct = Ef, Cholesky = Y, kernlab = y_p, sequential = ef) df <- melt(df, id = "x") ggplot(df)+ geom_jitter(aes(x, value, color = variable)) + geom_point(data = obs, aes(x,y)) ````
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.