diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 6c166db..6f2ce68 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -5,9 +5,8 @@ name: tests
 
 on:
   push:
-    branches: [ "main", "dev" ]
+    branches: [ "**" ]
   pull_request:
-    branches: [ "main", "dev" ]
 
 jobs:
   test:
@@ -16,7 +15,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12", "3.13", "3.14"]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13", "3.14"]
 
     steps:
     - uses: actions/checkout@v4
@@ -33,7 +32,7 @@ jobs:
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --exit-zero --select=E9,F63,F7,F82 --show-source --statistics
+        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
diff --git a/.gitignore b/.gitignore
index 8455419..c4898a8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -172,3 +172,5 @@ cython_debug/
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
 #.idea/
 
+# Pane worktrees (local only)
+worktrees/
diff --git a/PyNomaly/_compat.py b/PyNomaly/_compat.py
new file mode 100644
index 0000000..d6c11a5
--- /dev/null
+++ b/PyNomaly/_compat.py
@@ -0,0 +1,198 @@
+# Authors: Valentino Constantinou <vc@valentino.io>
+# License: Apache 2.0
+
+"""
+Optional scikit-learn integration.
+
+When scikit-learn is installed, LocalOutlierProbability inherits from
+``BaseEstimator`` and ``OutlierMixin`` for full ecosystem compatibility
+(pipelines, grid search, ``check_estimator``).
+
+When it is not installed, lightweight fallbacks provide the same public
+surface (``get_params``, ``set_params``, ``fit_predict``, ``check_is_fitted``)
+so PyNomaly works as a standalone package.
+
+Install the optional dependency with::
+
+    pip install PyNomaly[sklearn]
+"""
+
+from __future__ import annotations
+
+import inspect
+
+try:
+    from sklearn.base import BaseEstimator, OutlierMixin
+    from sklearn.exceptions import NotFittedError
+    from sklearn.utils.validation import check_array, check_is_fitted
+
+    _SKLEARN_AVAILABLE = True
+except ImportError:  # pragma: no cover - exercised only without sklearn
+    _SKLEARN_AVAILABLE = False
+
+    class NotFittedError(ValueError, AttributeError):
+        """Raised when an estimator method is used before fit."""
+
+    class BaseEstimator:
+        """Minimal BaseEstimator providing get_params / set_params."""
+
+        def get_params(self, deep=True):
+            params = {}
+            for name, param in inspect.signature(self.__init__).parameters.items():
+                if name == "self":
+                    continue
+                params[name] = getattr(self, name, param.default)
+            return params
+
+        def set_params(self, **params):
+            if not params:
+                return self
+            valid = self.get_params(deep=True)
+            for key, value in params.items():
+                if key not in valid:
+                    raise ValueError(
+                        f"Invalid parameter {key!r} for estimator "
+                        f"{type(self).__name__}. Valid parameters are: "
+                        f"{sorted(valid)!r}."
+                    )
+                setattr(self, key, value)
+            return self
+
+        def __repr__(self):
+            params = self.get_params(deep=False)
+            filtered = {
+                k: v
+                for k, v in params.items()
+                if v is not inspect.Parameter.empty
+            }
+            inner = ", ".join(f"{k}={v!r}" for k, v in filtered.items())
+            return f"{type(self).__name__}({inner})"
+
+    class OutlierMixin:
+        """Minimal OutlierMixin providing fit_predict."""
+
+        def fit_predict(self, X, y=None, **fit_params):
+            return self.fit(X, y, **fit_params).predict(X)
+
+    def check_is_fitted(estimator, attributes=None, *, msg=None):
+        if attributes is None:
+            fitted = [
+                a
+                for a in vars(estimator)
+                if a.endswith("_") and not a.startswith("__")
+            ]
+            if not fitted:
+                raise NotFittedError(
+                    msg
+                    or (
+                        f"This {type(estimator).__name__} instance is not "
+                        f"fitted yet. Call 'fit' with appropriate arguments "
+                        f"before using this estimator."
+                    )
+                )
+            return
+
+        if isinstance(attributes, str):
+            attributes = [attributes]
+        missing = [a for a in attributes if not hasattr(estimator, a)]
+        if missing:
+            raise NotFittedError(
+                msg
+                or (
+                    f"This {type(estimator).__name__} instance is not fitted "
+                    f"yet. Call 'fit' with appropriate arguments before using "
+                    f"this estimator."
+                )
+            )
+
+    def check_array(
+        array,
+        *,
+        accept_sparse=False,
+        dtype="numeric",
+        ensure_2d=True,
+        ensure_min_samples=1,
+        ensure_min_features=1,
+        copy=False,
+        **kwargs,
+    ):
+        """Fallback input validator approximating sklearn.utils.check_array."""
+        import numpy as np
+
+        try:
+            from scipy.sparse import issparse
+        except ImportError:
+
+            def issparse(x):
+                return False
+
+        if issparse(array):
+            if accept_sparse is False or accept_sparse == []:
+                raise TypeError(
+                    "A sparse matrix was passed, but dense data is required. "
+                    "Use X.toarray() to convert to a dense numpy array."
+                )
+            array = array.toarray()
+
+        if hasattr(array, "values") and hasattr(array, "columns"):
+            array = array.values
+
+        arr = np.asarray(array)
+        if copy:
+            arr = np.array(arr, copy=True)
+
+        if ensure_2d and arr.ndim == 1:
+            raise ValueError(
+                "Expected 2D array, got 1D array instead:\n"
+                f"array={arr}.\n"
+                "Reshape your data either using array.reshape(-1, 1) if "
+                "your data has a single feature or array.reshape(1, -1) if "
+                "it contains a single sample."
+            )
+        if ensure_2d and arr.ndim != 2:
+            raise ValueError(
+                f"Expected 2D array, got array with ndim={arr.ndim} instead."
+            )
+
+        if arr.size == 0:
+            n_samples = arr.shape[0] if arr.ndim >= 1 else 0
+            n_features = arr.shape[1] if arr.ndim >= 2 else 0
+            if n_samples < ensure_min_samples:
+                raise ValueError(
+                    f"Found array with {n_samples} sample(s) "
+                    f"(shape={arr.shape}) while a minimum of "
+                    f"{ensure_min_samples} is required."
+                )
+            if n_features < ensure_min_features:
+                raise ValueError(
+                    f"Found array with {n_features} feature(s) "
+                    f"(shape={arr.shape}) while a minimum of "
+                    f"{ensure_min_features} is required."
+                )
+
+        if np.iscomplexobj(arr):
+            raise ValueError("Complex data not supported")
+
+        if dtype == "numeric":
+            try:
+                arr = arr.astype(float, copy=False)
+            except (ValueError, TypeError) as exc:
+                raise ValueError(
+                    "Unable to convert array of strings to floating numbers."
+                ) from exc
+
+        if ensure_2d:
+            if arr.shape[0] < ensure_min_samples:
+                raise ValueError(
+                    f"Found array with {arr.shape[0]} sample(s) "
+                    f"(shape={arr.shape}) while a minimum of "
+                    f"{ensure_min_samples} is required."
+                )
+            if arr.shape[1] < ensure_min_features:
+                raise ValueError(
+                    f"Found array with {arr.shape[1]} feature(s) "
+                    f"(shape={arr.shape}) while a minimum of "
+                    f"{ensure_min_features} is required."
+                )
+
+        return arr
diff --git a/PyNomaly/_distance.py b/PyNomaly/_distance.py
index cc9ed58..98a1fcb 100644
--- a/PyNomaly/_distance.py
+++ b/PyNomaly/_distance.py
@@ -1,11 +1,12 @@
 # Authors: Valentino Constantinou <vc@valentino.io>
 # License: Apache 2.0
 
-import numpy as np
 import os
 import sys
-from typing import Tuple
 import warnings
+from typing import Tuple
+
+import numpy as np
 
 from PyNomaly._utils import Utils
 
@@ -30,7 +31,7 @@ def _numba_distance_kernel_parallel(clust_points_vector, n_neighbors):
                 for f in range(d_features):
                     diff_val = clust_points_vector[i, f] - clust_points_vector[j, f]
                     d += diff_val * diff_val
-                d = d ** 0.5
+                d = d**0.5
                 idx_max = 0
                 for k in range(1, n_neighbors):
                     if local_distances[i, k] > local_distances[i, idx_max]:
@@ -57,7 +58,7 @@ def _numba_distance_kernel_sequential(clust_points_vector, n_neighbors):
                 for f in range(d_features):
                     diff_val = clust_points_vector[i, f] - clust_points_vector[j, f]
                     d += diff_val * diff_val
-                d = d ** 0.5
+                d = d**0.5
                 idx_max = 0
                 for k in range(1, n_neighbors):
                     if local_distances[i, k] > local_distances[i, idx_max]:
@@ -91,21 +92,22 @@ def _euclidean(vector1: np.ndarray, vector2: np.ndarray) -> np.ndarray:
         return np.dot(diff, diff) ** 0.5
 
     def _assign_distances(self, data_store: np.ndarray) -> np.ndarray:
-        """
-        Takes a distance matrix, produced by _distances or provided through
-        user input, and assigns distances for each observation to the storage
-        matrix, data_store.
-        :param data_store: the storage matrix that collects information on
-        each observation.
-        :return: the updated storage matrix that collects information on
-        each observation.
-        """
-        for vec, cluster_id in zip(
-            range(self.distance_matrix.shape[0]), self._cluster_labels()
-        ):
+        _dist = getattr(
+            self, "distance_matrix_", getattr(self, "distance_matrix", None)
+        )
+        _neigh = getattr(
+            self, "neighbor_matrix_", getattr(self, "neighbor_matrix", None)
+        )
+
+        if _dist is None or _neigh is None:
+            raise ValueError(
+                "Distance or neighbor matrix is missing. Ensure inputs are valid."
+            )
+
+        for vec, cluster_id in zip(range(_dist.shape[0]), self._cluster_labels()):
             data_store[vec][0] = cluster_id
-            data_store[vec][1] = self.distance_matrix[vec]
-            data_store[vec][2] = self.neighbor_matrix[vec]
+            data_store[vec][1] = _dist[vec]
+            data_store[vec][2] = _neigh[vec]
         return data_store
 
     @staticmethod
