Skip to content

Instantly share code, notes, and snippets.

@thierrymoudiki
Last active October 26, 2023 16:18
Show Gist options
  • Save thierrymoudiki/8f9502397461b2a70df6729591f29fc4 to your computer and use it in GitHub Desktop.
Save thierrymoudiki/8f9502397461b2a70df6729591f29fc4 to your computer and use it in GitHub Desktop.
Conformal Kernel Ridge Regression on small data
# rm(list=ls())
install.packages('learningmachine',
repos = c('https://techtonique.r-universe.dev',
'https://cloud.r-project.org'))
# or
# install.packages("remotes", repos = c(CRAN="https://cloud.r-project.org"))
#remotes::install_github("Techtonique/learningmachine")
library(learningmachine)
## Data -----------------------------------------------------------------------------
X <- as.matrix(mtcars[,-1])
y <- mtcars$mpg
set.seed(123)
(index_train <- base::sample.int(n = nrow(X),
size = floor(0.8*nrow(X)),
replace = FALSE))
X_train <- X[index_train, ]
y_train <- y[index_train]
X_test <- X[-index_train, ]
y_test <- y[-index_train]
dim(X_train)
dim(X_test)
## Kernel Ridge Regressor (KRR) object -----------------------------------------------------------------------------
obj <- learningmachine::KernelRidgeRegressor$new()
obj$get_type()
obj$get_name()
## Fit KRR -----------------------------------------------------------------------------
t0 <- proc.time()[3]
obj$fit(X_train, y_train, lambda = 0.05)
cat("Elapsed: ", proc.time()[3] - t0, "s \n")
## Test RMSE -----------------------------------------------------------------------------
print(sqrt(mean((obj$predict(X_test) - y_test)^2)))
## Graph ------------------------------------------------------------
res <- obj$predict(X = X_test, level = 95,
method = "splitconformal")
res2 <- obj$predict(X = X_test, level = 95,
method = "jackknifeplus")
par(mfrow=c(1, 2))
plot(c(y_train, res$preds), type='l',
main="split conformal",
ylab="",
ylim = c(3, 33))
lines(c(y_train, res$upper), col="gray60")
lines(c(y_train, res$lower), col="gray60")
lines(c(y_train, res$preds), col = "red")
lines(c(y_train, y_test), col = "blue")
plot(c(y_train, res2$preds), type='l',
main="jackknife +",
ylab="",
ylim = c(3, 33))
lines(c(y_train, res2$upper), col="gray60")
lines(c(y_train, res2$lower), col="gray60")
lines(c(y_train, res2$preds), col = "red")
lines(c(y_train, y_test), col = "blue")
# coverage rate
mean((y_test >= as.numeric(res$lower)) * (y_test <= as.numeric(res$upper)))
mean((y_test >= as.numeric(res2$lower)) * (y_test <= as.numeric(res2$upper)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment