## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----message=FALSE, warning=FALSE---------------------------------------------
# Load the package
require(AddiVortes)

# Separate predictors (X) and the response variable (Y)
X_Boston <- as.matrix(Boston[, 1:13])
Y_Boston <- as.numeric(as.matrix(Boston[, 14]))

## -----------------------------------------------------------------------------
n <- length(Y_Boston)

# Set a seed for reproducibility
set.seed(1025)

# Create a training set containing 5/6 of the data
TrainSet <- sort(sample.int(n, 5 * n / 6))

# The remaining data will be our test set
TestSet <- setdiff(1:n, TrainSet)

## -----------------------------------------------------------------------------
# Run the AddiVortes algorithm on the training data
results <- AddiVortes(
  y = Y_Boston[TrainSet],
  x = X_Boston[TrainSet, ],
  m = 200,
  totalMCMCIter = 2000,
  mcmcBurnIn = 500,
  nu = 6,
  q = 0.85,
  k = 3,
  sd = 0.8,
  Omega = 3,
  LambdaRate = 25,
  InitialSigma = "Linear",
  showProgress = FALSE
)

## -----------------------------------------------------------------------------
# Generate predictions on the test set
preds <- predict(results,
  X_Boston[TestSet, ],
  showProgress = FALSE
)

# The RMSE is contained in the results object
cat("In-Sample RMSE:", results$inSampleRmse, "\n")

# Calculate the Root Mean Squared Error (RMSE) for the test set
rmse <- sqrt(mean((Y_Boston[TestSet] - preds)^2))
cat("Test Set RMSE:", rmse, "\n")

## ----fig.width=6, fig.height=6, fig.align='center'----------------------------
# Plot true values vs. predicted values
plot(Y_Boston[TestSet],
  preds,
  xlab = "True Values",
  ylab = "Predicted Values",
  main = "AddiVortes Predictions vs True Values",
  xlim = range(c(Y_Boston[TestSet], preds)),
  ylim = range(c(Y_Boston[TestSet], preds)),
  pch = 19, col = "darkblue"
)

# Add the line of equality (y = x) for reference
abline(a = 0, b = 1, col = "darkred", lwd = 2)

# Get quantile predictions to create error bars/intervals
preds_quantile <- predict(results,
  X_Boston[TestSet, ],
  "quantile",
  showProgress = FALSE
)

# Add error segments for each prediction
for (i in seq_len(nrow(preds_quantile))) {
  segments(Y_Boston[TestSet][i], preds_quantile[i, 1],
    Y_Boston[TestSet][i], preds_quantile[i, 2],
    col = "darkblue", lwd = 1
  )
}
legend("topleft",
  legend = c("Prediction", "y=x Line", "95% Interval"),
  col = c("darkblue", "darkred", "darkblue"),
  lty = 1, pch = c(19, NA, NA), lwd = 2
)