@@ -141,9 +143,7 @@ def _compute_distance_and_neighbor_matrix(
 
             yield distances, indexes, i
 
-    def _distances_vectorized(
-        self, clusters, distances, indexes, progress_bar
-    ) -> None:
+    def _distances_vectorized(self, clusters, distances, indexes, progress_bar) -> None:
         """Vectorized kNN distance computation with chunked progress."""
         progress = "="
         total_points = sum(cv.shape[0] for cv, _ in clusters)
@@ -161,21 +161,19 @@ def _distances_vectorized(
                 chunk = clust_points_vector[chunk_start:chunk_end]
 
                 if _scipy_cdist is not None:
-                    dist = _scipy_cdist(
-                        chunk, clust_points_vector, metric="euclidean"
-                    )
+                    dist = _scipy_cdist(chunk, clust_points_vector, metric="euclidean")
                 else:
                     diff = (
-                        chunk[:, np.newaxis, :]
-                        - clust_points_vector[np.newaxis, :, :]
+                        chunk[:, np.newaxis, :] - clust_points_vector[np.newaxis, :, :]
                     )
-                    dist = np.sqrt((diff ** 2).sum(axis=2))
+                    dist = np.sqrt((diff**2).sum(axis=2))
 
                 row_idx = np.arange(chunk_end - chunk_start)
                 dist[row_idx, row_idx + chunk_start] = np.inf
 
-                knn_idx = np.argpartition(dist, self.n_neighbors, axis=1)[
-                    :, : self.n_neighbors
+                n_neighbors = self._effective_n_neighbors()
+                knn_idx = np.argpartition(dist, n_neighbors, axis=1)[
+                    :, :n_neighbors
                 ]
                 knn_dists = np.take_along_axis(dist, knn_idx, axis=1)
 
@@ -201,16 +199,15 @@ def _distances_numba(
                 clust_points_vector = clust_points_vector.reshape(-1, 1)
 
             local_dists, local_idxs = kernel(
-                clust_points_vector.astype(np.float64), self.n_neighbors
+                clust_points_vector.astype(np.float64),
+                self._effective_n_neighbors(),
             )
 
             distances[global_indices] = local_dists
             indexes[global_indices] = global_indices[local_idxs]
 
             if progress_bar:
-                progress = Utils.emit_progress_bar(
-                    progress, idx + 1, len(clusters)
-                )
+                progress = Utils.emit_progress_bar(progress, idx + 1, len(clusters))
 
     def _distances(self, progress_bar: bool = False) -> None:
         """
@@ -221,33 +218,34 @@ def _distances(self, progress_bar: bool = False) -> None:
         :return: the updated storage matrix that collects information on
         each observation.
         """
+        n_neighbors = self._effective_n_neighbors()
         distances = np.full(
-            [self._n_observations(), self.n_neighbors], 9e10, dtype=float
-        )
-        indexes = np.full(
-            [self._n_observations(), self.n_neighbors], 9e10, dtype=float
+            [self._n_observations(), n_neighbors], 9e10, dtype=float
         )
-        self.points_vector = self._convert_to_array(self.data)
+        indexes = np.full([self._n_observations(), n_neighbors], 9e10, dtype=float)
+
+        _data = getattr(self, "data_", getattr(self, "data", None))
+        self.points_vector_ = np.asarray(_data, dtype=float)
 
         cluster_labels = self._cluster_labels()
         cluster_ids = sorted(set(cluster_labels))
 
         clusters = []
         for cluster_id in cluster_ids:
-            indices = np.where(cluster_labels == cluster_id)
+            indices = np.where(cluster_labels == cluster_id)[0]
             clust_points_vector = np.array(
-                self.points_vector.take(indices, axis=0)[0], dtype=np.float64
+                self.points_vector_[indices], dtype=np.float64
             )
-            clusters.append((clust_points_vector, indices[0]))
+            clusters.append((clust_points_vector, indices))
 
-        n_jobs = self.n_jobs
+        n_jobs = getattr(self, "_fit_n_jobs", self.n_jobs)
         if n_jobs == -1:
             n_jobs = os.cpu_count() or 1
 
-        if self.use_numba:
+        use_numba = getattr(self, "_fit_use_numba", self.use_numba)
+        if use_numba:
             self._distances_numba(
-                clusters, distances, indexes, progress_bar,
-                parallel=(n_jobs > 1)
+                clusters, distances, indexes, progress_bar, parallel=(n_jobs > 1)
             )
         else:
             if n_jobs > 1:
@@ -257,9 +255,8 @@ def _distances(self, progress_bar: bool = False) -> None:
                     "to enable parallelism. Falling back to sequential.",
                     UserWarning,
                 )
-            self._distances_vectorized(
-                clusters, distances, indexes, progress_bar
-            )
+            self._distances_vectorized(clusters, distances, indexes, progress_bar)
 
-        self.distance_matrix = distances
-        self.neighbor_matrix = indexes
+        # Assure that computed data is strictly stored in private fitted attributes
+        self.distance_matrix_ = distances
+        self.neighbor_matrix_ = indexes
diff --git a/PyNomaly/_pipeline.py b/PyNomaly/_pipeline.py
index 379d8fb..b38312c 100644
--- a/PyNomaly/_pipeline.py
+++ b/PyNomaly/_pipeline.py
@@ -2,6 +2,7 @@
 # License: Apache 2.0
 
 from math import erf, sqrt
+
 import numpy as np
 
 try:
@@ -94,11 +95,18 @@ def _local_outlier_probability(
     def _n_observations(self) -> int:
         """
         Calculates the number of observations in the data.
-        :return: the number of observations in the input data.
+        :return: the number of observations in the input data if data.
         """
-        if self.data is not None:
-            return len(self.data)
-        return len(self.distance_matrix)
+        _data = getattr(self, "data_", getattr(self, "data", None))
+        _dist = getattr(
+            self, "distance_matrix_", getattr(self, "distance_matrix", None)
+        )
+
+        if _data is not None:
+            return np.asarray(_data).shape[0]
+        if _dist is not None:
+            return np.asarray(_dist).shape[0]
+        return 0
 
     def _store(self) -> np.ndarray:
         """
@@ -115,11 +123,21 @@ def _cluster_labels(self) -> np.ndarray:
         points belong to the same cluster.
         :return: a numpy array of cluster labels.
         """
-        if self.cluster_labels is None:
-            if self.data is not None:
-                return np.array([0] * len(self.data))
-            return np.array([0] * len(self.distance_matrix))
-        return np.array(self.cluster_labels)
+        _labels = getattr(
+            self, "cluster_labels_", getattr(self, "cluster_labels", None)
+        )
+        _data = getattr(self, "data_", getattr(self, "data", None))
+        _dist = getattr(
+            self, "distance_matrix_", getattr(self, "distance_matrix", None)
+        )
+
+        if _labels is None:
+            if _data is not None:
+                return np.array([0] * np.asarray(_data).shape[0])
+            if _dist is not None:
+                return np.array([0] * np.asarray(_dist).shape[0])
+            return np.array([0])
+        return np.array(_labels)
 
     def _ssd(self, data_store: np.ndarray) -> np.ndarray:
         """
@@ -130,9 +148,9 @@ def _ssd(self, data_store: np.ndarray) -> np.ndarray:
         :return: the updated storage matrix that collects information on
         each observation.
         """
-        self.cluster_labels_u = np.unique(data_store[:, 0])
+        self._cluster_labels_u = np.unique(data_store[:, 0])
         ssd_array = np.empty([self._n_observations(), 1])
-        for cluster_id in self.cluster_labels_u:
+        for cluster_id in self._cluster_labels_u:
             indices = np.where(data_store[:, 0] == cluster_id)
             cluster_distances = np.take(data_store[:, 1], indices).tolist()
             ssd = np.power(cluster_distances[0], 2).sum(axis=1)
@@ -152,7 +170,7 @@ def _standard_distances(self, data_store: np.ndarray) -> np.ndarray:
         each observation.
         """
         ssd_vals = data_store[:, 3].astype(float)
-        std_distances = np.sqrt(ssd_vals / self.n_neighbors)
+        std_distances = np.sqrt(ssd_vals / self._effective_n_neighbors())
         return np.hstack((data_store, std_distances.reshape(-1, 1)))
 
     def _prob_distances(self, data_store: np.ndarray) -> np.ndarray:
@@ -178,7 +196,7 @@ def _prob_distances_ev(self, data_store) -> np.ndarray:
         each observation.
         """
         prob_set_distance_ev = np.empty([self._n_observations(), 1])
-        for cluster_id in self.cluster_labels_u:
+        for cluster_id in self._cluster_labels_u:
             indices = np.where(data_store[:, 0] == cluster_id)[0]
             for index in indices:
                 nbrhood = data_store[index][2].astype(int)
@@ -190,7 +208,7 @@ def _prob_distances_ev(self, data_store) -> np.ndarray:
                 ]
                 prob_set_distance_ev[index] = nbrhood_prob_distances_nonan.mean()
 
-        self.prob_distances_ev = prob_set_distance_ev
+        self.prob_distances_ev_ = prob_set_distance_ev
         return np.hstack((data_store, prob_set_distance_ev))
 
     def _prob_local_outlier_factors(self, data_store: np.ndarray) -> np.ndarray:
@@ -229,7 +247,7 @@ def _prob_local_outlier_factors_ev(self, data_store: np.ndarray) -> np.ndarray:
         each observation.
         """
         prob_local_outlier_factor_ev_dict = {}
-        for cluster_id in self.cluster_labels_u:
+        for cluster_id in self._cluster_labels_u:
             indices = np.where(data_store[:, 0] == cluster_id)
             prob_local_outlier_factors = np.take(data_store[:, 7], indices).astype(
                 float
diff --git a/PyNomaly/_validation.py b/PyNomaly/_validation.py
index b1f0a76..5d209a9 100644
--- a/PyNomaly/_validation.py
+++ b/PyNomaly/_validation.py
@@ -1,71 +1,211 @@
 # Authors: Valentino Constantinou <vc@valentino.io>
 # License: Apache 2.0
 
-import numpy as np
-from typing import Union
+"""Input validation helpers for LocalOutlierProbability."""
+
+from __future__ import annotations
+
 import warnings
+from functools import wraps
+from typing import TYPE_CHECKING, Union
+
+import numpy as np
 
+from PyNomaly._compat import check_array
 from PyNomaly.exceptions import ClusterSizeError, MissingValuesError
 
+if TYPE_CHECKING:  # pragma: no cover
+    import pandas as pd  # noqa: F401 (type hints only; pandas is optional)
+
+try:
+    from scipy.sparse import issparse as _issparse
+
+    _SCIPY_AVAILABLE = True
+except ImportError:  # pragma: no cover
+
+    _SCIPY_AVAILABLE = False
+
+    def _issparse(x):
+        return False
+
+
+def _is_sparse_like(obj) -> bool:
+    """Return True for scipy sparse matrices (with or without scipy installed)."""
+    if _issparse(obj):
+        return True
+    obj_type = type(obj)
+    module = getattr(obj_type, "__module__", "") or ""
+    return module.startswith("scipy.sparse")
+
+
+def _require_scipy_for_sparse(obj) -> None:
+    """Raise a clear error when sparse input is passed without scipy."""
+    if _is_sparse_like(obj) and not _SCIPY_AVAILABLE:
+        raise ImportError(
+            "Sparse matrix input requires scipy. "
+            "Install it with: pip install scipy"
+        )
+
 
 class ValidationMixin:
     """Mixin providing input validation methods for LocalOutlierProbability."""
 
+    def _validate_data_matrix(
+        self,
+        X,
+        *,
+        ensure_2d: bool = True,
+        reset: bool = False,
+    ) -> np.ndarray:
+        """
+        Validate and densify feature matrices for fit / predict / stream.
+
+        Sparse inputs are accepted and converted to dense arrays. This keeps
+        the public API sparse-capable while the LoOP kernels operate on dense
+        float64 data.
+        """
+        if X is None:
+            raise ValueError("Input data must be provided.")
+
+        # Preserve DataFrame column names before densifying / converting.
+        feature_names = None
+        if hasattr(X, "columns"):
+            try:
+                feature_names = np.asarray(X.columns, dtype=object)
+            except Exception:
+                feature_names = None
+
+        if _is_sparse_like(X):
+            _require_scipy_for_sparse(X)
+            X = X.toarray()
+
+        # Historical LoOP API accepts 1-d series and reshapes to a single feature.
+        arr_probe = np.asarray(X)
+        if getattr(arr_probe, "ndim", None) == 1:
+            if arr_probe.__class__.__name__ != "ndarray" or not isinstance(X, np.ndarray):
+                warnings.warn(
+                    "Provided data or distance matrix must be in ndarray or DataFrame.",
+                    UserWarning,
+                )
+            X = np.asarray(arr_probe, dtype=float).reshape(-1, 1)
+
+        # Prefer PyNomaly's MissingValuesError over sklearn's generic ValueError.
+        try:
+            probe = np.asarray(X, dtype=float)
+        except (ValueError, TypeError):
+            probe = None
+        if probe is not None and (np.any(np.isnan(probe)) or np.any(np.isinf(probe))):
+            raise MissingValuesError(
+                "Input contains NaN, infinity or a value too large for "
+                "dtype('float64')."
+            )
+
+        arr = check_array(
+            X,
+            accept_sparse=False,
+            dtype="numeric",
+            ensure_2d=ensure_2d,
+            ensure_min_samples=1,
+            ensure_min_features=1,
+        )
+        arr = np.asarray(arr, dtype=float)
+
+        if np.any(np.isnan(arr)) or np.any(np.isinf(arr)):
+            raise MissingValuesError(
+                "Input contains NaN, infinity or a value too large for "
+                "dtype('float64')."
+            )
+
+        if reset:
+            self.n_features_in_ = arr.shape[1] if arr.ndim == 2 else 1
+            if feature_names is not None and len(feature_names) == self.n_features_in_:
+                self.feature_names_in_ = feature_names
+
+        return arr
+
     @staticmethod
     def _convert_to_array(obj: Union["pd.DataFrame", np.ndarray]) -> np.ndarray:
         """
         Converts the input data to a numpy array if it is a Pandas DataFrame
         or validates it is already a numpy array.
-        :param obj: user-provided input data.
-        :return: a vector of values to be used in calculating the local
-        outlier probability.
         """
-        if obj.__class__.__name__ == "DataFrame":
-            points_vector = obj.values
-            return points_vector
+        if _is_sparse_like(obj):
+            _require_scipy_for_sparse(obj)
+            obj = obj.toarray()
+
+        if obj.__class__.__name__ in ("DataFrame", "Series"):
+            arr = np.asarray(obj.values)
         elif obj.__class__.__name__ == "ndarray":
-            points_vector = obj
-            return points_vector
+            arr = obj
         else:
             warnings.warn(
-                "Provided data or distance matrix must be in ndarray "
-                "or DataFrame.",
+                "Provided data or distance matrix must be in ndarray or DataFrame.",
                 UserWarning,
             )
-            if isinstance(obj, list):
-                points_vector = np.array(obj)
-                return points_vector
-            points_vector = np.array([obj])
-            return points_vector
+            arr = np.asarray(obj, dtype=float)
+            if arr.ndim == 0:
+                arr = np.array([obj], dtype=float)
+
+        if arr.size == 0:
+            if arr.ndim >= 2 and arr.shape[1] == 0:
+                raise ValueError(
+                    f"Found array with 0 feature(s) (shape={arr.shape}) while a "
+                    f"minimum of 1 is required."
+                )
+            raise ValueError(
+                f"Found array with 0 sample(s) (shape={arr.shape}) while a "
+                f"minimum of 1 is required."
+            )
+
+        if np.iscomplexobj(arr):
+            raise ValueError("Complex data not supported.")
+
+        return arr.astype(float, copy=False)
+
+    @staticmethod
+    def _convert_observation(obj: Union["pd.DataFrame", np.ndarray]) -> np.ndarray:
+        """
+        Convert a single observation to a 1-D float vector.
+
+        Accepts dense arrays, pandas objects, lists, and scipy sparse rows.
+        Sparse matrices are densified (requires scipy).
+        """
+        arr = ValidationMixin._convert_to_array(obj)
+        if arr.ndim == 2 and arr.shape[0] == 1:
+            return arr.ravel()
+        return arr
 
     def _validate_inputs(self):
         """
-        Validates the inputs provided during initialization to ensure
-        that the needed objects are provided.
-        :return: a tuple of (data, distance_matrix, neighbor_matrix) or
-        raises a warning for invalid inputs.
+        Validates that either data or a distance/neighbor matrix pair is set.
         """
-        if all(v is None for v in [self.data, self.distance_matrix]):
-            warnings.warn(
-                "Data or a distance matrix must be provided.", UserWarning
-            )
+        _data = getattr(self, "data_", getattr(self, "data", None))
+        _dist = getattr(
+            self, "distance_matrix_", getattr(self, "distance_matrix", None)
+        )
+        _neigh = getattr(
+            self, "neighbor_matrix_", getattr(self, "neighbor_matrix", None)
+        )
+
+        if all(v is None for v in [_data, _dist]):
+            warnings.warn("Data or a distance matrix must be provided.", UserWarning)
             return False
-        elif all(v is not None for v in [self.data, self.distance_matrix]):
+        if all(v is not None for v in [_data, _dist]):
             warnings.warn(
                 "Only one of the following may be provided: data or a "
                 "distance matrix (not both).",
                 UserWarning,
             )
             return False
-        if self.data is not None:
-            points_vector = self._convert_to_array(self.data)
-            return points_vector, self.distance_matrix, self.neighbor_matrix
-        if all(
-            matrix is not None
-            for matrix in [self.neighbor_matrix, self.distance_matrix]
-        ):
-            dist_vector = self._convert_to_array(self.distance_matrix)
-            neigh_vector = self._convert_to_array(self.neighbor_matrix)
+
+        if _data is not None:
+            points_vector = self._validate_data_matrix(_data, ensure_2d=True, reset=True)
+            self.data_ = points_vector
+            return points_vector, _dist, _neigh
+
+        if all(matrix is not None for matrix in [_neigh, _dist]):
+            dist_vector = self._convert_to_array(_dist)
+            neigh_vector = self._convert_to_array(_neigh)
         else:
             warnings.warn(
                 "A neighbor index matrix and distance matrix must both be "
@@ -73,37 +213,34 @@ def _validate_inputs(self):
                 UserWarning,
             )
             return False
-        if self.distance_matrix.shape != self.neighbor_matrix.shape:
+
+        if _dist.shape != _neigh.shape:
             warnings.warn(
-                "The shape of the distance and neighbor "
-                "index matrices must match.",
+                "The shape of the distance and neighbor index matrices must match.",
                 UserWarning,
             )
             return False
-        elif (self.distance_matrix.shape[1] != self.n_neighbors) or (
-            self.neighbor_matrix.shape[1] != self.n_neighbors
+        if (_dist.shape[1] != self.n_neighbors) or (
+            _neigh.shape[1] != self.n_neighbors
         ):
             warnings.warn(
-                "The shape of the distance or "
-                "neighbor index matrix does not "
-                "match the number of neighbors "
-                "specified.",
+                "The shape of the distance or neighbor index matrix does not "
+                "match the number of neighbors specified.",
                 UserWarning,
             )
             return False
-        return self.data, dist_vector, neigh_vector
+
+        self.distance_matrix_ = dist_vector
+        self.neighbor_matrix_ = neigh_vector
+        return _data, dist_vector, neigh_vector
 
     def _check_cluster_size(self) -> None:
-        """
-        Validates the cluster labels to ensure that the smallest cluster
-        size (number of observations in the cluster) is larger than the
-        specified number of neighbors.
-        :raises ClusterSizeError: if any cluster is too small.
-        """
+        """Raises ClusterSizeError if any cluster is too small for n_neighbors."""
+        n_neighbors = getattr(self, "n_neighbors_", self.n_neighbors)
         c_labels = self._cluster_labels()
         for cluster_id in set(c_labels):
             c_size = np.where(c_labels == cluster_id)[0].shape[0]
-            if c_size <= self.n_neighbors:
+            if c_size <= n_neighbors:
                 raise ClusterSizeError(
                     "Number of neighbors specified larger than smallest "
                     "cluster. Specify a number of neighbors smaller than "
@@ -113,59 +250,57 @@ def _check_cluster_size(self) -> None:
 
     def _check_n_neighbors(self) -> bool:
         """
-        Validates the specified number of neighbors to ensure that it is
-        greater than 0 and that the specified value is less than the total
-        number of observations.
-        :return: a boolean indicating whether validation has passed without
-        adjustment.
+        Validate n_neighbors without mutating the public parameter.
+
+        The effective value used during fit is stored on ``n_neighbors_``.
         """
-        if not self.n_neighbors > 0:
-            self.n_neighbors = 10
+        n_neighbors = self.n_neighbors
+        n_obs = self._n_observations()
+
+        if not n_neighbors > 0:
             warnings.warn(
-                "n_neighbors must be greater than 0."
-                " Fit with " + str(self.n_neighbors) + " instead.",
+                "n_neighbors must be greater than 0. Fit with 10 instead.",
                 UserWarning,
             )
-            return False
-        elif self.n_neighbors >= self._n_observations():
-            self.n_neighbors = self._n_observations() - 1
+            n_neighbors = 10
+            ok = False
+        else:
+            ok = True
+
+        if n_neighbors >= n_obs:
+            adjusted = max(n_obs - 1, 1)
             warnings.warn(
                 "n_neighbors must be less than the number of observations."
-                " Fit with " + str(self.n_neighbors) + " instead.",
+                f" Fit with {adjusted} instead.",
                 UserWarning,
             )
-        return True
+            n_neighbors = adjusted
+            ok = False
+
+        self.n_neighbors_ = n_neighbors
+        return ok
 
     def _check_extent(self) -> bool:
-        """
-        Validates the specified extent parameter to ensure it is either 1,
-        2, or 3.
-        :return: a boolean indicating whether validation has passed.
-        """
+        """Validate that extent is 1, 2, or 3."""
         if self.extent not in [1, 2, 3]:
-            warnings.warn(
-                "extent parameter (lambda) must be 1, 2, or 3.", UserWarning
-            )
+            warnings.warn("extent parameter (lambda) must be 1, 2, or 3.", UserWarning)
             return False
         return True
 
     def _check_missing_values(self) -> None:
-        """
-        Validates the provided data to ensure that it contains no
-        missing values.
-        :raises MissingValuesError: if data contains NaN values.
-        """
-        if np.any(np.isnan(self.data)):
-            raise MissingValuesError(
-                "Method does not support missing values in input data."
-            )
+        """Raise MissingValuesError if fitted data contains NaN/Inf."""
+        _data = getattr(self, "data_", getattr(self, "data", None))
+        if _data is not None:
+            arr = self._convert_to_array(_data)
+            if np.any(np.isnan(arr)) or np.any(np.isinf(arr)):
+                raise MissingValuesError(
+                    "Input contains NaN, infinity or a value too large for "
+                    "dtype('float64')."
+                )
 
     def _check_is_fit(self) -> bool:
-        """
-        Checks that the model was fit prior to calling the stream() method.
-        :return: a boolean indicating whether the model has been fit.
-        """
-        if self.is_fit is False:
+        """Return whether the estimator has been successfully fit."""
+        if getattr(self, "is_fit_", False) is False:
             warnings.warn(
                 "Must fit on historical data by calling fit() prior to "
                 "calling stream(x).",
@@ -175,13 +310,7 @@ def _check_is_fit(self) -> bool:
         return True
 
     def _check_no_cluster_labels(self) -> bool:
-        """
-        Checks to see if cluster labels are attempting to be used in
-        stream() and, if so, returns False. As PyNomaly does not accept
-        clustering algorithms as input, the stream approach does not
-        support clustering.
-        :return: a boolean indicating whether single cluster (no labels).
-        """
+        """Return False when multi-cluster labels are present (stream unsupported)."""
         if len(set(self._cluster_labels())) > 1:
             warnings.warn(
                 "Stream approach does not support clustered data. "
@@ -191,61 +320,42 @@ def _check_no_cluster_labels(self) -> bool:
             return False
         return True
 
+    def _effective_n_neighbors(self) -> int:
+        """Return the neighbor count used by the current fit."""
+        return getattr(self, "n_neighbors_", self.n_neighbors)
+
 
-def accepts(*types):
+def validate_init_types(*types):
     """
-    A decorator that facilitates a form of type checking for the inputs
-    which can be used in Python 3.4-3.7 in lieu of Python 3.5+'s type
-    hints.
-    :param types: the input types of the objects being passed as arguments
-    in __init__.
-    :return: a decorator.
+    Soft type checker used at fit-time for legacy constructor parameters.
     """
 
     def decorator(f):
-        assert len(types) == f.__code__.co_argcount
-
-        def new_f(*args, **kwds):
-            for a, t in zip(args, types):
-                if type(a).__name__ == "DataFrame":
-                    a = np.array(a)
-                if isinstance(a, t) is False:
-                    warnings.warn(
-                        "Argument %r is not of type %s" % (a, t), UserWarning
-                    )
-            opt_types = {
-                "extent": {"type": (int, np.integer)},
-                "n_neighbors": {"type": (int, np.integer)},
-                "use_numba": {"type": bool},
-                "n_jobs": {"type": (int, np.integer)},
-                "progress_bar": {"type": bool},
-                "data": {"type": np.ndarray},
-                "distance_matrix": {"type": np.ndarray},
-                "neighbor_matrix": {"type": np.ndarray},
-                "cluster_labels": {"type": (list, np.ndarray)},
-            }
-            for x in kwds:
-                if x in opt_types:
-                    v = kwds[x]
-                    if type(v).__name__ == "DataFrame":
-                        v = np.array(v)
-                    opt_types[x]["value"] = v
-            for k in opt_types:
-                try:
-                    if (
-                        isinstance(opt_types[k]["value"], opt_types[k]["type"])
-                        is False
-                    ):
+        @wraps(f)
+        def new_f(self, *args, **kwds):
+            attr_names = [
+                "extent",
+                "n_neighbors",
+                "use_numba",
+                "n_jobs",
+                "progress_bar",
+                "data",
+                "distance_matrix",
+                "neighbor_matrix",
+                "cluster_labels",
+            ]
+            for attr, expected_type in zip(attr_names, types[1:]):
+                val = getattr(self, attr, None)
+                if val is not None:
+                    if type(val).__name__ == "DataFrame":
+                        val = np.array(val)
+                    if not isinstance(val, expected_type):
                         warnings.warn(
-                            "Argument %r is not of type %s."
-                            % (k, opt_types[k]["type"]),
+                            "Argument %r is not of type %s." % (attr, expected_type),
                             UserWarning,
                         )
-                except KeyError:
-                    pass
-            return f(*args, **kwds)
+            return f(self, *args, **kwds)
 
-        new_f.__name__ = f.__name__
         return new_f
 
     return decorator
diff --git a/PyNomaly/exceptions.py b/PyNomaly/exceptions.py
index 7417ddb..7e04eff 100644
--- a/PyNomaly/exceptions.py
+++ b/PyNomaly/exceptions.py
@@ -4,19 +4,23 @@
 
 class PyNomalyError(Exception):
     """Base exception for PyNomaly."""
+
     pass
 
 
 class ValidationError(PyNomalyError):
     """Raised when input validation fails."""
+
     pass
 
 
 class ClusterSizeError(ValidationError):
     """Raised when cluster size is smaller than n_neighbors."""
+
     pass
 
 
-class MissingValuesError(ValidationError):
+class MissingValuesError(ValidationError, ValueError):
     """Raised when data contains missing values."""
+
     pass
diff --git a/PyNomaly/loop.py b/PyNomaly/loop.py
index a34c518..cabd820 100644
--- a/PyNomaly/loop.py
+++ b/PyNomaly/loop.py
@@ -1,78 +1,108 @@
 # Authors: Valentino Constantinou <vc@valentino.io>
 # License: Apache 2.0
 
-import numpy as np
+"""Local Outlier Probability (LoOP) estimator."""
+
+from __future__ import annotations
+
 import sys
 import warnings
 
-from PyNomaly.exceptions import (
-    PyNomalyError,
-    ValidationError,
-    ClusterSizeError,
-    MissingValuesError,
+import numpy as np
+
+from PyNomaly._compat import (
+    BaseEstimator,
+    OutlierMixin,
+    check_array,
+    check_is_fitted,
 )
-from PyNomaly._utils import Utils
-from PyNomaly._validation import ValidationMixin, accepts
 from PyNomaly._distance import DistanceMixin
 from PyNomaly._pipeline import PipelineMixin
+from PyNomaly._validation import ValidationMixin, validate_init_types
+from PyNomaly.exceptions import (  # noqa: F401 (re-exported for backward compatibility)
+    ClusterSizeError,
+    MissingValuesError,
+    PyNomalyError,
+    ValidationError,
+)
 
 __author__ = "Valentino Constantinou"
 __version__ = "1.0.0"
 __license__ = "Apache License, Version 2.0"
 
 
-class LocalOutlierProbability(ValidationMixin, DistanceMixin, PipelineMixin):
+class LocalOutlierProbability(
+    OutlierMixin, ValidationMixin, DistanceMixin, PipelineMixin, BaseEstimator
+):
     """Local Outlier Probability (LoOP) estimator.
 
-    :param extent: an integer value [1, 2, 3] that controls the statistical
-    extent, e.g. lambda times the standard deviation from the mean (optional,
-    default 3)
-    :param n_neighbors: the total number of neighbors to consider w.r.t. each
-    sample (optional, default 10)
-    :param use_numba: whether to use Numba JIT acceleration for distance
-    computation (optional, default False)
-    :param n_jobs: controls Numba thread-level parallelism via prange.
-    Use -1 to use all available CPU cores, or 1 for sequential processing.
-    Only effective when use_numba=True (optional, default 1)
-    :param progress_bar: whether to display a progress bar during distance
-    computation (optional, default False)
-
-    Based on the work of Kriegel, Kröger, Schubert, and Zimek (2009) in LoOP:
-    Local Outlier Probabilities.
+    Parameters
     ----------
+    extent : int, default=3
+        Statistical extent (lambda); must be 1, 2, or 3.
+    n_neighbors : int, default=10
+        Number of neighbors considered for each sample.
+    use_numba : bool, default=False
+        Whether to use Numba JIT acceleration for distance computation.
+    n_jobs : int, default=1
+        Numba parallelism. ``-1`` uses all cores. Only effective with
+        ``use_numba=True``.
+    progress_bar : bool, default=False
+        Whether to display a progress bar during distance computation.
+
+    Notes
+    -----
+    Fitted attributes end with ``_`` (scikit-learn convention). Legacy
+    read aliases without the trailing underscore remain available for
+    *output* attributes such as ``local_outlier_probabilities``. Init /
+    configuration attributes (e.g. ``distance_matrix``) are never aliased
+    to fitted state — that would break ``get_params`` / ``check_estimator``.
+
+    Sparse matrices (``scipy.sparse``) are accepted for ``fit``, ``predict``,
+    ``decision_function``, and ``stream``. They are converted to dense arrays
+    internally before distance computation. **scipy** is a soft dependency for
+    sparse input — install it with ``pip install scipy`` (included in the
+    ``PyNomaly[all]`` extra).
 
     References
     ----------
     .. [1] Breunig M., Kriegel H.-P., Ng R., Sander, J. LOF: Identifying
-           Density-based Local Outliers. ACM SIGMOD
-           International Conference on Management of Data (2000).
+           Density-based Local Outliers. ACM SIGMOD (2000).
     .. [2] Kriegel H.-P., Kröger P., Schubert E., Zimek A. LoOP: Local Outlier
-           Probabilities. 18th ACM conference on
-           Information and knowledge management, CIKM (2009).
-    .. [3] Goldstein M., Uchida S. A Comparative Evaluation of Unsupervised
-           Anomaly Detection Algorithms for Multivariate Data. PLoS ONE 11(4):
-           e0152173 (2016).
-    .. [4] Hamlet C., Straub J., Russell M., Kerlin S. An incremental and
-           approximate local outlier probability algorithm for intrusion
-           detection and its evaluation. Journal of Cyber Security Technology
-           (2016).
+           Probabilities. CIKM (2009).
     """
 
-    _DATA_PARAMS = ("data", "distance_matrix", "neighbor_matrix",
-                    "cluster_labels")
+    _DATA_PARAMS = ("data", "distance_matrix", "neighbor_matrix", "cluster_labels")
+
+    # ------------------------------------------------------------------
+    # Backward-compatible read aliases for *fitted output* attributes.
+    # These are NOT constructor parameters, so they do not affect
+    # get_params() / set_params() / check_estimator param comparisons.
+    # ------------------------------------------------------------------
+    @property
+    def points_vector(self):
+        return getattr(self, "points_vector_", None)
+
+    @property
+    def prob_distances(self):
+        return getattr(self, "prob_distances_", None)
+
+    @property
+    def prob_distances_ev(self):
+        return getattr(self, "prob_distances_ev_", None)
+
+    @property
+    def norm_prob_local_outlier_factor(self):
+        return getattr(self, "norm_prob_local_outlier_factor_", None)
+
+    @property
+    def local_outlier_probabilities(self):
+        return getattr(self, "local_outlier_probabilities_", None)
+
+    @property
+    def is_fit(self):
+        return getattr(self, "is_fit_", False)
 
-    @accepts(
-        object,
-        (int, np.integer),
-        (int, np.integer),
-        bool,
-        (int, np.integer),
-        bool,
-        np.ndarray,
-        np.ndarray,
-        np.ndarray,
-        list,
-    )
     def __init__(
         self,
         extent=3,
@@ -85,22 +115,26 @@ def __init__(
         neighbor_matrix=None,
         cluster_labels=None,
     ) -> None:
+        # Configuration parameters only — never mutate these in fit().
         self.extent = extent
         self.n_neighbors = n_neighbors
         self.use_numba = use_numba
         self.n_jobs = n_jobs
         self.progress_bar = progress_bar
 
-        # Emit deprecation warnings for data params passed to __init__
-        _locals = {"data": data, "distance_matrix": distance_matrix,
-                   "neighbor_matrix": neighbor_matrix,
-                   "cluster_labels": cluster_labels}
+        # Deprecated data-in-constructor path (backward compatible).
+        _locals = {
+            "data": data,
+            "distance_matrix": distance_matrix,
+            "neighbor_matrix": neighbor_matrix,
+            "cluster_labels": cluster_labels,
+        }
         for param in self._DATA_PARAMS:
             if _locals[param] is not None:
                 warnings.warn(
-                    "Passing '{}' to __init__ is deprecated. "
+                    f"Passing '{param}' to __init__ is deprecated. "
                     "Pass it to fit() instead. This will raise an error "
-                    "in a future version.".format(param),
+                    "in a future version.",
                     FutureWarning,
                     stacklevel=2,
                 )
@@ -110,166 +144,333 @@ def __init__(
         self.neighbor_matrix = neighbor_matrix
         self.cluster_labels = cluster_labels
 
-        self.points_vector = None
-        self.prob_distances = None
-        self.prob_distances_ev = None
-        self.norm_prob_local_outlier_factor = None
-        self.local_outlier_probabilities = None
-        self._objects = {}
-        self.is_fit = False
+    def __sklearn_tags__(self):
+        tags = super().__sklearn_tags__()
+        # Sparse inputs are densified in validation; tag as supported.
+        tags.input_tags.sparse = True
+        return tags
 
-        if self.use_numba is True and "numba" not in sys.modules:
-            self.use_numba = False
+    def _reset_state(self) -> None:
+        """Clear fitted state so the estimator can be re-fit."""
+        self.points_vector_ = None
+        self.prob_distances_ = None
+        self.prob_distances_ev_ = None
+        self.norm_prob_local_outlier_factor_ = None
+        self.local_outlier_probabilities_ = None
+        self.data_ = None
+        self.distance_matrix_ = None
+        self.neighbor_matrix_ = None
+        self.cluster_labels_ = None
+        self.n_neighbors_ = self.n_neighbors
+        self._objects_ = {}
+        self.is_fit_ = False
+
+    def _resolve_runtime_flags(self):
+        """Return (use_numba, n_jobs) without mutating public params."""
+        use_numba = self.use_numba
+        if use_numba is True and "numba" not in sys.modules:
             warnings.warn(
-                "Numba is not available, falling back to pure python mode.", UserWarning
+                "Numba is not available, falling back to pure python mode.",
+                UserWarning,
             )
+            use_numba = False
 
-        if self.n_jobs < -1 or self.n_jobs == 0:
+        n_jobs = self.n_jobs
+        if n_jobs < -1 or n_jobs == 0:
             warnings.warn(
                 "n_jobs must be -1 or a positive integer. Defaulting to 1.",
                 UserWarning,
             )
-            self.n_jobs = 1
-
-        self._check_extent()
-
-    def _reset_state(self) -> None:
-        """Resets computed state to allow re-fitting with new data."""
-        self.points_vector = None
-        self.prob_distances = None
-        self.prob_distances_ev = None
-        self.norm_prob_local_outlier_factor = None
-        self.local_outlier_probabilities = None
-        self._objects = {}
-        self.is_fit = False
+            n_jobs = 1
+        return use_numba, n_jobs
 
+    @validate_init_types(
+        object,
+        (int, np.integer),
+        (int, np.integer),
+        bool,
+        (int, np.integer),
+        bool,
+        np.ndarray,
+        np.ndarray,
+        np.ndarray,
+        list,
+    )
     def fit(
         self,
+        X=None,
+        y=None,
         data=None,
         distance_matrix=None,
         neighbor_matrix=None,
         cluster_labels=None,
     ) -> "LocalOutlierProbability":
         """
-        Calculates the local outlier probability for each observation in the
-        input data according to the input parameters extent, n_neighbors, and
-        cluster_labels.
-        :param data: a Pandas DataFrame or Numpy array of float data
-        (optional, default None)
-        :param distance_matrix: a precomputed distance matrix of shape
-        (n_observations, n_neighbors) (optional, default None)
-        :param neighbor_matrix: a precomputed neighbor index matrix of shape
-        (n_observations, n_neighbors) (optional, default None)
-        :param cluster_labels: a numpy array or list of cluster assignments
-        w.r.t. each sample (optional, default None)
-        :return: self, which contains the local outlier probabilities as
-        self.local_outlier_probabilities.
-        :raises ClusterSizeError: if any cluster is smaller than n_neighbors.
-        :raises MissingValuesError: if data contains missing values.
+        Fit the LoOP model.
+
+        Parameters
+        ----------
+        X : array-like of shape (n_samples, n_features), optional
+            Training data (scikit-learn style). Sparse matrices are densified.
+        y : ignored
+            Present for API consistency. For backward compatibility only,
+            a 2-d ``y`` with more than one column is treated as a deprecated
+            positional ``distance_matrix`` when no raw data (``X`` / ``data``)
+            is given.
+        data, distance_matrix, neighbor_matrix, cluster_labels
+            Explicit keyword alternatives to ``X`` / precomputed neighbors.
         """
+        if X is not None:
+            data = X
+
+        # Legacy positional fit(data, distance_matrix) support. Only applies
+        # when no raw data is given and y cannot be a label vector, so that
+        # sklearn-style fit(X, y) with (n, 1) labels is not hijacked.
+        if y is not None and distance_matrix is None and data is None:
+            y_array = np.asarray(y)
+            if y_array.ndim == 2 and y_array.shape[1] > 1:
+                distance_matrix = y
+                y = None
+                warnings.warn(
+                    "Passing 'distance_matrix' as the second positional "
+                    "argument is deprecated. Use keyword arguments instead.",
+                    FutureWarning,
+                )
 
+        use_numba, n_jobs = self._resolve_runtime_flags()
+        self._fit_use_numba = use_numba
+        self._fit_n_jobs = n_jobs
+
+        self._check_extent()
         self._reset_state()
 
-        if data is not None and distance_matrix is not None:
-            self.data = data
-            self.distance_matrix = distance_matrix
-            if neighbor_matrix is not None:
-                self.neighbor_matrix = neighbor_matrix
-        elif data is not None:
-            self.data = data
-            self.distance_matrix = None
-            self.neighbor_matrix = None
-        elif distance_matrix is not None:
-            self.distance_matrix = distance_matrix
-            self.data = None
-            if neighbor_matrix is not None:
-                self.neighbor_matrix = neighbor_matrix
-        elif neighbor_matrix is not None:
-            self.neighbor_matrix = neighbor_matrix
-        if cluster_labels is not None:
-            self.cluster_labels = cluster_labels
+        self.data_ = data if data is not None else self.data
+        self.distance_matrix_ = (
+            distance_matrix if distance_matrix is not None else self.distance_matrix
+        )
+        self.neighbor_matrix_ = (
+            neighbor_matrix if neighbor_matrix is not None else self.neighbor_matrix
+        )
+        self.cluster_labels_ = (
+            cluster_labels if cluster_labels is not None else self.cluster_labels
+        )
 
         if self._validate_inputs() is False:
             return self
+
+        n_obs = self._n_observations()
+        if n_obs < 2:
+            raise ValueError(
+                f"Found array with {n_obs} sample(s) while a minimum of 2 "
+                f"is required by LocalOutlierProbability."
+            )
+
         self._check_n_neighbors()
         self._check_cluster_size()
-        if self.data is not None:
+
+        _data = getattr(self, "data_", None)
+        if _data is not None:
             self._check_missing_values()
 
         store = self._store()
-        if self.data is not None:
+        if _data is not None:
             self._distances(progress_bar=self.progress_bar)
+
         store = self._assign_distances(store)
         store = self._ssd(store)
         store = self._standard_distances(store)
         store = self._prob_distances(store)
-        self.prob_distances = store[:, 5]
+        self.prob_distances_ = np.asarray(store[:, 5], dtype=float)
         store = self._prob_distances_ev(store)
         store = self._prob_local_outlier_factors(store)
         store = self._prob_local_outlier_factors_ev(store)
         store = self._norm_prob_local_outlier_factors(store)
-        self.norm_prob_local_outlier_factor = store[:, 9].max()
+        self.norm_prob_local_outlier_factor_ = float(
+            np.asarray(store[:, 9], dtype=float).max()
+        )
         store = self._local_outlier_probabilities(store)
-        self.local_outlier_probabilities = store[:, 10]
-
-        self.is_fit = True
+        self.local_outlier_probabilities_ = np.asarray(store[:, 10], dtype=float)
 
+        if _data is not None and self.points_vector_ is not None:
+            self.n_features_in_ = (
+                self.points_vector_.shape[1]
+                if self.points_vector_.ndim == 2
+                else 1
+            )
+        elif self.distance_matrix_ is not None:
+            self.n_features_in_ = 1
+
+        # sklearn outlier API: decision_function = score_samples - offset_
+        # score_samples = -LoOP probability (higher is more normal)
+        # offset_ = -0.5 so decision >= 0 iff probability <= 0.5
+        self.offset_ = -0.5
+        self.is_fit_ = True
         return self
 
-    def stream(self, x: np.ndarray) -> np.ndarray:
-        """
-        Calculates the local outlier probability for an individual sample
-        according to the input parameters extent, n_neighbors, and
-        cluster_labels after first calling fit(). Observations are assigned
-        a local outlier probability against the mean of expected values of
-        probabilistic distance and the normalized probabilistic outlier
-        factor from the earlier model, provided when calling fit().
-        distance
-        :param x: an observation to score for its local outlier probability.
-        :return: the local outlier probability of the input observation.
-        """
+    def stream(self, x: np.ndarray) -> float:
+        """Score a single observation against the fitted model."""
+        if self._check_is_fit() is False:
+            self.fit()
+
+        _data = getattr(self, "data_", getattr(self, "data", None))
+        _dist = getattr(
+            self, "distance_matrix_", getattr(self, "distance_matrix", None)
+        )
+        _clust = getattr(self, "cluster_labels_", getattr(self, "cluster_labels", None))
 
         orig_cluster_labels = None
         if self._check_no_cluster_labels() is False:
-            orig_cluster_labels = self.cluster_labels
-            self.cluster_labels = np.array([0] * len(self.data))
+            orig_cluster_labels = _clust
+            if getattr(self, "points_vector_", None) is not None:
+                self.cluster_labels_ = np.array([0] * self.points_vector_.shape[0])
 
-        if self._check_is_fit() is False:
-            self.fit()
+        point_vector = self._convert_observation(x)
+        n_neighbors = self._effective_n_neighbors()
+        distances = np.full([1, n_neighbors], 9e10, dtype=float)
 
-        point_vector = self._convert_to_array(x)
-        distances = np.full([1, self.n_neighbors], 9e10, dtype=float)
-        if self.data is not None:
-            matrix = self.points_vector
+        if _data is not None:
+            matrix = getattr(self, "points_vector_", None)
         else:
-            matrix = self.distance_matrix
-            # When using distance matrix mode, x is a scalar distance value.
-            # Extract scalar from array to avoid NumPy assignment errors.
+            matrix = _dist
             if point_vector.size == 1:
                 point_vector = float(point_vector.flat[0])
+
+        if matrix is None:
+            from PyNomaly._compat import NotFittedError
+
+            raise NotFittedError(
+                "This LocalOutlierProbability instance is not fitted yet. "
+                "Call 'fit' with appropriate arguments before using this "
+                "estimator."
+            )
+
         for p in range(0, matrix.shape[0]):
-            if self.data is not None:
-                d = self._euclidean(matrix[p, :], point_vector)
+            if _data is not None:
+                ref_point = matrix[p] if matrix.ndim == 1 else matrix[p, :]
+                d = float(self._euclidean(ref_point, point_vector))
             else:
-                d = point_vector
+                d = float(point_vector)
+
             idx_max = distances[0].argmax()
             if d < distances[0][idx_max]:
                 distances[0][idx_max] = d
 
         ssd = np.power(distances, 2).sum()
-        std_dist = np.sqrt(np.divide(ssd, self.n_neighbors))
+        std_dist = np.sqrt(np.divide(ssd, n_neighbors))
         prob_dist = self._prob_distance(self.extent, std_dist)
         plof = self._prob_outlier_factor(
-            np.array(prob_dist), np.array(self.prob_distances_ev.mean())
+            np.array(prob_dist), np.array(self.prob_distances_ev_.mean())
         )
-        loop = self._local_outlier_probability(
-            plof, self.norm_prob_local_outlier_factor
+        loop_score = self._local_outlier_probability(
+            plof, self.norm_prob_local_outlier_factor_
         )
 
         if orig_cluster_labels is not None:
-            self.cluster_labels = orig_cluster_labels
+            self.cluster_labels_ = orig_cluster_labels
+
+        return float(np.asarray(loop_score, dtype=float).ravel()[0])
+
+    def _score_matrix(self, X) -> np.ndarray:
+        """
+        Return LoOP probabilities for ``X``.
+
+        Training rows are matched to fitted scores so predict / decision
+        function stay order- and subset-invariant. Unseen rows use ``stream``.
+        """
+        check_is_fitted(self, ["is_fit_"])
+
+        if _issparse_safe(X):
+            from PyNomaly._validation import _require_scipy_for_sparse
+
+            _require_scipy_for_sparse(X)
+            X = X.toarray()
+
+        X_probe = np.asarray(X)
+        if X_probe.ndim == 1 and getattr(self, "n_features_in_", None) == 1:
+            X = X_probe.reshape(-1, 1)
+
+        X_array = check_array(
+            X,
+            accept_sparse=False,
+            dtype="numeric",
+            ensure_2d=True,
+        )
+        X_array = np.asarray(X_array, dtype=float)
+
+        if np.any(np.isnan(X_array)) or np.any(np.isinf(X_array)):
+            raise ValueError(
+                "Input contains NaN, infinity or a value too large for "
+                "dtype('float64')."
+            )
+
+        if hasattr(self, "n_features_in_"):
+            n_features = X_array.shape[1]
+            if n_features != self.n_features_in_:
+                raise ValueError(
+                    f"X has {n_features} features, but {type(self).__name__} "
+                    f"is expecting {self.n_features_in_} features as input."
+                )
+
+        train = getattr(self, "points_vector_", None)
+        fitted = getattr(self, "local_outlier_probabilities_", None)
+        scores = np.empty(X_array.shape[0], dtype=float)
 
-        return loop
+        if train is not None and fitted is not None and train.ndim == 2:
+            for i, row in enumerate(X_array):
+                matches = np.where(np.all(np.isclose(train, row), axis=1))[0]
+                if matches.size:
+                    scores[i] = fitted[matches[0]]
+                else:
+                    scores[i] = float(self.stream(row))
+            return scores
+
+        for i, row in enumerate(X_array):
+            scores[i] = float(self.stream(row))
+        return scores
+
+    def decision_function(self, X):
+        """
+        Signed outlier score (higher is more normal).
+
+        Equivalent to ``score_samples(X) - offset_``. Inliers satisfy
+        ``decision_function(X) >= 0``.
+        """
+        return self.score_samples(X) - self.offset_
+
+    def score_samples(self, X):
+        """
+        Opposite of the LoOP probability (higher is more normal).
+        """
+        return -1.0 * self._score_matrix(X)
+
+    def predict(self, X):
+        """
+        Predict inliers (1) and outliers (-1).
+
+        Labels are derived from ``decision_function(X) >= 0``.
+        """
+        decision = self.decision_function(X)
+        return np.where(decision >= 0, 1, -1)
+
+
+def _issparse_safe(x) -> bool:
+    try:
+        from scipy.sparse import issparse
+
+        return issparse(x)
+    except ImportError:
+        return False
 
 
 LoOP = LocalOutlierProbability
+
+__all__ = [
+    "LocalOutlierProbability",
+    "LoOP",
+    # Exceptions live in PyNomaly.exceptions; re-exported here so that
+    # ``from PyNomaly.loop import ClusterSizeError`` keeps working.
+    "PyNomalyError",
+    "ValidationError",
+    "ClusterSizeError",
+    "MissingValuesError",
+]
diff --git a/README.md b/README.md
index 7a5f7b6..5f4245d 100644
--- a/README.md
+++ b/README.md
@@ -60,7 +60,7 @@ LoOP scores are computed through a pipeline of steps that transform raw distance
 For full mathematical detail, see the [original paper](http://www.dbs.ifi.lmu.de/Publikationen/Papers/LoOP1649.pdf) by Kriegel et al. (2009). For the full API and additional documentation, see the [PyNomaly documentation](https://vc1492a.github.io/PyNomaly/).
 
 ## Dependencies
-- Python 3.8 - 3.14
+- Python 3.9 - 3.14
 - numpy >= 1.16.3
 - python-utils >= 2.3.0
 - (optional) numba >= 0.45.1
@@ -72,6 +72,9 @@ distance between observations, providing a reduction in computation time
 requirement and PyNomaly may still be used solely with numpy if desired
 (details below). When scipy is available, PyNomaly uses its optimized distance 
 computation and error function implementations for additional performance gains.
+Sparse matrices (`scipy.sparse`) are also accepted as input and converted to
+dense arrays internally; **scipy is required for sparse input** (`pip install scipy`
+or `pip install PyNomaly[sparse]`).
 
 ## Quick Start
 
@@ -394,11 +397,13 @@ neigh.fit(data)
 d, idx = neigh.kneighbors(data, return_distance=True)
 
 # Remove self-distances - you MUST do this to preserve the same results as intended by the definition of LoOP
-indices = np.delete(indices, 0, 1)
-distances = np.delete(distances, 0, 1)
+idx = np.delete(idx, 0, 1)
+d = np.delete(d, 0, 1)
 
 # Fit and return scores
-m = loop.LocalOutlierProbability(n_neighbors=n_neighbors+1).fit(distance_matrix=d, neighbor_matrix=idx)
+m = loop.LocalOutlierProbability(n_neighbors=n_neighbors).fit(
+    distance_matrix=d, neighbor_matrix=idx
+)
 scores = m.local_outlier_probabilities
 ```
 
diff --git a/changelog.md b/changelog.md
index 80f809a..ada9ae2 100644
--- a/changelog.md
+++ b/changelog.md
@@ -16,11 +16,24 @@ and adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
   `FutureWarning` and will be removed in a future version.
 - Calling `fit()` multiple times with different data is now supported 
   (re-fitting resets internal state).
+- `stream()` now returns a Python `float` instead of a NumPy array.
+- Exception classes (`PyNomalyError`, `ValidationError`, `ClusterSizeError`,
+  `MissingValuesError`) now live in the new `PyNomaly.exceptions` module.
+  They remain importable from `PyNomaly.loop` and from `PyNomaly` for
+  backward compatibility.
+- Passing data as a positional argument to the constructor
+  (e.g. `LocalOutlierProbability(data)`) no longer works: the first
+  positional parameter is now `extent`. Pass data to `fit()` instead
+  (e.g. `LocalOutlierProbability().fit(data)`). Keyword form
+  `LocalOutlierProbability(data=data)` still works with a `FutureWarning`.
 ### Added
 - `LoOP` class alias for `LocalOutlierProbability`, enabling 
   `from PyNomaly import LoOP` 
   ([Issue #5](https://github.com/vc1492a/PyNomaly/issues/5)).
 - `__version__` is now exported from the package (`from PyNomaly import __version__`).
+### Removed
+- Support for Python 3.8. PyNomaly 1.0.0 requires Python 3.9 or newer
+  (`python_requires='>=3.9'`).
 
 ## 0.4.0
 ### Added
diff --git a/docs/api.md b/docs/api.md
index a37d3a9..79b5c50 100644
--- a/docs/api.md
+++ b/docs/api.md
@@ -22,7 +22,7 @@ clf = loop.LocalOutlierProbability(
 
 | Parameter | Type | Default | Description |
 |---|---|---|---|
-| `data` | `np.ndarray` or `pd.DataFrame` | `None` | Input data as a 2D array with shape (n_observations, n_features). Mutually exclusive with `distance_matrix`. |
+| `data` | `np.ndarray`, `pd.DataFrame`, or `scipy.sparse` matrix | `None` | Input data as a 2D array with shape (n_observations, n_features). Sparse matrices are accepted and densified internally (requires scipy). Mutually exclusive with `distance_matrix`. |
 | `distance_matrix` | `np.ndarray` or `pd.DataFrame` | `None` | Pre-computed distance matrix with shape (n_observations, n_neighbors). Must be provided together with `neighbor_matrix`. Mutually exclusive with `data`. |
 | `neighbor_matrix` | `np.ndarray` or `pd.DataFrame` | `None` | Pre-computed neighbor index matrix with shape (n_observations, n_neighbors). Required when `distance_matrix` is provided. |
 | `extent` | `int` | `3` | Controls scoring sensitivity. Must be 1, 2, or 3. Corresponds to lambda times the standard deviation from the mean (1 = ~68%, 2 = ~95%, 3 = ~99.7%). |
@@ -47,7 +47,7 @@ clf.fit() -> LocalOutlierProbability
 
 Calculates the Local Outlier Probability for each observation in the input data.
 
-**Returns**: `self` -- the fitted model instance. Access scores via `clf.local_outlier_probabilities`.
+**Returns**: `self` -- the fitted model instance. Access scores via `clf.local_outlier_probabilities_`.
 
 **Raises**:
 
@@ -59,7 +59,7 @@ Calculates the Local Outlier Probability for each observation in the input data.
 #### `stream(x)`
 
 ```python
-clf.stream(x) -> np.ndarray
+clf.stream(x) -> float
 ```
 
 Calculates the Local Outlier Probability for an individual observation against the fitted model. Must be called after `fit()`.
@@ -70,7 +70,7 @@ Calculates the Local Outlier Probability for an individual observation against t
 |---|---|---|
 | `x` | `np.ndarray` | A single observation to score. When using raw data mode, this should be a 1D array with the same number of features as the training data. When using distance matrix mode, this should be a scalar distance value. |
 
-**Returns**: `np.ndarray` -- the Local Outlier Probability of the input observation (a value in [0, 1]).
+**Returns**: `float` -- the Local Outlier Probability of the input observation (a value in [0, 1]).
 
 !!! warning
     The stream approach does **not** support clustered data. If `cluster_labels` were provided during `fit()`, PyNomaly will automatically refit using a single cluster and issue a `UserWarning`.
@@ -83,22 +83,25 @@ Attributes available after calling `fit()`:
 
 | Attribute | Type | Description |
 |---|---|---|
-| `local_outlier_probabilities` | `np.ndarray` | Array of LoOP scores for each observation, with values in [0, 1]. |
-| `prob_distances` | `np.ndarray` | Probabilistic distances for each observation. |
-| `prob_distances_ev` | `np.ndarray` | Expected values of probabilistic distances for each observation's neighborhood. |
-| `norm_prob_local_outlier_factor` | `float` | Maximum normalized probabilistic local outlier factor across all observations. Used internally by `stream()`. |
-| `is_fit` | `bool` | Whether the model has been fit. |
+| `local_outlier_probabilities_` | `np.ndarray` | Array of LoOP scores for each observation, with values in [0, 1]. |
+| `prob_distances_` | `np.ndarray` | Probabilistic distances for each observation. |
+| `prob_distances_ev_` | `np.ndarray` | Expected values of probabilistic distances for each observation's neighborhood. |
+| `norm_prob_local_outlier_factor_` | `float` | Maximum normalized probabilistic local outlier factor across all observations. Used internally by `stream()`. |
+| `is_fit_` | `bool` | Whether the model has been fit. |
 | `n_neighbors` | `int` | The number of neighbors used (may differ from the value passed to the constructor if it was adjusted). |
 
 ---
 
 ## Exceptions
 
-All exceptions are importable from `PyNomaly.loop` or directly from `PyNomaly`:
+All exceptions are defined in `PyNomaly.exceptions`. They are also importable
+directly from `PyNomaly`, and from `PyNomaly.loop` for backward compatibility:
 
 ```python
-from PyNomaly import ClusterSizeError, MissingValuesError
+from PyNomaly.exceptions import PyNomalyError, ValidationError, ClusterSizeError, MissingValuesError
 # or
+from PyNomaly import ClusterSizeError, MissingValuesError
+# or (backward compatible)
 from PyNomaly.loop import PyNomalyError, ValidationError, ClusterSizeError, MissingValuesError
 ```
 
diff --git a/docs/changelog.md b/docs/changelog.md
index 06b2624..c217d79 100644
--- a/docs/changelog.md
+++ b/docs/changelog.md
@@ -17,11 +17,24 @@ and adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
   `FutureWarning` and will be removed in a future version.
 - Calling `fit()` multiple times with different data is now supported 
   (re-fitting resets internal state).
+- `stream()` now returns a Python `float` instead of a NumPy array.
+- Exception classes (`PyNomalyError`, `ValidationError`, `ClusterSizeError`,
+  `MissingValuesError`) now live in the new `PyNomaly.exceptions` module.
+  They remain importable from `PyNomaly.loop` and from `PyNomaly` for
+  backward compatibility.
+- Passing data as a positional argument to the constructor
+  (e.g. `LocalOutlierProbability(data)`) no longer works: the first
+  positional parameter is now `extent`. Pass data to `fit()` instead
+  (e.g. `LocalOutlierProbability().fit(data)`). Keyword form
+  `LocalOutlierProbability(data=data)` still works with a `FutureWarning`.
 ### Added
 - `LoOP` class alias for `LocalOutlierProbability`, enabling 
   `from PyNomaly import LoOP` 
   ([Issue #5](https://github.com/vc1492a/PyNomaly/issues/5)).
 - `__version__` is now exported from the package (`from PyNomaly import __version__`).
+### Removed
+- Support for Python 3.8. PyNomaly 1.0.0 requires Python 3.9 or newer
+  (`python_requires='>=3.9'`).
 
 ## 0.4.0
 ### Added
diff --git a/docs/examples.md b/docs/examples.md
index 827d6ff..672f989 100644
--- a/docs/examples.md
+++ b/docs/examples.md
@@ -30,9 +30,9 @@ Cluster the data using DBSCAN and generate two sets of scores. In both cases, we
 
 ```python
 db = DBSCAN(eps=0.9, min_samples=10).fit(iris)
-m = loop.LocalOutlierProbability(iris).fit()
+m = loop.LocalOutlierProbability().fit(iris)
 scores_noclust = m.local_outlier_probabilities
-m_clust = loop.LocalOutlierProbability(iris, cluster_labels=list(db.labels_)).fit()
+m_clust = loop.LocalOutlierProbability().fit(iris, cluster_labels=list(db.labels_))
 scores_clust = m_clust.local_outlier_probabilities
 ```
 
diff --git a/docs/getting-started.md b/docs/getting-started.md
index e400ec2..f98adf6 100644
--- a/docs/getting-started.md
+++ b/docs/getting-started.md
@@ -2,11 +2,14 @@
 
 ## Dependencies
 
-- Python 3.8 - 3.14
+- Python 3.9 - 3.14
 - numpy >= 1.16.3
 - python-utils >= 2.3.0
 - (optional) numba >= 0.45.1
-- (optional) scipy >= 1.3.0
+- (optional) scipy >= 1.3.0 — performance optimizations and **sparse matrix input**
+  (`scipy.sparse` CSR/CSC matrices are accepted and densified internally)
+- (optional) scikit-learn >= 1.0 — full estimator API / pipeline integration
+  (`pip install PyNomaly[sklearn]`)
 
 Numba just-in-time (JIT) compiles the function which calculates the Euclidean
 distance between observations, providing a reduction in computation time
@@ -17,6 +20,17 @@ When scipy is available, PyNomaly uses its optimized distance
 computation (`scipy.spatial.distance.cdist`) and error function (`scipy.special.erf`)
 implementations for additional performance gains.
 
+Sparse matrices from `scipy.sparse` (e.g. CSR, CSC) can be passed to `fit()`,
+`predict()`, `decision_function()`, and `stream()`. They are converted to dense
+NumPy arrays before distance computation. **scipy must be installed** to use
+sparse input:
+
+```shell
+pip install scipy
+# or
+pip install PyNomaly[sparse]
+```
+
 ## Installation
 
 Install from the Python Package Index:
diff --git a/docs/user-guide.md b/docs/user-guide.md
index 496c1cf..7da1dc4 100644
--- a/docs/user-guide.md
+++ b/docs/user-guide.md
@@ -54,6 +54,39 @@ This provides **2-3x speedups** on multi-core machines (benchmarked on 8 cores).
 
 Numba must be installed to use JIT compilation. PyNomaly has been tested with Numba versions 0.45.1 through 0.65.1.
 
+## Sparse Matrix Input
+
+LoOP accepts `scipy.sparse` matrices (CSR, CSC, etc.) wherever raw feature data
+is expected — `fit(X)`, `predict(X)`, `decision_function(X)`, and `stream(x)`.
+Sparse inputs are **densified internally** before distance computation; there is
+no sparse-native distance kernel yet.
+
+**scipy is a soft dependency for sparse input.** If you pass a sparse matrix
+without scipy installed, PyNomaly raises a clear `ImportError` with install
+instructions. Install scipy with:
+
+```shell
+pip install scipy
+# or
+pip install PyNomaly[sparse]
+```
+
+Example:
+
+```python
+from scipy.sparse import csr_matrix
+from PyNomaly import LoOP
+
+X = csr_matrix([[1, 0, 2], [0, 3, 0], [4, 0, 5]], dtype=float)
+clf = LoOP(n_neighbors=2).fit(X)
+scores = clf.local_outlier_probabilities_
+```
+
+!!! note
+    Densifying very large sparse matrices can use significant memory. For
+    extremely high-dimensional sparse data, consider dimensionality reduction
+    or sampling before fitting LoOP.
+
 ## Progress Bars
 
 You may choose to print progress bars _with or without_ the use of Numba by passing `progress_bar=True` to `LocalOutlierProbability()`:
@@ -119,10 +152,10 @@ neigh.fit(data)
 d, idx = neigh.kneighbors(data, return_distance=True)
 
 # Remove self-distances
-indices = np.delete(indices, 0, 1)
-distances = np.delete(distances, 0, 1)
+idx = np.delete(idx, 0, 1)
+d = np.delete(d, 0, 1)
 
-m = loop.LocalOutlierProbability(n_neighbors=n_neighbors+1).fit(
+m = loop.LocalOutlierProbability(n_neighbors=n_neighbors).fit(
     distance_matrix=d, neighbor_matrix=idx
 )
 scores = m.local_outlier_probabilities
diff --git a/examples/iris.py b/examples/iris.py
index 24b71d6..d095778 100644
--- a/examples/iris.py
+++ b/examples/iris.py
@@ -11,9 +11,9 @@
 
 
 db = DBSCAN(eps=0.9, min_samples=10).fit(iris)
-m = loop.LocalOutlierProbability(iris).fit()
+m = loop.LocalOutlierProbability().fit(iris)
 scores_noclust = m.local_outlier_probabilities
-m_clust = loop.LocalOutlierProbability(iris, cluster_labels=list(db.labels_)).fit()
+m_clust = loop.LocalOutlierProbability().fit(iris, cluster_labels=list(db.labels_))
 scores_clust = m_clust.local_outlier_probabilities
 
 
diff --git a/examples/iris_dist_grid.py b/examples/iris_dist_grid.py
index 78092d2..1257aba 100644
--- a/examples/iris_dist_grid.py
+++ b/examples/iris_dist_grid.py
@@ -29,8 +29,9 @@
     neigh.fit(iris)
     d, idx = neigh.kneighbors(iris, return_distance=True)
 
-    m = loop.LocalOutlierProbability(distance_matrix=d,
-                                     neighbor_matrix=idx).fit()
+    m = loop.LocalOutlierProbability().fit(
+        distance_matrix=d, neighbor_matrix=idx
+    )
     iris['scores'] = m.local_outlier_probabilities
 
     ax = fig.add_subplot(3, 3, i, projection='3d')
diff --git a/examples/multiple_gaussian_2d.py b/examples/multiple_gaussian_2d.py
index 56fc7a0..e50ba28 100644
--- a/examples/multiple_gaussian_2d.py
+++ b/examples/multiple_gaussian_2d.py
@@ -8,7 +8,7 @@
 print(df)
 
 # fit LoOP according to the original settings outlined in the paper #
-m = loop.LocalOutlierProbability(df[['x', 'y']], n_neighbors=20, extent=3).fit()
+m = loop.LocalOutlierProbability(n_neighbors=20, extent=3).fit(df[['x', 'y']])
 scores = m.local_outlier_probabilities
 print(scores)
 
@@ -25,4 +25,3 @@
         plt.text(df['x'].loc[i] * (1 + 0.01), df['y'].loc[i] * (1 + 0.01), round(scores[i], 2), fontsize=8)
 
 plt.show()
-
diff --git a/examples/numba_speed_diff.py b/examples/numba_speed_diff.py
index 25260ee..0384a9e 100644
--- a/examples/numba_speed_diff.py
+++ b/examples/numba_speed_diff.py
@@ -10,11 +10,10 @@ def main():
     # Vectorized NumPy (default, no Numba)
     t1 = time.time()
     scores_numpy = loop.LocalOutlierProbability(
-        data,
         n_neighbors=3,
         use_numba=False,
         progress_bar=True
-    ).fit().local_outlier_probabilities
+    ).fit(data).local_outlier_probabilities
     t2 = time.time()
     seconds_no_numba = t2 - t1
     print("\nComputation took " + str(seconds_no_numba) + " seconds without Numba JIT.")
@@ -22,11 +21,10 @@ def main():
     # Numba JIT (sequential)
     t3 = time.time()
     scores_numba = loop.LocalOutlierProbability(
-        data,
         n_neighbors=3,
         use_numba=True,
         progress_bar=True
-    ).fit().local_outlier_probabilities
+    ).fit(data).local_outlier_probabilities
     t4 = time.time()
     seconds_numba = t4 - t3
     print("\nComputation took " + str(seconds_numba) + " seconds with Numba JIT (sequential).")
@@ -34,12 +32,11 @@ def main():
     # Numba JIT (parallel, n_jobs=-1)
     t5 = time.time()
     scores_numba_par = loop.LocalOutlierProbability(
-        data,
         n_neighbors=3,
         use_numba=True,
         n_jobs=-1,
         progress_bar=True
-    ).fit().local_outlier_probabilities
+    ).fit(data).local_outlier_probabilities
     t6 = time.time()
     seconds_numba_par = t6 - t5
     print("\nComputation took " + str(seconds_numba_par) + " seconds with Numba JIT (parallel, n_jobs=-1).")
diff --git a/examples/parallel_benchmark.py b/examples/parallel_benchmark.py
index 432a165..638631f 100644
--- a/examples/parallel_benchmark.py
+++ b/examples/parallel_benchmark.py
@@ -42,13 +42,12 @@ def generate_clustered_data(n_per_cluster, n_clusters, n_features=N_FEATURES, se
 
 def time_fit(data, cluster_labels, n_neighbors, use_numba, n_jobs):
     clf = loop.LocalOutlierProbability(
-        data, n_neighbors=n_neighbors,
-        cluster_labels=cluster_labels,
+        n_neighbors=n_neighbors,
         use_numba=use_numba,
-        n_jobs=n_jobs
+        n_jobs=n_jobs,
     )
     t0 = time.perf_counter()
-    clf.fit()
+    clf.fit(data, cluster_labels=cluster_labels)
     return time.perf_counter() - t0
 
 
diff --git a/examples/stream.py b/examples/stream.py
index c00d90f..3ecfb3e 100644
--- a/examples/stream.py
+++ b/examples/stream.py
@@ -12,12 +12,12 @@
 iris_train = iris.iloc[:, 0:4].head(120)
 iris_test = iris.iloc[:, 0:4].tail(30)
 
-m = loop.LocalOutlierProbability(iris).fit()
+m = loop.LocalOutlierProbability().fit(iris)
 scores_noclust = m.local_outlier_probabilities
 iris['scores'] = scores_noclust
 
-m_train = loop.LocalOutlierProbability(iris_train, n_neighbors=10)
-m_train.fit()
+m_train = loop.LocalOutlierProbability(n_neighbors=10)
+m_train.fit(iris_train)
 iris_train_scores = m_train.local_outlier_probabilities
 
 iris_test_scores = []
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..61e1c76
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,12 @@
+[tool:pytest]
+norecursedirs = .* build dist venv worktrees *.egg-info
+
+[flake8]
+exclude =
+    .git,
+    __pycache__,
+    build,
+    dist,
+    venv,
+    worktrees,
+    *.egg-info
diff --git a/setup.py b/setup.py
index e701fc2..8895c72 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,26 @@
     download_url='https://github.com/vc1492a/PyNomaly/archive/1.0.0.tar.gz',
     keywords=['outlier', 'anomaly', 'detection', 'machine', 'learning',
               'probability'],
-    classifiers=[],
+    classifiers=[
+        'Intended Audience :: Science/Research',
+        'License :: OSI Approved :: Apache Software License',
+        'Operating System :: OS Independent',
+        'Programming Language :: Python :: 3',
+        'Programming Language :: Python :: 3.9',
+        'Programming Language :: Python :: 3.10',
+        'Programming Language :: Python :: 3.11',
+        'Programming Language :: Python :: 3.12',
+        'Programming Language :: Python :: 3.13',
+        'Programming Language :: Python :: 3.14',
+        'Topic :: Scientific/Engineering',
+    ],
+    python_requires='>=3.9',
     license='Apache License, Version 2.0',
-    install_requires=['numpy', 'python-utils']
+    install_requires=['numpy', 'python-utils'],
+    extras_require={
+        'sklearn': ['scikit-learn>=1.0', 'scipy>=1.3.0'],
+        'sparse': ['scipy>=1.3.0'],
+        'numba': ['numba'],
+        'all': ['scikit-learn>=1.0', 'numba', 'scipy>=1.3.0'],
+    },
 )
diff --git a/tests/test_loop.py b/tests/test_loop.py
index 0ce1be3..48522d2 100644
--- a/tests/test_loop.py
+++ b/tests/test_loop.py
@@ -1,23 +1,23 @@
 # Authors: Valentino Constantinou <vc@valentino.io>
 # License: Apache 2.0
 
-from PyNomaly import loop
-from PyNomaly.loop import ClusterSizeError, MissingValuesError
-
 import logging
-from typing import Tuple
+import re
+import sys
 import warnings
+from typing import Tuple
+
 import numpy as np
-from numpy.testing import assert_array_equal, assert_array_almost_equal
 import pandas as pd
 import pytest
+from numpy.testing import assert_array_almost_equal, assert_array_equal
 from sklearn.datasets import load_iris
 from sklearn.metrics import roc_auc_score
 from sklearn.neighbors import NearestNeighbors
 from sklearn.utils import check_random_state
 from sklearn.utils.estimator_checks import check_estimator
 
-import sys
+from PyNomaly import ClusterSizeError, MissingValuesError, loop
 
 logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
 
@@ -162,8 +162,24 @@ def test_sklearn_api_compliance():
     """
     Check scikit-learn compatibility as referenced in the developer guide:
     https://scikit-learn.org/stable/developers/develop.html#rolling-your-own-estimator
-    """
-    check_estimator(loop.LoOP())
+
+    Requires the optional sklearn extra: ``pip install PyNomaly[sklearn]``.
+
+    ``check_fit1d`` is an expected failure: LoOP historically accepts 1-d
+    inputs (e.g. univariate series) by reshaping to ``(-1, 1)`` rather than
+    raising.
+    """
+    pytest.importorskip("sklearn")
+    check_estimator(
+        loop.LoOP(),
+        expected_failed_checks={
+            "check_fit1d": (
+                "LoOP accepts 1-d inputs and reshapes them to (-1, 1) for "
+                "univariate / time-series use cases."
+            ),
+        },
+    )
+
 
 def test_loop(X_n8) -> None:
     """
@@ -175,7 +191,7 @@ def test_loop(X_n8) -> None:
     """
     # Test LocalOutlierProbability:
     clf = loop.LocalOutlierProbability(n_neighbors=5, use_numba=NUMBA)
-    score = clf.fit(X_n8).local_outlier_probabilities
+    score = clf.fit(X_n8).local_outlier_probabilities_
     share_outlier = 2.0 / 8.0
     predictions = [-1 if s > share_outlier else 1 for s in score]
     assert_array_equal(predictions, 6 * [1] + 2 * [-1])
@@ -188,7 +204,7 @@ def test_loop(X_n8) -> None:
 
     # Test LocalOutlierProbability:
     clf = loop.LocalOutlierProbability(n_neighbors=5, use_numba=NUMBA)
-    score = clf.fit(X_df).local_outlier_probabilities
+    score = clf.fit(X_df).local_outlier_probabilities_
     share_outlier = 2.0 / 8.0
     predictions = [-1 if s > share_outlier else 1 for s in score]
     assert_array_equal(predictions, 6 * [1] + 2 * [-1])
@@ -205,7 +221,7 @@ def test_regression(X_n20_scores) -> None:
     """
     input_data, expected_scores = X_n20_scores
     clf = loop.LocalOutlierProbability().fit(input_data)
-    scores = clf.local_outlier_probabilities
+    scores = clf.local_outlier_probabilities_
     assert_array_almost_equal(scores, expected_scores, 6)
 
 
@@ -230,7 +246,7 @@ def test_loop_performance(X_n120) -> None:
     )
 
     # predict scores (the lower, the more normal)
-    score = clf.fit(X_test).local_outlier_probabilities
+    score = clf.fit(X_test).local_outlier_probabilities_
     share_outlier = X_outliers.shape[0] / X_test.shape[0]
     X_pred = [-1 if s > share_outlier else 1 for s in score]
 
@@ -245,18 +261,13 @@ def test_input_nodata(X_n140_outliers) -> None:
     :param X_n140_outliers: A pytest Fixture that generates 140 observations.
     :return: None
     """
-    with pytest.warns(UserWarning) as record:
+    with pytest.warns(UserWarning, match="Data or a distance matrix must be provided.") as record:
         # attempt to fit loop without data or a distance matrix
         clf = loop.LocalOutlierProbability(
             n_neighbors=X_n140_outliers.shape[0] - 1, use_numba=NUMBA
         )
         clf.fit()
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert record[0].message.args[0] == "Data or a distance matrix must be provided."
-
 
 def test_input_incorrect_type(X_n140_outliers) -> None:
     """
@@ -265,21 +276,15 @@ def test_input_incorrect_type(X_n140_outliers) -> None:
     :param X_n140_outliers: A pytest Fixture that generates 140 observations.
     :return: None
     """
-    with pytest.warns(UserWarning) as record:
+    
+    match_str = re.escape("Argument 'n_neighbors' is not of type (<class 'int'>, <class 'numpy.integer'>).")
+    with pytest.warns(UserWarning, match=match_str):
         # attempt to fit loop with a string input for n_neighbors
-        loop.LocalOutlierProbability(
+        clf = loop.LocalOutlierProbability(
             n_neighbors=str(X_n140_outliers.shape[0] - 1),
             use_numba=NUMBA,
         )
-
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0]
-        == "Argument 'n_neighbors' is not of type (<class 'int'>, "
-        "<class 'numpy.integer'>)."
-    )
+        clf.fit()   
 
 
 def test_input_neighbor_zero(X_n120) -> None:
@@ -289,19 +294,12 @@ def test_input_neighbor_zero(X_n120) -> None:
     :param X_n120: A pytest Fixture that generates 120 observations.
     :return: None
     """
-    clf = loop.LocalOutlierProbability(n_neighbors=0, use_numba=NUMBA)
 
-    with pytest.warns(UserWarning) as record:
+    clf = loop.LocalOutlierProbability(n_neighbors=0, use_numba=NUMBA)
+    with pytest.warns(UserWarning, match="n_neighbors must be greater than 0. Fit with 10 instead."):
         # attempt to fit loop with a 0 neighbor count
         clf.fit(X_n120)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0]
-        == "n_neighbors must be greater than 0. Fit with 10 instead."
-    )
 
 
 def test_input_distonly(X_n120) -> None:
@@ -316,19 +314,12 @@ def test_input_distonly(X_n120) -> None:
     neigh.fit(X_n120)
     d, idx = neigh.kneighbors(X_n120, n_neighbors=10, return_distance=True)
 
-    with pytest.warns(UserWarning) as record:
+    clf = loop.LocalOutlierProbability(use_numba=NUMBA)
+    match_str = re.escape("A neighbor index matrix and distance matrix must both be provided when not using raw input data.")
+    with pytest.warns(UserWarning, match=match_str):
         # attempt to fit loop with a distance matrix and no neighbor matrix
-        clf = loop.LocalOutlierProbability(use_numba=NUMBA)
         clf.fit(distance_matrix=d)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0]
-        == "A neighbor index matrix and distance matrix must both "
-        "be provided when not using raw input data."
-    )
 
 
 def test_input_neighboronly(X_n120) -> None:
@@ -343,16 +334,12 @@ def test_input_neighboronly(X_n120) -> None:
     neigh.fit(X_n120)
     d, idx = neigh.kneighbors(X_n120, n_neighbors=10, return_distance=True)
 
-    with pytest.warns(UserWarning) as record:
+    clf = loop.LocalOutlierProbability(use_numba=NUMBA)
+    
+    with pytest.warns(UserWarning, match="Data or a distance matrix must be provided.") as record:
         # attempt to fit loop with a neighbor matrix and no distance matrix
-        clf = loop.LocalOutlierProbability(use_numba=NUMBA)
         clf.fit(neighbor_matrix=idx)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert record[0].message.args[0] == "Data or a distance matrix must be provided."
-
 
 def test_input_too_many(X_n120) -> None:
     """
@@ -365,20 +352,12 @@ def test_input_too_many(X_n120) -> None:
     neigh = NearestNeighbors(metric="euclidean")
     neigh.fit(X_n120)
     d, idx = neigh.kneighbors(X_n120, n_neighbors=10, return_distance=True)
-
-    with pytest.warns(UserWarning) as record:
+    clf = loop.LocalOutlierProbability(use_numba=NUMBA)
+    match_str = re.escape("Only one of the following may be provided: data or a distance matrix (not both).")
+    with pytest.warns(UserWarning, match=match_str):
         # attempt to fit loop with data and a distance matrix
-        clf = loop.LocalOutlierProbability(use_numba=NUMBA)
         clf.fit(X_n120, distance_matrix=d, neighbor_matrix=idx)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0]
-        == "Only one of the following may be provided: data or a "
-        "distance matrix (not both)."
-    )
 
 
 def test_distance_neighbor_shape_mismatch(X_n120) -> None:
