Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
rls <- function(frml, lambda, data, k) {
# R build-in function for setting up linear regression model
mf <- model.frame(as.formula(frml), data, na.action = na.pass)
# The model output
y <- mf[ ,1]
# The design matrix
X <- model.matrix(as.formula(frml), mf)
# The number of observations
n <- nrow(X)
# The number of parameters
p <- ncol(X)
# Parameter matrix
Theta <- matrix(as.numeric(NA), nrow = n, ncol = p)
# The predictions
yhat <- as.numeric(rep(NA,n))
# The parameter vector
theta <- matrix(rep(0,p), ncol = 1)
# Start value for the parameter covariance P
P <- 10000 * diag(1, p)
# Use the inverse in the RLS
R <- solve(P)
# Iterate through and estimate the parameters
for (i in 1:(n-k)) {
x <- matrix(X[i, ])
# Check for NAs in inputs and output
if(all(!is.na(x)) & !is.na(y[i])){
# Update
R <- lambda * R + x %*% t(x)
theta <- theta + solve(R, x) %*% (y[i] - t(x) %*% theta)
Theta[i, ] <- t(theta)
}
# Predict
x <- matrix(X[i+k, ])
if(all(!is.na(x))){
yhat[i+k] <- t(x) %*% theta
}
}
# Return a list
L <- list()
L$residuals <- y - yhat
L$X <- X
L$y <- y
L$Theta <- Theta
return(L)
}