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)

dtype: object

y.dtype

dtype('float64')

X

# 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 ▂▂▇▁▂