@@ -397,20 +376,13 @@ def test_distance_neighbor_shape_mismatch(X_n120) -> None:
     neigh_2 = NearestNeighbors(metric="euclidean")
     neigh_2.fit(X_n120)
     d_2, idx_2 = neigh.kneighbors(X_n120, n_neighbors=5, return_distance=True)
+    clf = loop.LocalOutlierProbability(n_neighbors=5, use_numba=NUMBA)
 
-    with pytest.warns(UserWarning) as record:
+    with pytest.warns(UserWarning, match="The shape of the distance and neighbor "
+        "index matrices must match.") as record:
         # attempt to fit loop with a mismatch in shapes
-        clf = loop.LocalOutlierProbability(n_neighbors=5, use_numba=NUMBA)
         clf.fit(distance_matrix=d, neighbor_matrix=idx_2)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0] == "The shape of the distance and neighbor "
-        "index matrices must match."
-    )
-
 
 def test_input_neighbor_mismatch(X_n120) -> None:
     """
@@ -423,21 +395,15 @@ def test_input_neighbor_mismatch(X_n120) -> None:
     neigh = NearestNeighbors(metric="euclidean")
     neigh.fit(X_n120)
     d, idx = neigh.kneighbors(X_n120, n_neighbors=5, return_distance=True)
+    clf = loop.LocalOutlierProbability(n_neighbors=10, use_numba=NUMBA)
 
-    with pytest.warns(UserWarning) as record:
+    with pytest.warns(UserWarning, match=("The shape of the distance or "
+        "neighbor index matrix does not "
+        "match the number of neighbors "
+        "specified.")) as record:
         # attempt to fit loop with a neighbor size mismatch
-        clf = loop.LocalOutlierProbability(n_neighbors=10, use_numba=NUMBA)
         clf.fit(distance_matrix=d, neighbor_matrix=idx)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0] == "The shape of the distance or "
-        "neighbor index matrix does not "
-        "match the number of neighbors "
-        "specified."
-    )
 
 
 def test_loop_dist_matrix(X_n120) -> None:
@@ -455,8 +421,8 @@ def test_loop_dist_matrix(X_n120) -> None:
     # fit loop using data and distance matrix
     clf1 = loop.LocalOutlierProbability(use_numba=NUMBA)
     clf2 = loop.LocalOutlierProbability(use_numba=NUMBA)
-    scores1 = clf1.fit(X_n120).local_outlier_probabilities
-    scores2 = clf2.fit(distance_matrix=d, neighbor_matrix=idx).local_outlier_probabilities
+    scores1 = clf1.fit(X_n120).local_outlier_probabilities_
+    scores2 = clf2.fit(distance_matrix=d, neighbor_matrix=idx).local_outlier_probabilities_
 
     # compare the agreement between the results
     assert np.abs(scores2 - scores1).all() <= 0.1
@@ -476,9 +442,9 @@ def test_lambda_values(X_n140_outliers) -> None:
     clf3 = loop.LocalOutlierProbability(extent=3, use_numba=NUMBA)
 
     # predict scores (the lower, the more normal)
-    score1 = clf1.fit(X_n140_outliers).local_outlier_probabilities
-    score2 = clf2.fit(X_n140_outliers).local_outlier_probabilities
-    score3 = clf3.fit(X_n140_outliers).local_outlier_probabilities
+    score1 = clf1.fit(X_n140_outliers).local_outlier_probabilities_
+    score2 = clf2.fit(X_n140_outliers).local_outlier_probabilities_
+    score3 = clf3.fit(X_n140_outliers).local_outlier_probabilities_
 
     # Get the mean of all the scores
     score_mean1 = np.mean(score1)
@@ -505,28 +471,195 @@ def test_parameters(X_n120) -> None:
     assert hasattr(clf, "n_neighbors") and clf.n_neighbors is not None
     assert hasattr(clf, "extent") and clf.extent is not None
     assert hasattr(clf, "cluster_labels") and clf._cluster_labels() is not None
-    assert hasattr(clf, "prob_distances") and clf.prob_distances is not None
-    assert hasattr(clf, "prob_distances_ev") and clf.prob_distances_ev is not None
+    assert hasattr(clf, "prob_distances") and clf.prob_distances_ is not None
+    assert hasattr(clf, "prob_distances_ev") and clf.prob_distances_ev_ is not None
     assert (
         hasattr(clf, "norm_prob_local_outlier_factor")
-        and clf.norm_prob_local_outlier_factor is not None
+        and clf.norm_prob_local_outlier_factor_ is not None
     )
     assert (
         hasattr(clf, "local_outlier_probabilities")
-        and clf.local_outlier_probabilities is not None
+        and clf.local_outlier_probabilities_ is not None
     )
 
 
+def test_public_api_contract(X_n120) -> None:
+    """
+    Guards against unintentional breaking changes to the public API of
+    LocalOutlierProbability. Exercises every non-underscored constructor
+    keyword argument, every fit() keyword argument, public method, and
+    public attribute so that a rename or removal of any of them causes this
+    test to fail. Additions to the API are unaffected.
+
+    Written against the 1.0.0 scikit-learn-style API: the constructor takes
+    configuration only, data goes to ``fit(X, y=None, ...)``, and fitted
+    state lives on ``_``-suffixed attributes (with read-only legacy aliases).
+    :param X_n120: A pytest Fixture that generates 120 observations.
+    :return: None
+    """
+    # exercise every public constructor keyword argument by name; a rename
+    # or removal raises a TypeError here. Data parameters are no longer
+    # constructor arguments in 1.0.0 -- they are exercised via fit() below.
+    clf = loop.LocalOutlierProbability(
+        extent=3,
+        n_neighbors=10,
+        use_numba=NUMBA,
+        n_jobs=1,
+        progress_bar=False,
+    )
+
+    # sklearn contract: get_params() exposes exactly the constructor's
+    # parameters -- the configuration parameters plus the deprecated data
+    # keywords, which stay constructor arguments (with a FutureWarning)
+    # until they are removed in a future version
+    assert set(clf.get_params()) == {
+        "extent",
+        "n_neighbors",
+        "use_numba",
+        "n_jobs",
+        "progress_bar",
+        "data",
+        "distance_matrix",
+        "neighbor_matrix",
+        "cluster_labels",
+    }
+
+    # exercise every public method by name, and every fit() keyword
+    # argument by name (a rename or removal raises a TypeError here)
+    assert hasattr(clf, "fit") and callable(clf.fit)
+    fitted = clf.fit(
+        X_n120,
+        y=None,
+        distance_matrix=None,
+        neighbor_matrix=None,
+        cluster_labels=None,
+    )
+    assert fitted is clf
+    for method_name in (
+        "fit",
+        "fit_predict",
+        "predict",
+        "decision_function",
+        "score_samples",
+        "stream",
+        "get_params",
+        "set_params",
+    ):
+        assert hasattr(fitted, method_name) and callable(
+            getattr(fitted, method_name)
+        ), f"missing public method: {method_name}"
+    fitted.stream(X_n120[0])
+    fitted.predict(X_n120)
+    fitted.decision_function(X_n120)
+    fitted.score_samples(X_n120)
+
+    # every public (non-underscored) attribute expected to exist post-fit:
+    # configuration parameters, the deprecated constructor data inputs,
+    # the `_`-suffixed fitted state, and the read-only legacy aliases
+    expected_public_attributes = [
+        # configuration parameters (constructor)
+        "extent",
+        "n_neighbors",
+        "use_numba",
+        "n_jobs",
+        "progress_bar",
+        # deprecated constructor data inputs (always present, None here)
+        "data",
+        "distance_matrix",
+        "neighbor_matrix",
+        "cluster_labels",
+        # fitted state (scikit-learn convention)
+        "data_",
+        "distance_matrix_",
+        "neighbor_matrix_",
+        "cluster_labels_",
+        "n_neighbors_",
+        "n_features_in_",
+        "offset_",
+        "points_vector_",
+        "prob_distances_",
+        "prob_distances_ev_",
+        "norm_prob_local_outlier_factor_",
+        "local_outlier_probabilities_",
+        "is_fit_",
+        # read-only legacy aliases for fitted output attributes
+        "points_vector",
+        "prob_distances",
+        "prob_distances_ev",
+        "norm_prob_local_outlier_factor",
+        "local_outlier_probabilities",
+        "is_fit",
+    ]
+    for attr_name in expected_public_attributes:
+        assert hasattr(fitted, attr_name), f"missing public attribute: {attr_name}"
+
+    # the legacy aliases must read the canonical fitted attributes
+    assert fitted.local_outlier_probabilities is fitted.local_outlier_probabilities_
+    assert fitted.is_fit is fitted.is_fit_
+
+    # exercise every public constructor argument positionally, in the
+    # documented 1.0.0 order, using a distinct value per slot. Swapping two
+    # parameters of the same type (e.g. extent/n_neighbors/n_jobs, or
+    # use_numba/progress_bar) doesn't raise an error -- it silently
+    # assigns the wrong value -- so each assigned attribute is checked
+    # against the value passed at its expected position.
+    positional_clf = loop.LocalOutlierProbability(2, 15, False, 1, True)
+    assert positional_clf.extent == 2
+    assert positional_clf.n_neighbors == 15
+    assert positional_clf.use_numba is False
+    assert positional_clf.n_jobs == 1
+    assert positional_clf.progress_bar is True
+
+    # data is positional in fit(), sklearn-style: fit(X, y)
+    positional_fit = loop.LocalOutlierProbability(n_neighbors=10).fit(
+        X_n120, None
+    )
+    assert positional_fit.data_ is X_n120
+    assert positional_fit.is_fit_ is True
+
+    # the `LoOP` alias is part of the public API
+    from PyNomaly import LoOP
+
+    assert LoOP is loop.LocalOutlierProbability
+    assert loop.LoOP is loop.LocalOutlierProbability
+
+    # the exception hierarchy is part of the public API: customers import
+    # these names directly (as this test file does at the top) and may
+    # catch on the base classes, so both the import paths and the
+    # hierarchy itself must remain stable. The canonical home is
+    # PyNomaly.exceptions; PyNomaly.loop and PyNomaly re-export them.
+    from PyNomaly import exceptions
+    from PyNomaly.loop import (
+        PyNomalyError,
+        ValidationError,
+        ClusterSizeError as ImportedClusterSizeError,
+        MissingValuesError as ImportedMissingValuesError,
+    )
+
+    assert PyNomalyError is exceptions.PyNomalyError
+    assert ValidationError is exceptions.ValidationError
+    assert ImportedClusterSizeError is exceptions.ClusterSizeError
+    assert ImportedMissingValuesError is exceptions.MissingValuesError
+    assert issubclass(ImportedClusterSizeError, ValidationError)
+    assert issubclass(ImportedMissingValuesError, ValidationError)
+    assert issubclass(ValidationError, PyNomalyError)
+    assert issubclass(PyNomalyError, Exception)
+
+
 def test_n_neighbors() -> None:
     """
     Tests the functionality of providing a large number of neighbors that
     is greater than the number of observations (software defaults to the
     data input size and provides a UserWarning).
