The model presented here is a frequentist – conformalized – version of the Bayesian one presented in #152. It is implemented in learningmachine
, both in Python and R, and is updated as new observations arrive, using Polyak averaging. Model explanations are given as sensitivity analyses.
1 - R version
%load_ext rpy2.ipython
%%R
utils::install.packages("bayesianrvfl", repos = c("https://techtonique.r-universe.dev", "https://cloud.r-project.org"))
utils::install.packages("learningmachine", repos = c("https://techtonique.r-universe.dev", "https://cloud.r-project.org"))
%%R
library(learningmachine)
X <- as.matrix(mtcars[,-1])
y <- mtcars$mpg
set.seed(123)
(index_train <- base::sample.int(n = nrow(X),
size = floor(0.6*nrow(X)),
replace = FALSE))
## [1] 31 15 19 14 3 10 18 22 11 5 20 29 23 30 9 28 8 27 7
X_train <- X[index_train, ]
y_train <- y[index_train]
X_test <- X[-index_train, ]
y_test <- y[-index_train]
dim(X_train)
## [1] 19 10
dim(X_test)
[1] 13 10
%%R
obj <- learningmachine::Regressor$new(method = "bayesianrvfl",
nb_hidden = 5L)
obj$get_type()
[1] "regression"
%%R
obj_GCV <- bayesianrvfl::fit_rvfl(x = X_train, y = y_train)
(best_lambda <- obj_GCV$lambda[which.min(obj_GCV$GCV)])
[1] 12.9155
%%R
t0 <- proc.time()[3]
obj$fit(X_train, y_train, reg_lambda = best_lambda)
cat("Elapsed: ", proc.time()[3] - t0, "s \n")
Elapsed: 0.01 s
%%R
previous_coefs <- drop(obj$model$coef)
newx <- X_test[1, ]
newy <- y_test[1]
new_X_test <- X_test[-1, ]
new_y_test <- y_test[-1]
t0 <- proc.time()[3]
obj$update(newx, newy, method = "polyak", alpha = 0.6)
cat("Elapsed: ", proc.time()[3] - t0, "s \n")
print(summary(previous_coefs))
print(summary(drop(obj$model$coef) - previous_coefs))
plot(drop(obj$model$coef) - previous_coefs, type='l')
abline(h = 0, lty=2, col="red")
Elapsed: 0.003 s
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.96778 -0.51401 -0.16335 -0.05234 0.31900 0.98482
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.065436 -0.002152 0.027994 0.015974 0.040033 0.058892
%%R
print(obj$summary(new_X_test, y=new_y_test, show_progress=FALSE))
$R_squared
[1] 0.6692541
$R_squared_adj
[1] -2.638205
$Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.5014 -2.2111 -0.5532 -0.3928 1.3495 3.9206
$Coverage_rate
[1] 100
$citests
estimate lower upper p-value signif
cyl -41.4815528 -43.6039915 -39.3591140 1.306085e-13 ***
disp -0.5937584 -0.7014857 -0.4860311 1.040246e-07 ***
hp -1.0226867 -1.2175471 -0.8278263 1.719172e-07 ***
drat 84.5859637 73.2987057 95.8732217 4.178658e-09 ***
wt -169.1047879 -189.5595154 -148.6500603 1.469605e-09 ***
qsec 22.3026258 15.1341951 29.4710566 2.772362e-05 ***
vs 113.3209911 88.3101728 138.3318093 7.599984e-07 ***
am 175.1639102 139.5755741 210.7522464 3.304560e-07 ***
gear 44.3270639 36.1456398 52.5084881 1.240722e-07 ***
carb -59.6511203 -69.8576126 -49.4446280 5.677270e-08 ***
$effects
── Data Summary ────────────────────────
Values
Name effects
Number of rows 12
Number of columns 10
_______________________
Column type frequency:
numeric 10
________________________
Group variables None
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable mean sd p0 p25 p50 p75 p100
1 cyl -41.5 3.34 -43.4 -43.4 -43.3 -41.7 -34.5
2 disp -0.594 0.170 -0.916 -0.635 -0.505 -0.505 -0.356
3 hp -1.02 0.307 -1.44 -1.40 -0.877 -0.768 -0.768
4 drat 84.6 17.8 59.5 76.7 89.5 89.5 128.
5 wt -169. 32.2 -204. -199. -166. -138. -138.
6 qsec 22.3 11.3 13.3 13.3 17.4 29.2 40.1
7 vs 113. 39.4 59.6 94.4 94.4 117. 191.
8 am 175. 56.0 124. 124. 153. 226. 245.
9 gear 44.3 12.9 26.3 38.7 47.9 47.9 76.0
10 carb -59.7 16.1 -77.3 -74.6 -58.2 -44.4 -44.4
hist
1 ▇▁▁▁▂
2 ▂▁▃▇▁
3 ▅▁▁▂▇
4 ▂▃▇▁▁
5 ▇▁▁▁▇
6 ▇▂▁▁▃
7 ▁▇▃▁▂
8 ▇▁▁▂▃
9 ▂▃▇▁▁
10 ▇▁▁▁▇
%%R
res <- obj$predict(X = new_X_test)
new_y_train <- c(y_train, newy)
plot(c(new_y_train, res$preds), type='l',
main="",
ylab="",
ylim = c(min(c(res$upper, res$lower, y)),
max(c(res$upper, res$lower, y))))
lines(c(new_y_train, res$upper), col="gray60")
lines(c(new_y_train, res$lower), col="gray60")
lines(c(new_y_train, res$preds), col = "red")
lines(c(new_y_train, new_y_test), col = "blue")
abline(v = length(y_train), lty=2, col="black")
%%R
newx <- X_test[2, ]
newy <- y_test[2]
new_X_test <- X_test[-c(1, 2), ]
new_y_test <- y_test[-c(1, 2)]
t0 <- proc.time()[3]
obj$update(newx, newy, method = "polyak", alpha = 0.9)
cat("Elapsed: ", proc.time()[3] - t0, "s \n")
print(obj$summary(new_X_test, y=new_y_test, show_progress=FALSE))
res <- obj$predict(X = new_X_test)
new_y_train <- c(y_train, y_test[c(1, 2)])
plot(c(new_y_train, res$preds), type='l',
main="",
ylab="",
ylim = c(min(c(res$upper, res$lower, y)),
max(c(res$upper, res$lower, y))))
lines(c(new_y_train, res$upper), col="gray60")
lines(c(new_y_train, res$lower), col="gray60")
lines(c(new_y_train, res$preds), col = "red")
lines(c(new_y_train, new_y_test), col = "blue")
abline(v = length(y_train), lty=2, col="black")
Elapsed: 0.003 s
$R_squared
[1] 0.6426871
$R_squared_adj
[1] -Inf
$Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.5686 -2.4084 -1.0397 -0.3897 1.5507 4.0215
$Coverage_rate
[1] 100
$citests
estimate lower upper p-value signif
cyl -42.1261096 -44.5327541 -39.7194651 2.932516e-12 ***
disp -0.6256505 -0.7347381 -0.5165629 1.613495e-07 ***
hp -1.0139634 -1.2198651 -0.8080617 6.747693e-07 ***
drat 82.8645391 74.8033348 90.9257434 5.680663e-10 ***
wt -170.7891742 -193.1932631 -148.3850853 1.053193e-08 ***
qsec 22.2365552 13.9564091 30.5167012 1.350094e-04 ***
vs 119.1784891 94.0163626 144.3406157 9.681321e-07 ***
am 174.2138307 134.1390652 214.2885963 2.127371e-06 ***
gear 42.7943293 36.9622907 48.6263678 1.523695e-08 ***
carb -59.4034661 -70.5135723 -48.2933599 3.127231e-07 ***
$effects
── Data Summary ────────────────────────
Values
Name effects
Number of rows 11
Number of columns 10
_______________________
Column type frequency:
numeric 10
________________________
Group variables None
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable mean sd p0 p25 p50 p75 p100
1 cyl -42.1 3.58 -44.1 -44.1 -44.1 -43.0 -34.9
2 disp -0.626 0.162 -0.933 -0.643 -0.514 -0.514 -0.514
3 hp -1.01 0.306 -1.47 -1.24 -0.787 -0.787 -0.787
4 drat 82.9 12.0 61.2 79.6 91.7 91.7 91.7
5 wt -171. 33.3 -210. -204. -142. -142. -142.
6 qsec 22.2 12.3 13.2 13.2 13.2 30.7 41.2
7 vs 119. 37.5 96.0 96.0 96.0 117. 193.
8 am 174. 59.7 123. 123. 123. 233. 247.
9 gear 42.8 8.68 27.1 40.4 49.2 49.2 49.2
10 carb -59.4 16.5 -78.8 -76.0 -45.1 -45.1 -45.1
hist
1 ▇▁▁▁▂
2 ▂▁▁▃▇
3 ▃▁▁▂▇
4 ▂▁▁▃▇
5 ▇▁▁▁▇
6 ▇▂▁▁▃
7 ▇▃▁▁▂
8 ▇▁▁▂▃
9 ▂▁▁▃▇
10 ▇▁▁▁▇
2 - Python version
!pip install git+https://github.com/Techtonique/learningmachine_python.git --verbose
import pandas as pd
import numpy as np
import warnings
import learningmachine as lm
# Load the mtcars dataset
data = pd.read_csv("https://raw.githubusercontent.com/plotly/datasets/master/mtcars.csv")
X = data.drop("mpg", axis=1).values
X = pd.DataFrame(X).iloc[:,1:]
X = X.astype(np.float16)
X.columns = ["cyl","disp","hp","drat","wt","qsec","vs","am","gear","carb"]
y = data["mpg"].values
display(X.describe())
display(X.head())
display(X.dtypes)
cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | |
---|---|---|---|---|---|---|---|---|---|---|
count | 32.000000 | 32.000000 | 32.0000 | 32.000000 | 32.000000 | 32.000000 | 32.000000 | 32.000000 | 32.000000 | 32.000000 |
mean | 6.187500 | 230.750000 | 146.7500 | 3.595703 | 3.216797 | 17.843750 | 0.437500 | 0.406250 | 3.687500 | 2.812500 |
std | 1.786133 | 123.937500 | 68.5625 | 0.534668 | 0.978516 | 1.788086 | 0.503906 | 0.499023 | 0.737793 | 1.615234 |
min | 4.000000 | 71.125000 | 52.0000 | 2.759766 | 1.512695 | 14.500000 | 0.000000 | 0.000000 | 3.000000 | 1.000000 |
25% | 4.000000 | 120.828125 | 96.5000 | 3.080078 | 2.580566 | 16.898438 | 0.000000 | 0.000000 | 3.000000 | 2.000000 |
50% | 6.000000 | 196.312500 | 123.0000 | 3.694336 | 3.325195 | 17.703125 | 0.000000 | 0.000000 | 4.000000 | 2.000000 |
75% | 8.000000 | 326.000000 | 180.0000 | 3.919922 | 3.610352 | 18.906250 | 1.000000 | 1.000000 | 4.000000 | 4.000000 |
max | 8.000000 | 472.000000 | 335.0000 | 4.929688 | 5.425781 | 22.906250 | 1.000000 | 1.000000 | 5.000000 | 8.000000 |
cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | |
---|---|---|---|---|---|---|---|---|---|---|
0 | 6.0 | 160.0 | 110.0 | 3.900391 | 2.619141 | 16.453125 | 0.0 | 1.0 | 4.0 | 4.0 |
1 | 6.0 | 160.0 | 110.0 | 3.900391 | 2.875000 | 17.015625 | 0.0 | 1.0 | 4.0 | 4.0 |
2 | 4.0 | 108.0 | 93.0 | 3.849609 | 2.320312 | 18.609375 | 1.0 | 1.0 | 4.0 | 1.0 |
3 | 6.0 | 258.0 | 110.0 | 3.080078 | 3.214844 | 19.437500 | 1.0 | 0.0 | 3.0 | 1.0 |
4 | 8.0 | 360.0 | 175.0 | 3.150391 | 3.439453 | 17.015625 | 0.0 | 0.0 | 3.0 | 2.0 |
0 | |
---|---|
cyl | float16 |
disp | float16 |
hp | float16 |
drat | float16 |
wt | float16 |
qsec | float16 |
vs | float16 |
am | float16 |
gear | float16 |
carb | float16 |
y.dtype
dtype('float64')
X
cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | |
---|---|---|---|---|---|---|---|---|---|---|
0 | 6.0 | 160.0000 | 110.0 | 3.900391 | 2.619141 | 16.453125 | 0.0 | 1.0 | 4.0 | 4.0 |
1 | 6.0 | 160.0000 | 110.0 | 3.900391 | 2.875000 | 17.015625 | 0.0 | 1.0 | 4.0 | 4.0 |
2 | 4.0 | 108.0000 | 93.0 | 3.849609 | 2.320312 | 18.609375 | 1.0 | 1.0 | 4.0 | 1.0 |
3 | 6.0 | 258.0000 | 110.0 | 3.080078 | 3.214844 | 19.437500 | 1.0 | 0.0 | 3.0 | 1.0 |
4 | 8.0 | 360.0000 | 175.0 | 3.150391 | 3.439453 | 17.015625 | 0.0 | 0.0 | 3.0 | 2.0 |
5 | 6.0 | 225.0000 | 105.0 | 2.759766 | 3.460938 | 20.218750 | 1.0 | 0.0 | 3.0 | 1.0 |
6 | 8.0 | 360.0000 | 245.0 | 3.210938 | 3.570312 | 15.843750 | 0.0 | 0.0 | 3.0 | 4.0 |
7 | 4.0 | 146.7500 | 62.0 | 3.689453 | 3.189453 | 20.000000 | 1.0 | 0.0 | 4.0 | 2.0 |
8 | 4.0 | 140.7500 | 95.0 | 3.919922 | 3.150391 | 22.906250 | 1.0 | 0.0 | 4.0 | 2.0 |
9 | 6.0 | 167.6250 | 123.0 | 3.919922 | 3.439453 | 18.296875 | 1.0 | 0.0 | 4.0 | 4.0 |
10 | 6.0 | 167.6250 | 123.0 | 3.919922 | 3.439453 | 18.906250 | 1.0 | 0.0 | 4.0 | 4.0 |
11 | 8.0 | 275.7500 | 180.0 | 3.070312 | 4.070312 | 17.406250 | 0.0 | 0.0 | 3.0 | 3.0 |
12 | 8.0 | 275.7500 | 180.0 | 3.070312 | 3.730469 | 17.593750 | 0.0 | 0.0 | 3.0 | 3.0 |
13 | 8.0 | 275.7500 | 180.0 | 3.070312 | 3.779297 | 18.000000 | 0.0 | 0.0 | 3.0 | 3.0 |
14 | 8.0 | 472.0000 | 205.0 | 2.929688 | 5.250000 | 17.984375 | 0.0 | 0.0 | 3.0 | 4.0 |
15 | 8.0 | 460.0000 | 215.0 | 3.000000 | 5.425781 | 17.812500 | 0.0 | 0.0 | 3.0 | 4.0 |
16 | 8.0 | 440.0000 | 230.0 | 3.230469 | 5.343750 | 17.421875 | 0.0 | 0.0 | 3.0 | 4.0 |
17 | 4.0 | 78.6875 | 66.0 | 4.078125 | 2.199219 | 19.468750 | 1.0 | 1.0 | 4.0 | 1.0 |
18 | 4.0 | 75.6875 | 52.0 | 4.929688 | 1.615234 | 18.515625 | 1.0 | 1.0 | 4.0 | 2.0 |
19 | 4.0 | 71.1250 | 65.0 | 4.218750 | 1.834961 | 19.906250 | 1.0 | 1.0 | 4.0 | 1.0 |
20 | 4.0 | 120.1250 | 97.0 | 3.699219 | 2.464844 | 20.015625 | 1.0 | 0.0 | 3.0 | 1.0 |
21 | 8.0 | 318.0000 | 150.0 | 2.759766 | 3.519531 | 16.875000 | 0.0 | 0.0 | 3.0 | 2.0 |
22 | 8.0 | 304.0000 | 150.0 | 3.150391 | 3.435547 | 17.296875 | 0.0 | 0.0 | 3.0 | 2.0 |
23 | 8.0 | 350.0000 | 245.0 | 3.730469 | 3.839844 | 15.406250 | 0.0 | 0.0 | 3.0 | 4.0 |
24 | 8.0 | 400.0000 | 175.0 | 3.080078 | 3.845703 | 17.046875 | 0.0 | 0.0 | 3.0 | 2.0 |
25 | 4.0 | 79.0000 | 66.0 | 4.078125 | 1.934570 | 18.906250 | 1.0 | 1.0 | 4.0 | 1.0 |
26 | 4.0 | 120.3125 | 91.0 | 4.429688 | 2.140625 | 16.703125 | 0.0 | 1.0 | 5.0 | 2.0 |
27 | 4.0 | 95.1250 | 113.0 | 3.769531 | 1.512695 | 16.906250 | 1.0 | 1.0 | 5.0 | 2.0 |
28 | 8.0 | 351.0000 | 264.0 | 4.218750 | 3.169922 | 14.500000 | 0.0 | 1.0 | 5.0 | 4.0 |
29 | 6.0 | 145.0000 | 175.0 | 3.619141 | 2.769531 | 15.500000 | 0.0 | 1.0 | 5.0 | 6.0 |
30 | 8.0 | 301.0000 | 335.0 | 3.539062 | 3.570312 | 14.601562 | 0.0 | 1.0 | 5.0 | 8.0 |
31 | 4.0 | 121.0000 | 109.0 | 4.109375 | 2.779297 | 18.593750 | 1.0 | 1.0 | 4.0 | 2.0 |
# Split the data into training and testing sets
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)
# Create a Bayesian RVFL regressor object
obj = lm.Regressor(method = "bayesianrvfl", nb_hidden = 5)
# Fit the model using the training data
obj.fit(X_train, y_train, reg_lambda=12.9155)
# Print the summary of the model
print(obj.summary(X_test, y=y_test, show_progress=False))
$R_squared
[1] 0.6416309
$R_squared_adj
[1] 1.537554
$Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.0724 -2.0122 -0.1018 -0.1941 1.4361 3.9676
$Coverage_rate
[1] 100
$citests
estimate lower upper p-value signif
cyl -24.5943583 -40.407994 -8.7807230 8.909365e-03 **
disp -0.2419797 -0.370835 -0.1131245 3.711077e-03 **
hp -1.5734483 -1.722903 -1.4239939 2.255640e-07 ***
drat 142.5646192 124.575179 160.5540599 1.217808e-06 ***
wt -144.8871352 -158.911143 -130.8631275 2.523441e-07 ***
qsec 46.8290859 27.829411 65.8287611 9.388045e-04 ***
vs 75.0555146 30.645127 119.4659017 6.110043e-03 **
am 207.5935234 133.205572 281.9814744 4.843095e-04 ***
gear 73.6892658 60.186232 87.1922995 1.091470e-05 ***
carb -71.2974988 -79.480400 -63.1145974 6.944475e-07 ***
$effects
── Data Summary ────────────────────────
Values
Name effects
Number of rows 7
Number of columns 10
_______________________
Column type frequency:
numeric 10
________________________
Group variables None
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable mean sd p0 p25 p50 p75 p100
1 cyl -24.6 17.1 -38.5 -38.5 -33.4 -12.4 1.66
2 disp -0.242 0.139 -0.351 -0.351 -0.285 -0.181 0.00546
3 hp -1.57 0.162 -1.90 -1.61 -1.48 -1.47 -1.47
4 drat 143. 19.5 125. 125. 141. 154. 174.
5 wt -145. 15.2 -167. -152. -142. -142. -117.
6 qsec 46.8 20.5 14.1 35.3 55.7 62.4 62.4
7 vs 75.1 48.0 37.2 37.2 58.7 93.9 167.
8 am 208. 80.4 64.3 168. 250. 267. 267.
9 gear 73.7 14.6 60.6 60.6 72.7 82.1 96.9
10 carb -71.3 8.85 -84.4 -75.2 -69.7 -69.7 -55.2
hist
1 ▇▁▂▁▃
2 ▇▂▁▂▂
3 ▂▁▁▃▇
4 ▇▂▂▂▂
5 ▂▅▇▁▂
6 ▃▁▁▂▇
7 ▇▁▃▁▂
8 ▂▂▁▂▇
9 ▇▂▂▂▂
10 ▂▅▇▁▂
# Select the first test sample
newx = X_test.iloc[0,:]
newy = y_test[0]
# Update the model with the new sample
new_X_test = X_test[1:]
new_y_test = y_test[1:]
obj.update(newx, newy, method="polyak", alpha=0.9)
# Print the summary of the model
print(obj.summary(new_X_test, y=new_y_test, show_progress=False))
$R_squared
[1] 0.6051442
$R_squared_adj
[1] 1.394856
$Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.6214 -2.5055 -1.5003 -0.8308 1.2738 3.2794
$Coverage_rate
[1] 100
$citests
estimate lower upper p-value signif
cyl -30.0502823 -48.4171958 -11.683369 8.442658e-03 **
disp -0.2958477 -0.4386085 -0.153087 3.121989e-03 **
hp -1.6053302 -1.6789750 -1.531685 3.424156e-08 ***
drat 153.7968829 131.5239191 176.069847 1.041460e-05 ***
wt -155.1954135 -174.4144275 -135.976399 4.804729e-06 ***
qsec 49.8967685 26.9993778 72.794159 2.504905e-03 **
vs 87.4170764 32.4599776 142.374175 9.457226e-03 **
am 214.5918910 119.8712855 309.312496 2.108825e-03 **
gear 83.1355825 65.3159018 100.955263 7.110354e-05 ***
carb -77.1384645 -88.2087477 -66.068181 9.958425e-06 ***
$effects
── Data Summary ────────────────────────
Values
Name effects
Number of rows 6
Number of columns 10
_______________________
Column type frequency:
numeric 10
________________________
Group variables None
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable mean sd p0 p25 p50 p75 p100
1 cyl -30.1 17.5 -42.5 -42.5 -39.9 -17.7 -4.35
2 disp -0.296 0.136 -0.377 -0.377 -0.343 -0.308 -0.0269
3 hp -1.61 0.0702 -1.70 -1.66 -1.57 -1.57 -1.53
4 drat 154. 21.2 137. 137. 144. 169. 185.
5 wt -155. 18.3 -182. -160. -154. -154. -125.
6 qsec 49.9 21.8 18.8 33.3 60.6 66.6 66.6
7 vs 87.4 52.4 46.7 46.7 73.7 105. 178.
8 am 215. 90.3 65.6 169. 266. 275. 275.
9 gear 83.1 17.0 70.0 70.0 75.5 94.1 109.
10 carb -77.1 10.5 -92.5 -79.9 -76.6 -76.6 -59.7
hist
1 ▇▁▁▁▃
2 ▇▂▁▁▂
3 ▅▁▁▇▂
4 ▇▂▁▁▅
5 ▂▂▇▁▂
6 ▅▁▁▂▇
7 ▇▁▅▁▂
8 ▂▂▁▁▇
9 ▇▂▁▂▂
10 ▂▂▇▁▂
Comments powered by Talkyard.