PRML p152-p154

PRML.R

# 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]