Example 1: MPG Prediction (mtcars dataset)
Load and prepare data
data(mtcars)
set.seed(1243)
train_idx <- sample(nrow(mtcars), size = floor(0.8 * nrow(mtcars)))
train_data <- mtcars[train_idx, ]
test_data <- mtcars[-train_idx, -1]
Fit models
## elapsed
## 0.017
##
## Call:
## lm(formula = mpg ~ ., data = train_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5211 -0.9792 -0.0324 1.1808 4.9814
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.054416 25.456900 -0.199 0.8455
## cyl 0.695392 1.396506 0.498 0.6262
## disp 0.005254 0.017342 0.303 0.7664
## hp -0.007610 0.027723 -0.274 0.7877
## drat 4.128157 2.724353 1.515 0.1520
## wt -1.621396 2.139071 -0.758 0.4610
## qsec 0.064356 0.932144 0.069 0.9459
## vs 0.138716 3.421183 0.041 0.9682
## am -0.498476 2.956568 -0.169 0.8685
## gear 4.402648 2.287816 1.924 0.0749 .
## carb -1.999389 1.299580 -1.538 0.1462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.464 on 14 degrees of freedom
## Multiple R-squared: 0.8938, Adjusted R-squared: 0.818
## F-statistic: 11.79 on 10 and 14 DF, p-value: 3.4e-05
## 2.5 % 97.5 %
## (Intercept) -59.65403559 49.54520296
## cyl -2.29981561 3.69060001
## disp -0.03194096 0.04244882
## hp -0.06707095 0.05185084
## drat -1.71500030 9.97131342
## wt -6.20924769 2.96645550
## qsec -1.93489537 2.06360651
## vs -7.19899241 7.47642359
## am -6.83968216 5.84273112
## gear -0.50422869 9.30952400
## carb -4.78671119 0.78793282
## elapsed
## 0.193
##
## Call:
## engine(formula = y ~ . - 1, data = df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4586 -1.7751 -0.4056 1.5360 5.2494
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## cyl -0.2934 0.6109 -0.480 0.638
## disp -0.4021 0.6033 -0.667 0.516
## hp -0.4266 0.6034 -0.707 0.491
## drat 0.5285 0.5736 0.921 0.372
## wt -0.5200 0.5947 -0.874 0.397
## qsec 0.1419 0.5949 0.238 0.815
## vs 0.2738 0.5787 0.473 0.643
## am 0.4077 0.5711 0.714 0.487
## gear 0.2922 0.5785 0.505 0.621
## carb -0.3778 0.5421 -0.697 0.497
## h1 -0.2444 0.6039 -0.405 0.692
## h2 0.6964 0.5848 1.191 0.253
## h3 -0.3387 0.5802 -0.584 0.569
## h4 -0.1628 0.5448 -0.299 0.769
## h5 0.5906 0.5822 1.014 0.328
##
## Residual standard error: 2.925 on 14 degrees of freedom
## Multiple R-squared: 0.7534, Adjusted R-squared: 0.4892
## F-statistic: 2.851 on 15 and 14 DF, p-value: 0.02863
## 2.5 % 97.5 %
## cyl -1.6037724 1.0168960
## disp -1.6959495 0.8917476
## hp -1.7208611 0.8676077
## drat -0.7018175 1.7587974
## wt -1.7955969 0.7555437
## qsec -1.1340563 1.4177867
## vs -0.9673542 1.5149112
## am -0.8172659 1.6326674
## gear -0.9485594 1.5330304
## carb -1.5405137 0.7848791
## h1 -1.5396642 1.0508179
## h2 -0.5577916 1.9506325
## h3 -1.5831708 0.9057347
## h4 -1.3312796 1.0055891
## h5 -0.6580768 1.8392040
#print(simulate(ridge_model, newdata = test_data))
Make predictions
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")
Compare predictions
results <- data.frame(
Actual = mtcars[-train_idx, ]$mpg,
LM_Pred = lm_pred[,"fit"],
LM_Lower = lm_pred[,"lwr"],
LM_Upper = lm_pred[,"upr"],
Ridge_Pred = ridge_pred[,"fit"],
Ridge_Lower = ridge_pred[,"lwr"],
Ridge_Upper = ridge_pred[,"upr"]
)
# Print results
print("Prediction Intervals Comparison:")
## [1] "Prediction Intervals Comparison:"
## Actual LM_Pred LM_Lower LM_Upper Ridge_Pred Ridge_Lower
## Valiant 18.1 17.93324 10.149847 25.71663 18.96396 15.81971
## Merc 280C 17.8 20.63530 13.636618 27.63398 19.75734 15.72996
## Toyota Corolla 33.9 28.58373 22.379666 34.78779 26.85749 23.43556
## Camaro Z28 13.3 15.85710 8.140858 23.57335 15.29984 11.44900
## Porsche 914-2 26.0 31.07535 18.988702 43.16201 25.59502 22.15517
## Ford Pantera L 15.8 27.07516 14.930150 39.22016 19.74423 16.25442
## Ridge_Upper
## Valiant 22.78275
## Merc 280C 23.97086
## Toyota Corolla 32.21672
## Camaro Z28 18.58237
## Porsche 914-2 28.98960
## Ford Pantera L 23.25351
# Calculate coverage and Winkler scores
lm_coverage <- mean(mtcars[-train_idx, ]$mpg >= results$LM_Lower &
mtcars[-train_idx, ]$mpg <= results$LM_Upper)
ridge_coverage <- mean(mtcars[-train_idx, ]$mpg >= results$Ridge_Lower &
mtcars[-train_idx, ]$mpg <= results$Ridge_Upper)
lm_winkler <- misc::winkler_score(mtcars[-train_idx, ]$mpg, results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(mtcars[-train_idx, ]$mpg, results$Ridge_Lower, results$Ridge_Upper)
print(sprintf("\nPrediction interval metrics:"))
## [1] "\nPrediction interval metrics:"
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score",
100 * lm_coverage, mean(lm_winkler)))
## [1] "Linear Model: 100.0% coverage, 18.226 Winkler score"
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score",
100 * ridge_coverage, mean(ridge_winkler)))
## [1] "Calibrated Model: 71.4% coverage, 19.857 Winkler score"
# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
results$Ridge_Lower, results$Ridge_Upper))
# Plot prediction intervals
par(mfrow=c(1,2))
# Linear Model Plot
plot(results$Actual, results$LM_Pred,
main="Linear Model Predictions",
xlab="Actual MPG", ylab="Predicted MPG",
ylim=y_limits)
# Add shaded prediction intervals
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred) # Replot points over shading
abline(0, 1, col="red", lty=2) # Add diagonal line
# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
main="Ridge Model Predictions",
xlab="Actual MPG", ylab="Predicted MPG",
ylim=y_limits)
# Add shaded prediction intervals
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred) # Replot points over shading
abline(0, 1, col="red", lty=2) # Add diagonal line

# Add simulation plot
par(mfrow=c(1,1))
# Generate 100 simulations
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
# Plot simulations
matplot(sims, type = "l",
col = rgb(0, 0, 1, 0.1), lty = 1,
xlab = "obs. #", ylab = "Simulated MPG",
main = "Ridge Model Simulations")
lines(mtcars[-train_idx, ]$mpg, col = "red")
Example 2: Boston Housing Price Prediction
Load and prepare data
library(MASS)
data(Boston)
set.seed(1243)
train_idx <- sample(nrow(Boston), size = floor(0.8 * nrow(Boston)))
train_data <- Boston[train_idx, ]
test_data <- Boston[-train_idx, -14] # -14 removes 'medv' (target variable)
Fit models
# Fit regular linear model
start <- proc.time()[3]
lm_model <- rvfl::rvfl(x = as.matrix(train_data[,-14]), y = train_data$medv, activ = "linear")#lm(medv ~ ., data = train_data)
print(proc.time()[3] - start)
## elapsed
## 0.097
##
## Call:
## engine(formula = y ~ . - 1, data = df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8740 -2.8709 -0.5828 1.8772 26.0796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## crim -0.97711 1.03548 -0.944 0.3465
## zn 1.05329 1.00842 1.044 0.2975
## indus 0.25561 1.09993 0.232 0.8165
## chas 0.69443 0.85692 0.810 0.4187
## nox -1.65544 1.11978 -1.478 0.1409
## rm 2.54157 1.01833 2.496 0.0134 *
## age 0.59156 1.15721 0.511 0.6098
## dis -2.50321 1.19731 -2.091 0.0378 *
## rad 2.05723 1.09839 1.873 0.0625 .
## tax -1.12820 1.30239 -0.866 0.3874
## ptratio -1.49972 1.00988 -1.485 0.1391
## black 0.88341 0.90808 0.973 0.3318
## lstat -3.85354 0.94186 -4.091 6.18e-05 ***
## h1 -0.12471 1.90545 -0.065 0.9479
## h2 0.12471 1.90545 0.065 0.9479
## h3 0.12541 1.60619 0.078 0.9378
## h4 -0.06691 1.61069 -0.042 0.9669
## h5 0.68743 1.65118 0.416 0.6776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.062 on 203 degrees of freedom
## Multiple R-squared: 0.69, Adjusted R-squared: 0.6625
## F-statistic: 25.1 on 18 and 203 DF, p-value: < 2.2e-16
## 2.5 % 97.5 %
## crim -3.0187805 1.0645547
## zn -0.9350311 3.0416020
## indus -1.9131387 2.4243638
## chas -0.9951706 2.3840263
## nox -3.8633281 0.5524414
## rm 0.5337147 4.5494251
## age -1.6901406 2.8732508
## dis -4.8639566 -0.1424538
## rad -0.1084805 4.2229500
## tax -3.6961423 1.4397357
## ptratio -3.4909239 0.4914810
## black -0.9070666 2.6738816
## lstat -5.7106232 -1.9964644
## h1 -3.8817198 3.6323082
## h2 -3.6323082 3.8817198
## h3 -3.0415556 3.2923664
## h4 -3.2427373 3.1089125
## h5 -2.5682466 3.9430974
## elapsed
## 0.135
##
## Call:
## engine(formula = y ~ . - 1, data = df_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6436 -2.5860 -0.4639 2.0094 24.4992
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## crim -3.0623 1.1491 -2.665 0.008317 **
## zn 0.2952 0.9853 0.300 0.764763
## indus 0.3479 0.9903 0.351 0.725768
## chas 0.7550 0.7475 1.010 0.313689
## nox -2.4144 0.9770 -2.471 0.014291 *
## rm 1.6957 0.8527 1.989 0.048091 *
## age 1.0275 0.9888 1.039 0.299986
## dis -3.8911 1.0732 -3.626 0.000364 ***
## rad 2.9026 1.1001 2.639 0.008970 **
## tax -1.8431 1.2394 -1.487 0.138527
## ptratio -0.7089 0.9333 -0.760 0.448393
## black 1.3562 1.0188 1.331 0.184607
## lstat -3.6097 0.8140 -4.434 1.51e-05 ***
## h1 1.4467 1.5982 0.905 0.366441
## h2 2.2793 1.6075 1.418 0.157757
## h3 2.0981 0.6087 3.447 0.000689 ***
## h4 1.8865 1.3423 1.405 0.161431
## h5 1.8774 0.7984 2.352 0.019650 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.403 on 203 degrees of freedom
## Multiple R-squared: 0.7655, Adjusted R-squared: 0.7447
## F-statistic: 36.81 on 18 and 203 DF, p-value: < 2.2e-16
## 2.5 % 97.5 %
## crim -5.32791535 -0.7966590
## zn -1.64742341 2.2378519
## indus -1.60482173 2.3005293
## chas -0.71890031 2.2289375
## nox -4.34089360 -0.4879895
## rm 0.01437801 3.3770093
## age -0.92219782 2.9772528
## dis -6.00709823 -1.7751782
## rad 0.73360783 5.0715876
## tax -4.28680612 0.6005612
## ptratio -2.54904493 1.1312701
## black -0.65250165 3.3648953
## lstat -5.21467118 -2.0046819
## h1 -1.70458059 4.5979892
## h2 -0.89028309 5.4487882
## h3 0.89797558 3.2982400
## h4 -0.76016904 4.5331138
## h5 0.30329787 3.4515435
#print(simulate(ridge_model, newdata = test_data))
Make predictions and compare
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")
results <- data.frame(
Actual = Boston[-train_idx, ]$medv,
LM_Pred = lm_pred[,"fit"],
LM_Lower = lm_pred[,"lwr"],
LM_Upper = lm_pred[,"upr"],
Ridge_Pred = ridge_pred[,"fit"],
Ridge_Lower = ridge_pred[,"lwr"],
Ridge_Upper = ridge_pred[,"upr"]
)
# Print results
print("Prediction Intervals Comparison:")
## [1] "Prediction Intervals Comparison:"
## Actual LM_Pred LM_Lower LM_Upper Ridge_Pred Ridge_Lower Ridge_Upper
## 1 24.0 29.79046 19.608073 39.97285 28.75770 19.568995 37.96624
## 3 34.7 30.96586 20.872948 41.05878 34.14903 25.222195 40.80861
## 4 33.4 29.08560 18.927692 39.24351 34.41607 25.955182 42.07926
## 18 17.5 17.88680 7.769896 28.00369 14.73369 7.990620 25.30735
## 21 13.6 13.55881 3.342978 23.77464 13.65759 6.005290 22.40761
## 24 14.5 15.04997 4.841070 25.25886 13.86386 5.700152 22.48673
# Calculate coverage and Winkler scores
lm_coverage <- mean(Boston[-train_idx, ]$medv >= results$LM_Lower &
Boston[-train_idx, ]$medv <= results$LM_Upper)
ridge_coverage <- mean(Boston[-train_idx, ]$medv >= results$Ridge_Lower &
Boston[-train_idx, ]$medv <= results$Ridge_Upper)
lm_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv, results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv, results$Ridge_Lower, results$Ridge_Upper)
print(sprintf("\nPrediction interval metrics:"))
## [1] "\nPrediction interval metrics:"
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score",
100 * lm_coverage, mean(lm_winkler)))
## [1] "Linear Model: 94.1% coverage, 26.946 Winkler score"
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score",
100 * ridge_coverage, mean(ridge_winkler)))
## [1] "Calibrated Model: 95.1% coverage, 25.419 Winkler score"
# Visualization
# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
results$Ridge_Lower, results$Ridge_Upper))
par(mfrow=c(1,2))
# Linear Model Plot
plot(results$Actual, results$LM_Pred,
main="Linear Model Predictions",
xlab="Actual Median Value", ylab="Predicted Median Value",
ylim=y_limits)
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred)
abline(0, 1, col="red", lty=2)
# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
main="Ridge Model Predictions",
xlab="Actual Median Value", ylab="Predicted Median Value",
ylim=y_limits)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred)
abline(0, 1, col="red", lty=2)

