This project is an example of Bayesian linear regression. To experiment more with this, please see my app, Bayesian Linear Regression Example.

First we generate a set of sample points a sine curve and apply some noise to them. These will be used to train the model.

# Generates a noisy sine curve.

N <- 10
sig2 <- 0.2

x <- as.matrix( runif(N, 0, 1) )

f <- function(x) { sin(2*pi*x) }

t <- f(x)  + sig2 * runif(N, -1, 1) 

sample_points <- data.frame(x = x, t = t)

sample_points %>% ggplot(aes(x = x, y = t)) + geom_point() + labs(x = "x", y = "y") +
      stat_function(fun = f, color = "green")

Now we will use these points to construct our model.

# First we define the basis functions
# with M (the number of basis functions) and s (the Gaussian width parameter).

M <- 12
dx <- 1/(M-1)
mu <- 0:(M-1)
mu <- mu / (M-1)

s <- 0.05

# Next we define the kernel we'll use.
phi <- function(x,mu,s) {
  exp(-(x-mu)^2/(2*s^2))
}

Phi <- matrix(nrow = N, ncol = M)

for (j in 1:M) Phi[,j] <- phi(x, mu[j],s)  

# Now we'll compute the regularized least squares solution with penalty lambda.

lambda <- 0.1

w <- solve( t(Phi) %*% Phi + lambda*diag(M) ) %*% t(Phi) %*%  t

# Next we compute the residual.

binv <- norm(t - Phi %*% w)^2/N

# Computes alpha from lambda and beta.

alpha <- lambda/binv

# Finally we compute the estimate at x (xs) and the variance of the uncertainty (sigN) at each x.

SN <- solve(alpha * diag(M) + t(Phi) %*% Phi / binv )
  
xs <- 0:100
xs <- xs / 100

tm <- matrix(nrow = 1, ncol = 101)

sigN <- matrix(nrow = 1, ncol = 101)

for (k in 1:length(xs) ) {
    
  p <- t( phi(xs[k],t(mu),s) )

  tm[k] <- t(w) %*% p
  sigN[k] <- binv + t(p) %*% SN %*% p
}

Finally we will plot our results. Note the gray error bars.

mod <- data.frame(xs = xs, tm = t( tm ), stdev = sqrt( t( sigN ) ) )

mod %>% ggplot(aes(x = xs, y = tm)) + geom_line(color = "red") + 
  geom_segment(aes(x = xs, xend = xs, y = tm + stdev, yend = tm - stdev), color = "gray") + 
  geom_point(data = sample_points, aes(x = x, y = t)) +
  labs(x = "x", y = "y") +
  stat_function(fun = f, color = "green")