Skip to content

Instantly share code, notes, and snippets.

@edvakf edvakf/PRML.R
Created Nov 26, 2013

Embed
What would you like to do?
PRML p152-p154
# PRML p152-p154
library(ggplot2)
library(MASS)
# a0 and a1 are parameters to predict
a0 <- -0.3
a1 <- 0.5
f <- function(x) {a0 + a1 * x}
# prepare test data
N <- 20
prepare_test_set <- function() {
x <- runif(N, min=-1, max=1)
t <- f(x) + rnorm(N, mean=0, sd=0.2)
return (data.frame(x=x, t=t))
}
test_set <- prepare_test_set()
# plot test data and the line to generate the test data
p1 <- qplot(x, t, data=test_set) +
xlim(-1, 1) +
ylim(-1, 1)
p1 <- p1 + stat_function(fun = f)
#print(p1)
#par(ask=TRUE)
# want to find w0 and w1 that fit the "f" most likely
# prior probability distribution
M <- 100
w0 <- seq(-1,1,length=M)
w1 <- seq(-1,1,length=M)
multivariate_normal_distribution_2 <- function(x, mu, Sigma) {
# x: 1 x 2 column vector
# mu: 1 x 2 column vector
# Sigma: 2 x 2 matrix
return ((1/2/pi) / sqrt(det(Sigma)) * exp(-1/2 * t(x - mu) %*% solve(Sigma) %*% (x - mu)))
# x: n x 2 matrix
# mu: 1 x 2 column vector
# Sigma: 2 x 2 matrix
#return ((1/2/pi) / sqrt(det(Sigma)) * exp(-1/2 * t(x - mu) %*% solve(Sigma) %*% (x - mu)))
}
alpha <- 2.0
mu <- c(0,0)
Sigma <- rbind(
c(alpha, 0),
c(0, alpha)
)
density <- matrix(1:(M*M), nrow=M)
for (i in 1:M) {
for (j in 1:M) {
density[j,i] <- multivariate_normal_distribution_2( c(w0[i], w1[j]), mu, Sigma )
}
}
#prior$density <- multivariate_normal_distribution_2(mu, Sigma, rbind(prior$w0, prior$w1))
#filled.contour(w0, w1, density)
#par(ask=TRUE)
hypotheses <- mvrnorm(20, mu, Sigma)
p2 <- qplot(x, t, data=test_set) +
xlim(-1, 1) +
ylim(-1, 1)
for (i in 1:20) {
p2 <- p2 + geom_abline(slope=hypotheses[i,1], intercept=hypotheses[i,2], binwidth=0.1)
}
#print(p2)
Phi <- rbind()
w0 = hypotheses[1,1]
w1 = hypotheses[1,2]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.