+
+    The public ``n_neighbors`` parameter is left unchanged (sklearn
+    immutability); the effective value used during fit is ``n_neighbors_``.
     :return: None
     """
     X = iris.data
     clf = loop.LocalOutlierProbability(n_neighbors=500, use_numba=NUMBA).fit(X)
-    assert clf.n_neighbors == X.shape[0] - 1
+    assert clf.n_neighbors == 500
+    assert clf.n_neighbors_ == X.shape[0] - 1
 
     clf = loop.LocalOutlierProbability(n_neighbors=500, use_numba=NUMBA)
 
@@ -536,7 +669,8 @@ def test_n_neighbors() -> None:
     # check that only one warning was raised
     assert len(record) == 1
 
-    assert clf.n_neighbors == X.shape[0] - 1
+    assert clf.n_neighbors == 500
+    assert clf.n_neighbors_ == X.shape[0] - 1
 
 
 def test_extent() -> None:
@@ -547,13 +681,10 @@ def test_extent() -> None:
     """
     X = np.array([[1, 1], [1, 0]])
     clf = loop.LocalOutlierProbability(n_neighbors=2, extent=4, use_numba=NUMBA)
-
-    with pytest.warns(UserWarning) as record:
+    match_str = re.escape("extent parameter (lambda) must be 1, 2, or 3.")
+    with pytest.warns(UserWarning, match=match_str):
         clf.fit(X)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-
 
 def test_data_format() -> None:
     """
@@ -565,14 +696,9 @@ def test_data_format() -> None:
     X = [1.3, 1.1, 0.9, 1.4, 1.5, 3.2]
     clf = loop.LocalOutlierProbability(n_neighbors=3, use_numba=NUMBA)
 
-    with pytest.warns(UserWarning) as record:
+    with pytest.warns(UserWarning, match="Provided data or distance matrix must be in ndarray") as record:
         clf.fit(X)
 
-    assert any(
-        "Provided data or distance matrix must be in ndarray" in r.message.args[0]
-        for r in record
-    )
-
 
 def test_missing_values() -> None:
     """
