Sequential Updating

Fixed the sequential updating algorithm (avoids numerical instabilities associated with matrix inversion). Currently defined as the function recursion:

for the covariance:

C_seq <- function(X, X_prime, i){
if(i <= 1)
cov(X, X_prime) - mmult(cov(X,x[i]), cov(x[i], X_prime)) / as.numeric( cov(x[i], x[i]) + sigma_n^2)
else
C_seq(X, X_prime,   i-1) - mmult(C_seq(X,x[i],   i-1), C_seq(x[i], X_prime,   i-1)) / as.numeric( C_seq(x[i], x[i],   i-1)  + sigma_n^2  )
}

for the mean:

mu_seq <- function(X, i){
if(i <= 1)
cov(x[i], X) * (y[i]-mu[i]) / as.numeric( cov(x[i], x[i]) + sigma_n^2)
else
mu_seq(X, i-1) + C_seq(x[i], X,  i-1) * (y[i]-mu[i]) / as.numeric( C_seq(x[i], x[i], i-1)  + sigma_n^2 )
}

Note that this needed a custom matrix multiplication to handle scalar by vector cases:

mmult <- function(x,y){
if(length(x) == 1){
x <- as.numeric(x)
x * y
} else if(length(y) == 1){
y <- as.numeric(y)
x * y
}  else
x %*% y
}

And probably needs to be written in C at least to have any hope of scaling. See implementation in sequential-method

Match is perfect (jittered to avoid overplotting the lazy way):