diff --git a/pyproject.toml b/pyproject.toml
index 7e415c9..3f71cdf 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,5 +4,7 @@ requires = [
     "setuptools>=42",
     "wheel",
     "pybind11>=2.6.0",
+    "scipy",
+    "sklearn",
     "numpy"
 ]
diff --git a/sortedl1/__init__.py b/sortedl1/__init__.py
index 6203bcb..386b4c4 100644
--- a/sortedl1/__init__.py
+++ b/sortedl1/__init__.py
@@ -1 +1 @@
-from .slope import slope
+from .estimators import Slope
diff --git a/sortedl1/estimators.py b/sortedl1/estimators.py
new file mode 100644
index 0000000..bb019e8
--- /dev/null
+++ b/sortedl1/estimators.py
@@ -0,0 +1,73 @@
+import _sortedl1 as sl1
+import numpy as np
+from sklearn.base import BaseEstimator, RegressorMixin
+from sklearn.utils.validation import check_array, check_is_fitted, check_X_y
+
+
+class Slope(BaseEstimator, RegressorMixin):
+    """
+    Sorted L-One Penalized Estimation
+
+    Parameters
+    ----------
+    lam : array_like
+        The lambda parameter vector for the Sorted L1 Penalty
+    alph : array_like
+        A multiplier for the Sorted L1 Penalty
+
+    Attributes
+    ----------
+    coef_ : ndarray
+        The estimated regression coefficients.
+    """
+
+    def __init__(self, lam=None, alph=1.0):
+        self.lam = lam
+        self.alph = alph
+
+    def fit(self, X, y):
+        """
+        Fit the model according to the given training data.
+
+        Parameters
+        ----------
+        X : array_like, shape (n_samples, n_features)
+            Training vector, where n_samples is the number of samples and
+            n_features is the number of features.
+        y : array_like, shape (n_samples,)
+            Target vector relative to X.
+
+        Returns
+        -------
+        self : object
+        """
+        X, y = check_X_y(X, y, accept_sparse=True, order="F", y_numeric=True)
+
+        self.X_ = X
+        self.y_ = y
+
+        self.coef_ = sl1.fit_slope(
+            self.X_, self.y_, np.array(self.lam), np.array(self.alph)
+        )
+
+        return self
+
+    def predict(self, X):
+        """
+        Predict using the linear model.
+
+        Parameters
+        ----------
+        X : array_like, shape (n_samples, n_features)
+            Samples.
+
+        Returns
+        -------
+        C : array, shape (n_samples,)
+            Returns predicted values.
+        """
+        check_is_fitted(self)
+
+        X = check_array(X, accept_sparse=True, order="F")
+
+        return np.dot(X, self.coef_)
diff --git a/sortedl1/slope.py b/sortedl1/slope.py
deleted file mode 100644
index e9c0f21..0000000
--- a/sortedl1/slope.py
+++ /dev/null
@@ -1,34 +0,0 @@
-import _sortedl1 as sl1
-import numpy as np
-
-
-def slope(x, y, lam, alph):
-    """
-    Sorted L-One Penalized Estimation
-
-    Parameters
-    ----------
-    x : array_like
-        The design matrix
-    y : array_like
-        The response vector
-    lam : array_like
-        The lambda paramter vector for the Sorted L1 Penalty
-    alph : array_like
-        A multiplied for the Sorted L1 Penalty
-
-    Returns
-    -------
-    res : ndarray
-        An object containing the estimated regression coefficients.
-
-    """
-    x = np.array(x)
-    y = np.array(y)
-
-    lam = np.array(lam)
-    alpha = np.array(alph)
-
-    res = sl1.fit_slope(x, y, lam, alpha)
-
-    return res