@@ -589,7 +715,7 @@ def test_missing_values() -> None:
     # check that the message matches
     assert (
         str(record.value)
-        == "Method does not support missing values in input data."
+        == "Input contains NaN, infinity or a value too large for dtype('float64')."
     )
 
 
@@ -634,17 +760,11 @@ def test_stream_fit(X_n140_outliers) -> None:
     X_test = X_n140_outliers[139]
     clf = loop.LocalOutlierProbability(use_numba=NUMBA)
     clf.data = X_train
+    match_str = re.escape("Must fit on historical data by calling fit() prior to calling stream(x).")
 
-    with pytest.warns(UserWarning) as record:
+    with pytest.warns(UserWarning, match=match_str):
         clf.stream(X_test)
 
-    # check that the message matches
-    messages = [i.message.args[0] for i in record]
-    assert (
-        "Must fit on historical data by calling fit() prior to "
-        "calling stream(x)." in messages
-    )
-
 
 def test_stream_distance(X_n140_outliers) -> None:
     """
@@ -705,16 +825,10 @@ def test_stream_cluster(X_n140_outliers) -> None:
         use_numba=NUMBA
     ).fit(X_train, cluster_labels=cluster_labels)
 
-    with pytest.warns(UserWarning) as record:
+    with pytest.warns(UserWarning, match=("Stream approach does not support clustered data. "
+        "Automatically refit using single cluster of points.")) as record:
         clf.stream(X_test)
 
-    # check that only one warning was raised
-    assert len(record) == 1
-    # check that the message matches
-    assert (
-        record[0].message.args[0] == "Stream approach does not support clustered data. "
-        "Automatically refit using single cluster of points."
-    )
 
 
 def test_stream_performance(X_n140_outliers) -> None:
@@ -731,11 +845,11 @@ def test_stream_performance(X_n140_outliers) -> None:
 
     # Fit the models in standard and stream form
     m = loop.LocalOutlierProbability(use_numba=NUMBA).fit(X_n140_outliers)
-    scores_noclust = m.local_outlier_probabilities
+    scores_noclust = m.local_outlier_probabilities_
 
     m_train = loop.LocalOutlierProbability(use_numba=NUMBA)
     m_train.fit(X_train)
-    X_train_scores = m_train.local_outlier_probabilities
+    X_train_scores = m_train.local_outlier_probabilities_
 
     X_test_scores = []
     for idx in range(X_test.shape[0]):
@@ -775,18 +889,18 @@ def test_data_flipping() -> None:
     clus = np.append(np.ones(n), 2 * np.ones(n)).tolist()
     model = loop.LocalOutlierProbability(n_neighbors=5)
     fit = model.fit(data, cluster_labels=clus)
-    res = fit.local_outlier_probabilities
+    res = fit.local_outlier_probabilities_
 
     data_flipped = np.flipud(data)
     clus_flipped = np.flipud(clus).tolist()
     model2 = loop.LocalOutlierProbability(n_neighbors=5)
     fit2 = model2.fit(data_flipped, cluster_labels=clus_flipped)
-    res2 = np.flipud(fit2.local_outlier_probabilities)
+    res2 = np.flipud(fit2.local_outlier_probabilities_)
 
     assert_array_almost_equal(res, res2, decimal=6)
     assert_array_almost_equal(
-        fit.norm_prob_local_outlier_factor,
-        fit2.norm_prob_local_outlier_factor,
+        fit.norm_prob_local_outlier_factor_,
+        fit2.norm_prob_local_outlier_factor_,
         decimal=6,
     )
 
@@ -811,8 +925,8 @@ def test_distance_matrix_consistency(X_n120) -> None:
     clf_dist = loop.LocalOutlierProbability(n_neighbors=10)
 
     # Attempt to retrieve scores and check types
-    scores_data = clf_data.fit(X_n120).local_outlier_probabilities
-    scores_dist = clf_dist.fit(distance_matrix=distances, neighbor_matrix=indices).local_outlier_probabilities
+    scores_data = clf_data.fit(X_n120).local_outlier_probabilities_
+    scores_dist = clf_dist.fit(distance_matrix=distances, neighbor_matrix=indices).local_outlier_probabilities_
 
     # Debugging prints to investigate types and contents
     print("Type of scores_data:", type(scores_data))
@@ -842,7 +956,7 @@ def test_vectorized_1d_data() -> None:
     """
     X = np.array([1.0, 2.0, 3.0, 10.0, 11.0, 12.0, 50.0])
     clf = loop.LocalOutlierProbability(n_neighbors=3)