# Add simulation plot
par(mfrow=c(1,1))
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
matplot(sims, type = "l",
col = rgb(0, 0, 1, 0.1), lty = 1,
xlab = "obs. #", ylab = "Simulated Median Value",
main = "Ridge Model Simulations")
lines(Boston[-train_idx, ]$medv, col = "red")
Compare different activation functions
# Adjust margins to fit plots
par(mfrow=c(2,2), mar=c(4,4,2,1)) # Bottom, left, top, right margins
# Fit models with different activation functions
ridge_relu <- rvfl::rvfl(x = as.matrix(train_data[,-14]), y = train_data$medv, activ = "relu")
ridge_sigmoid <- rvfl::rvfl(x = as.matrix(train_data[,-14]), y = train_data$medv, activ = "sigmoid")
ridge_tanh <- rvfl::rvfl(x = as.matrix(train_data[,-14]), y = train_data$medv, activ = "tanh")
# Make predictions
relu_pred <- predict(ridge_relu, newdata = as.matrix(test_data), method = "gaussian")
sigmoid_pred <- predict(ridge_sigmoid, newdata = as.matrix(test_data), method = "gaussian")
tanh_pred <- predict(ridge_tanh, newdata = as.matrix(test_data), method = "gaussian")
# Combine results
results_all <- data.frame(
Actual = Boston[-train_idx, ]$medv,
LM_Pred = lm_pred[,"fit"],
LM_Lower = lm_pred[,"lwr"],
LM_Upper = lm_pred[,"upr"],
ReLU_Pred = relu_pred[,"fit"],
ReLU_Lower = relu_pred[,"lwr"],
ReLU_Upper = relu_pred[,"upr"],
Sigmoid_Pred = sigmoid_pred[,"fit"],
Sigmoid_Lower = sigmoid_pred[,"lwr"],
Sigmoid_Upper = sigmoid_pred[,"upr"],
Tanh_Pred = tanh_pred[,"fit"],
Tanh_Lower = tanh_pred[,"lwr"],
Tanh_Upper = tanh_pred[,"upr"]
)
# Calculate coverage and Winkler scores for each model
lm_coverage <- mean(Boston[-train_idx, ]$medv >= results_all$LM_Lower &
Boston[-train_idx, ]$medv <= results_all$LM_Upper)
relu_coverage <- mean(Boston[-train_idx, ]$medv >= results_all$ReLU_Lower &
Boston[-train_idx, ]$medv <= results_all$ReLU_Upper)
sigmoid_coverage <- mean(Boston[-train_idx, ]$medv >= results_all$Sigmoid_Lower &
Boston[-train_idx, ]$medv <= results_all$Sigmoid_Upper)
tanh_coverage <- mean(Boston[-train_idx, ]$medv >= results_all$Tanh_Lower &
Boston[-train_idx, ]$medv <= results_all$Tanh_Upper)
lm_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv,
results_all$LM_Lower, results_all$LM_Upper)
relu_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv,
results_all$ReLU_Lower, results_all$ReLU_Upper)
sigmoid_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv,
results_all$Sigmoid_Lower, results_all$Sigmoid_Upper)
tanh_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv,
results_all$Tanh_Lower, results_all$Tanh_Upper)
print(sprintf("\nPrediction interval metrics:"))
## [1] "\nPrediction interval metrics:"
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score",
100 * lm_coverage, mean(lm_winkler)))
## [1] "Linear Model: 94.1% coverage, 26.946 Winkler score"
print(sprintf("ReLU Model: %.1f%% coverage, %.3f Winkler score",
100 * relu_coverage, mean(relu_winkler)))
## [1] "ReLU Model: 95.1% coverage, 25.419 Winkler score"
print(sprintf("Sigmoid Model: %.1f%% coverage, %.3f Winkler score",
100 * sigmoid_coverage, mean(sigmoid_winkler)))
## [1] "Sigmoid Model: 93.1% coverage, 31.566 Winkler score"
print(sprintf("Tanh Model: %.1f%% coverage, %.3f Winkler score",
100 * tanh_coverage, mean(tanh_winkler)))
## [1] "Tanh Model: 93.1% coverage, 30.259 Winkler score"
# Visualization
par(mfrow=c(2,2))
# Linear Model Plot
plot(results_all$Actual, results_all$LM_Pred,
main="Linear Model",
xlab="Actual", ylab="Predicted",
ylim=range(results_all[,c("LM_Lower","LM_Upper")]))
x_ordered <- order(results_all$Actual)
polygon(c(results_all$Actual[x_ordered], rev(results_all$Actual[x_ordered])),
c(results_all$LM_Lower[x_ordered], rev(results_all$LM_Upper[x_ordered])),
col=rgb(0,0,1,0.2), border=NA)
points(results_all$Actual, results_all$LM_Pred)
abline(0,1, col="red", lty=2)
# ReLU Model Plot
plot(results_all$Actual, results_all$ReLU_Pred,
main="ReLU Model",
xlab="Actual", ylab="Predicted",
ylim=range(results_all[,c("ReLU_Lower","ReLU_Upper")]))
polygon(c(results_all$Actual[x_ordered], rev(results_all$Actual[x_ordered])),
c(results_all$ReLU_Lower[x_ordered], rev(results_all$ReLU_Upper[x_ordered])),
col=rgb(0,0,1,0.2), border=NA)
points(results_all$Actual, results_all$ReLU_Pred)
abline(0,1, col="red", lty=2)
# Sigmoid Model Plot
plot(results_all$Actual, results_all$Sigmoid_Pred,
main="Sigmoid Model",
xlab="Actual", ylab="Predicted",
ylim=range(results_all[,c("Sigmoid_Lower","Sigmoid_Upper")]))
polygon(c(results_all$Actual[x_ordered], rev(results_all$Actual[x_ordered])),
c(results_all$Sigmoid_Lower[x_ordered], rev(results_all$Sigmoid_Upper[x_ordered])),
col=rgb(0,0,1,0.2), border=NA)
points(results_all$Actual, results_all$Sigmoid_Pred)
abline(0,1, col="red", lty=2)
# Tanh Model Plot
plot(results_all$Actual, results_all$Tanh_Pred,
main="Tanh Model",
xlab="Actual", ylab="Predicted",
ylim=range(results_all[,c("Tanh_Lower","Tanh_Upper")]))
polygon(c(results_all$Actual[x_ordered], rev(results_all$Actual[x_ordered])),
c(results_all$Tanh_Lower[x_ordered], rev(results_all$Tanh_Upper[x_ordered])),
col=rgb(0,0,1,0.2), border=NA)
points(results_all$Actual, results_all$Tanh_Pred)
abline(0,1, col="red", lty=2)

