diff --git a/R/checks.R b/R/checks.R
index 6bf7e2a..2cc9cdb 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -120,6 +120,13 @@ check.centers <- function(nu, centers, ncenters) {
   centers
 }
 
+check.intercept <- function(intercept) {
+  if (!is.logical(intercept)) {
+    stop("'intercept' must be either 'TRUE' or 'FALSE'")
+  }
+  intercept
+}
+
 check.symmetric <- function(nu, centers) {
   if(isTRUE(all.equal(nu, centers))) {
     symmetric <- TRUE
diff --git a/R/predict.R b/R/predict.R
index 8f6ff1f..0532804 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -23,13 +23,18 @@ predict.ulsif <- function(object, newdata = NULL, sigma = c("sigmaopt", "all"),
   newlambda <- check.lambda.predict(object, lambda)
   newdata <- check.newdata(object, newdata)
 
-  alpha   <- extract.alpha(object, newsigma, newlambda)
-  nsigma  <- dim(alpha)[2]
-  nlambda <- dim(alpha)[3]
-  dratio  <- array(0, c(nrow(newdata), nlambda, nsigma))
+  alpha     <- extract.alpha(object, newsigma, newlambda)
+  nsigma    <- length(newsigma)
+  nlambda   <- length(newlambda)
+  dratio    <- array(0, c(nrow(newdata), nlambda, nsigma))
+  intercept <- nrow(object$alpha) > nrow(object$centers)
 
   for (i in 1:nsigma) {
-    K <- distance(newdata, object$centers) |> kernel_gaussian(newsigma[i])
+    if (intercept) {
+      K <- cbind(0, distance(newdata, object$centers)) |> kernel_gaussian(newsigma[i])
+    } else {
+      K <- distance(newdata, object$centers) |> kernel_gaussian(newsigma[i])
+    }
     for (j in 1:nlambda) {
       dratio[ , i, j] <- K %*% alpha[, i, j]
     }
diff --git a/R/ulsif.R b/R/ulsif.R
index 2737e58..be69c57 100644
--- a/R/ulsif.R
+++ b/R/ulsif.R
@@ -5,6 +5,8 @@
 #' @param df_denominator \code{data.frame} with exclusively numeric variables
 #' with the denominator samples (must have the same variables as
 #' \code{df_denominator})
+#' @param intercept \code{logical} Indicating whether to include an intercept
+#' term in the model. Defaults to \code{TRUE}.
 #' @param nsigma Integer indicating the number of sigma values (bandwidth
 #' parameter of the Gaussian kernel gram matrix) to use in cross-validation.
 #' @param sigma_quantile \code{NULL} or numeric vector with probabilities to
@@ -41,10 +43,10 @@
 #' ulsif(x, y)
 #' ulsif(x, y, sigma = 2, lambda = 2)
 
-ulsif <- function(df_numerator, df_denominator, nsigma = 10, sigma_quantile = NULL,
-                  sigma = NULL, nlambda = 20, lambda = NULL, ncenters = 200,
-                  centers = NULL, parallel = FALSE, nthreads = NULL,
-                  progressbar = TRUE) {
+ulsif <- function(df_numerator, df_denominator, intercept = TRUE, nsigma = 10,
+                  sigma_quantile = NULL, sigma = NULL, nlambda = 20,
+                  lambda = NULL, ncenters = 200, centers = NULL,
+                  parallel = FALSE, nthreads = NULL, progressbar = TRUE) {
 
   cl <- match.call()
   nu <- as.matrix(df_numerator)
@@ -55,9 +57,14 @@ ulsif <- function(df_numerator, df_denominator, nsigma = 10, sigma_quantile = NU
   symmetric <- check.symmetric(nu, centers)
   parallel  <- check.parallel(parallel, nthreads, sigma, lambda)
   nthreads  <- check.threads(parallel, nthreads)
+  intercept <- check.intercept(intercept)
 
   dist_nu <- distance(nu, centers, symmetric)
   dist_de <- distance(de, centers)
+  if (intercept) {
+    dist_nu <- cbind(0, dist_nu)
+    dist_de <- cbind(0, dist_de)
+  }
 
   sigma  <- check.sigma(nsigma, sigma_quantile, sigma, dist_nu)
   lambda <- check.lambda(nlambda, lambda)
diff --git a/man/ulsif.Rd b/man/ulsif.Rd
index b880ee4..08a5576 100644
--- a/man/ulsif.Rd
+++ b/man/ulsif.Rd
@@ -7,6 +7,7 @@
 ulsif(
   df_numerator,
   df_denominator,
+  intercept = TRUE,
   nsigma = 10,
   sigma_quantile = NULL,
   sigma = NULL,
@@ -27,6 +28,9 @@ the numerator samples}
 with the denominator samples (must have the same variables as
 \code{df_denominator})}
 
+\item{intercept}{\code{logical} Indicating whether to include an intercept
+term in the model. Defaults to \code{TRUE}.}
+
 \item{nsigma}{Integer indicating the number of sigma values (bandwidth
 parameter of the Gaussian kernel gram matrix) to use in cross-validation.}