-    scores = clf.fit(X).local_outlier_probabilities
+    scores = clf.fit(X).local_outlier_probabilities_
     assert scores is not None
     assert len(scores) == len(X)
     assert scores[-1] > 0
@@ -857,7 +971,7 @@ def test_n_jobs_without_numba_warns(X_n120) -> None:
         clf = loop.LocalOutlierProbability(
             n_neighbors=10, n_jobs=2
         )
-        scores = clf.fit(X_n120).local_outlier_probabilities
+        scores = clf.fit(X_n120).local_outlier_probabilities_
 
     assert scores is not None
     assert len(scores) == len(X_n120)
@@ -865,14 +979,15 @@ def test_n_jobs_without_numba_warns(X_n120) -> None:
 
 def test_n_jobs_negative_two() -> None:
     """
-    Tests that n_jobs=-2 (invalid) produces a warning and defaults to 1.
+    Tests that n_jobs=-2 (invalid) produces a warning and uses an effective
+    value of 1 without mutating the public parameter (sklearn immutability).
     """
-    with pytest.warns(UserWarning) as record:
-        clf = loop.LocalOutlierProbability(n_neighbors=2, n_jobs=-2)
+    clf = loop.LocalOutlierProbability(n_neighbors=2, n_jobs=-2)
+    with pytest.warns(UserWarning, match="n_jobs must be -1 or a positive integer") as record:
+        clf.fit()
 