Simulation plots
# Adjust margins to fit plots
par(mfrow=c(2,2), mar=c(4,4,2,1)) # Bottom, left, top, right margins
# Linear Model simulations
sims_lm <- simulate(lm_model, newdata = test_data, nsim = 500)
matplot(sims_lm, type = "l",
col = rgb(0,0,1,0.1), lty = 1,
xlab = "obs. #", ylab = "Simulated Value",
main = "Linear Model Simulations")
lines(Boston[-train_idx, ]$medv, col = "red")
# ReLU Model simulations
sims_relu <- simulate(ridge_relu, newdata = as.matrix(test_data), nsim = 500)
matplot(sims_relu, type = "l",
col = rgb(0,0,1,0.1), lty = 1,
xlab = "obs. #", ylab = "Simulated Value",
main = "ReLU Model Simulations")
lines(Boston[-train_idx, ]$medv, col = "red")
# Sigmoid Model simulations
sims_sigmoid <- simulate(ridge_sigmoid, newdata = as.matrix(test_data), nsim = 500)
matplot(sims_sigmoid, type = "l",
col = rgb(0,0,1,0.1), lty = 1,
xlab = "obs. #", ylab = "Simulated Value",
main = "Sigmoid Model Simulations")
lines(Boston[-train_idx, ]$medv, col = "red")
# Tanh Model simulations
sims_tanh <- simulate(ridge_tanh, newdata = as.matrix(test_data), nsim = 500)
matplot(sims_tanh, type = "l",
col = rgb(0,0,1,0.1), lty = 1,
xlab = "obs. #", ylab = "Simulated Value",
main = "Tanh Model Simulations")
lines(Boston[-train_idx, ]$medv, col = "red")