-    messages = [r.message.args[0] for r in record]
-    assert any("n_jobs must be -1 or a positive integer" in m for m in messages)
-    assert clf.n_jobs == 1
+    assert clf.n_jobs == -2
+    assert clf._fit_n_jobs == 1
 
 
 def test_vectorized_progress_bar_single_cluster(X_n120) -> None:
@@ -883,21 +998,21 @@ def test_vectorized_progress_bar_single_cluster(X_n120) -> None:
     clf = loop.LocalOutlierProbability(
         n_neighbors=10, n_jobs=1, progress_bar=True
     )
-    scores = clf.fit(X_n120).local_outlier_probabilities
+    scores = clf.fit(X_n120).local_outlier_probabilities_
     assert scores is not None
     assert len(scores) == len(X_n120)
 
 
 def test_n_jobs_invalid() -> None:
     """
-    Tests that invalid n_jobs values produce a warning and default to 1.
+    Tests that invalid n_jobs values produce a warning and use an effective
+    value of 1 without mutating the public parameter (sklearn immutability).
     """
-    with pytest.warns(UserWarning) as record:
-        clf = loop.LocalOutlierProbability(n_neighbors=2, n_jobs=0)
-
-    messages = [r.message.args[0] for r in record]
-    assert any("n_jobs must be -1 or a positive integer" in m for m in messages)
-    assert clf.n_jobs == 1
+    clf = loop.LocalOutlierProbability(n_neighbors=2, n_jobs=0)
+    with pytest.warns(UserWarning, match="n_jobs must be -1 or a positive integer") as record:
+        clf.fit()
+    assert clf.n_jobs == 0
+    assert clf._fit_n_jobs == 1
 
 
 # --- Numba-specific tests ---
@@ -918,10 +1033,10 @@ def test_numba_sequential_equivalence(X_n8) -> None:
     the default vectorized path.
     """
     clf_vec = loop.LocalOutlierProbability(n_neighbors=5, use_numba=False)
-    scores_vec = clf_vec.fit(X_n8).local_outlier_probabilities
+    scores_vec = clf_vec.fit(X_n8).local_outlier_probabilities_
 
     clf_numba = loop.LocalOutlierProbability(n_neighbors=5, use_numba=True)
-    scores_numba = clf_numba.fit(X_n8).local_outlier_probabilities
+    scores_numba = clf_numba.fit(X_n8).local_outlier_probabilities_
 
     assert_array_almost_equal(scores_vec, scores_numba, decimal=6)
 
@@ -941,14 +1056,14 @@ def test_numba_parallel_equivalence(X_n140_outliers) -> None:
     )
     scores_seq = clf_seq.fit(
         X_n140_outliers, cluster_labels=cluster_labels
-    ).local_outlier_probabilities
+    ).local_outlier_probabilities_
 
     clf_par = loop.LocalOutlierProbability(
         n_neighbors=10, use_numba=True, n_jobs=2
     )
     scores_par = clf_par.fit(
         X_n140_outliers, cluster_labels=cluster_labels
-    ).local_outlier_probabilities
+    ).local_outlier_probabilities_
 
     assert_array_almost_equal(scores_seq, scores_par, decimal=10)
 
@@ -961,7 +1076,7 @@ def test_numba_with_progress_bar(X_n120) -> None:
     clf = loop.LocalOutlierProbability(
         n_neighbors=10, use_numba=True, progress_bar=True
     )
-    scores = clf.fit(X_n120).local_outlier_probabilities
+    scores = clf.fit(X_n120).local_outlier_probabilities_
     assert scores is not None
     assert len(scores) == len(X_n120)
 
@@ -975,12 +1090,12 @@ def test_numba_prange_single_cluster(X_n120) -> None:
     clf_seq = loop.LocalOutlierProbability(
         n_neighbors=10, use_numba=True, n_jobs=1
     )
-    scores_seq = clf_seq.fit(X_n120).local_outlier_probabilities
+    scores_seq = clf_seq.fit(X_n120).local_outlier_probabilities_
 
     clf_par = loop.LocalOutlierProbability(
         n_neighbors=10, use_numba=True, n_jobs=-1
     )
-    scores_par = clf_par.fit(X_n120).local_outlier_probabilities
+    scores_par = clf_par.fit(X_n120).local_outlier_probabilities_
 
     assert_array_almost_equal(scores_seq, scores_par, decimal=6)
 
@@ -988,6 +1103,33 @@ def test_numba_prange_single_cluster(X_n120) -> None:
 # --- scikit-learn API convention tests (1.0.0) ---
 
 
+def test_exceptions_importable_from_loop() -> None:
+    """
+    Tests that the exception classes, which moved to ``PyNomaly.exceptions``
+    in 1.0.0, remain importable from ``PyNomaly.loop`` for backward
+    compatibility and are the same objects as the canonical ones.
+    """
+    from PyNomaly import exceptions
+    from PyNomaly.loop import (
+        ClusterSizeError as LoopClusterSizeError,
+        MissingValuesError as LoopMissingValuesError,
+        PyNomalyError as LoopPyNomalyError,
+        ValidationError as LoopValidationError,
+    )
+
+    assert LoopPyNomalyError is exceptions.PyNomalyError
+    assert LoopValidationError is exceptions.ValidationError
+    assert LoopClusterSizeError is exceptions.ClusterSizeError
+    assert LoopMissingValuesError is exceptions.MissingValuesError
+    for name in (
+        "PyNomalyError",
+        "ValidationError",
+        "ClusterSizeError",
+        "MissingValuesError",
+    ):
+        assert name in loop.__all__
+
+
 def test_numerical_equivalence_old_style(X_n8) -> None:
     """
     Tests that old-style (deprecated) and new-style API produce identical
@@ -996,10 +1138,10 @@ def test_numerical_equivalence_old_style(X_n8) -> None:
     with warnings.catch_warnings():
         warnings.simplefilter("ignore", FutureWarning)
         clf_old = loop.LocalOutlierProbability(data=X_n8, n_neighbors=5)
-        scores_old = clf_old.fit().local_outlier_probabilities
+        scores_old = clf_old.fit().local_outlier_probabilities_
 
     clf_new = loop.LocalOutlierProbability(n_neighbors=5)
-    scores_new = clf_new.fit(X_n8).local_outlier_probabilities
+    scores_new = clf_new.fit(X_n8).local_outlier_probabilities_
 
     assert_array_equal(scores_old, scores_new)
 
@@ -1022,10 +1164,51 @@ def test_fit_overrides_init(X_n8, X_n120) -> None:
         clf = loop.LocalOutlierProbability(data=X_n8, n_neighbors=5)
 
     clf.fit(X_n120)
-    scores = clf.local_outlier_probabilities
+    scores = clf.local_outlier_probabilities_
     assert len(scores) == 120
 
 
+def test_fit_with_2d_label_y_is_ignored(X_n120) -> None:
+    """
+    Tests that a (n, 1) label array passed as y (sklearn-style fit(X, y))
+    is ignored rather than being hijacked as a legacy positional
+    distance_matrix, which previously left the model unfitted.
+    """
+    y = np.ones((X_n120.shape[0], 1))
+    clf = loop.LocalOutlierProbability(n_neighbors=10).fit(X_n120, y)
+    assert clf.is_fit_ is True
+    assert clf.local_outlier_probabilities_ is not None
+    assert len(clf.local_outlier_probabilities_) == len(X_n120)
+
+
+@pytest.mark.parametrize(
+    "y_shape", [(120,), (120, 1), (120, 2)], ids=["1d", "n_by_1", "n_by_2"]
+)
+def test_fit_ignores_y_when_X_is_given(X_n120, y_shape) -> None:
+    """
+    Regression test for the fit(X, y) review finding on PR #92: when raw
+    data X is supplied, any label array y (1-d, a single column, or a
+    multi-column 2-d array) must be ignored without triggering the legacy
+    positional distance_matrix fallback, and the result must equal fit(X).
+    """
+    expected = (
+        loop.LocalOutlierProbability(n_neighbors=10)
+        .fit(X_n120)
+        .local_outlier_probabilities_
+    )
+
+    y = np.ones(y_shape)
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", FutureWarning)
+        clf = loop.LocalOutlierProbability(n_neighbors=10).fit(X_n120, y)
+
+    assert clf.is_fit_ is True
+    # raw-data path was used: y did not become the distance matrix
+    assert clf.data_ is X_n120
+    assert clf.distance_matrix is None
+    assert_array_equal(clf.local_outlier_probabilities_, expected)
+
+
 def test_refit(X_n8, X_n120) -> None:
     """
     Tests that calling fit() twice with different data produces correct
@@ -1034,19 +1217,144 @@ def test_refit(X_n8, X_n120) -> None:
     clf = loop.LocalOutlierProbability(n_neighbors=5)
 
     clf.fit(X_n8)
-    scores_first = clf.local_outlier_probabilities.copy()
+    scores_first = clf.local_outlier_probabilities_.copy()
     assert len(scores_first) == 8
 
     clf.fit(X_n120)
-    scores_second = clf.local_outlier_probabilities
+    scores_second = clf.local_outlier_probabilities_
     assert len(scores_second) == 120
-    assert clf.is_fit is True
+    assert clf.is_fit_ is True
 
 
 def test_loop_alias_import() -> None:
     """
     Tests that LoOP alias can be imported and is the same class.
     """
-    from PyNomaly import LoOP
-    from PyNomaly import LocalOutlierProbability
+    from PyNomaly import LocalOutlierProbability, LoOP
     assert LoOP is LocalOutlierProbability
+
+def test_legacy_positional_arguments(X_n120) -> None:
+    """
+    Tests that passing a distance matrix as the second positional argument (y)
+    is correctly routed to distance_matrix and raises a FutureWarning.
+    """
+    neigh = NearestNeighbors(metric="euclidean")
+    neigh.fit(X_n120)
+    d, idx = neigh.kneighbors(X_n120, n_neighbors=10, return_distance=True)
+
+    clf = loop.LocalOutlierProbability(n_neighbors=10)
+
+    match_str = re.escape("Passing 'distance_matrix' as the second positional argument is deprecated")
+    with pytest.warns(FutureWarning, match=match_str):
+        clf.fit(None, d, neighbor_matrix=idx)
+
+    assert clf.distance_matrix_ is not None
+    assert clf.local_outlier_probabilities_ is not None
+
+
+def test_standard_y_argument_ignored(X_n120) -> None:
+    """
+    Tests that passing standard 1D target labels to y does not trigger
+    the distance_matrix routing logic.
+    """
+    y_labels = np.zeros(120)
+
+    clf = loop.LocalOutlierProbability(n_neighbors=10)
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("error", FutureWarning)
+        clf.fit(X_n120, y_labels)
+
+    assert clf.local_outlier_probabilities_ is not None
+
+
+def test_decision_function(X_n20_scores) -> None:
+    """
+    Tests that decision_function follows the sklearn outlier convention:
+    ``score_samples(X) - offset_`` with ``offset_ = -0.5``, so inliers satisfy
+    ``decision_function(X) >= 0`` iff LoOP probability <= 0.5.
+    """
+    input_data, expected_scores = X_n20_scores
+    clf = loop.LocalOutlierProbability().fit(input_data)
+
+    expected_decision_scores = -1.0 * expected_scores - clf.offset_
+
+    decision_scores = clf.decision_function(input_data)
+
+    assert_array_almost_equal(decision_scores, expected_decision_scores, decimal=6)
+
+
+def test_predict(X_n20_scores) -> None:
+    """
+    Tests that predict returns exact 1 and -1 labels based on the 0.5 
+    probability threshold applied to the known historical fixture results.
+    """
+    input_data, expected_scores = X_n20_scores
+    clf = loop.LocalOutlierProbability().fit(input_data)
+    
+    # predict returns -1 for outliers (prob >= 0.5), 1 for inliers (prob < 0.5)
+    expected_predictions = np.where(expected_scores >= 0.5, -1, 1)
+    
+    predictions = clf.predict(input_data)
+    
+    assert_array_equal(predictions, expected_predictions)
+
+
+# --- Sparse matrix tests ---
+# Require scipy for sparse input. Skipped when scipy is not installed.
+
+_scipy = pytest.importorskip("scipy", reason="scipy required for sparse matrix tests")
+from scipy.sparse import csr_matrix, csc_matrix
+
+
+def test_sparse_fit_csr(X_n8) -> None:
+    """CSR sparse matrices can be passed to fit()."""
+    X_sparse = csr_matrix(X_n8, dtype=float)
+    clf = loop.LocalOutlierProbability(n_neighbors=5)
+    scores = clf.fit(X_sparse).local_outlier_probabilities_
+    assert scores is not None
+    assert len(scores) == X_n8.shape[0]
+
+
+def test_sparse_fit_csc_predict(X_n8) -> None:
+    """CSC sparse matrices work for fit() and predict()."""
+    X_sparse = csc_matrix(X_n8, dtype=float)
+    clf = loop.LocalOutlierProbability(n_neighbors=5).fit(X_sparse)
+    labels = clf.predict(X_sparse)
+    assert labels.shape == (X_n8.shape[0],)
+    assert set(labels).issubset({-1, 1})
+
+
+def test_sparse_equivalent_to_dense(X_n8) -> None:
+    """Sparse and dense inputs produce identical LoOP scores."""
+    clf_dense = loop.LocalOutlierProbability(n_neighbors=5).fit(X_n8)
+    clf_sparse = loop.LocalOutlierProbability(n_neighbors=5).fit(
+        csr_matrix(X_n8, dtype=float)
+    )
+    assert_array_almost_equal(
+        clf_dense.local_outlier_probabilities_,
+        clf_sparse.local_outlier_probabilities_,
+        decimal=6,
+    )
+
+
+def test_sparse_stream_row(X_n8) -> None:
+    """stream() accepts a single sparse observation row."""
+    clf = loop.LocalOutlierProbability(n_neighbors=5).fit(
+        csr_matrix(X_n8, dtype=float)
+    )
+    row_sparse = csr_matrix(X_n8[0:1], dtype=float)
+    score = clf.stream(row_sparse)
+    assert isinstance(score, (float, np.floating))
+    assert 0.0 <= score <= 1.0
+
+
+def test_sparse_requires_scipy(X_n8) -> None:
+    """Sparse input without scipy installed raises a clear ImportError."""
+    from unittest.mock import patch
+
+    X_sparse = csr_matrix(X_n8, dtype=float)
+    clf = loop.LocalOutlierProbability(n_neighbors=5)
+    with patch("PyNomaly._validation._SCIPY_AVAILABLE", False):
+        with pytest.raises(ImportError, match="Sparse matrix input requires scipy"):
+            clf.fit(X_sparse)