Repository: GAA-UAM/scikit-fda

Issue:
Unify scalar and functional regressors
I was mentioned in a meeting that it makes no sense to differentiate between scalar and functional regressors.  First, because there are not scalar, but multivariate, and second, because there is no confusion: the output of `predict` depends on the input to `fit`.

I forgot to mention this on the review of #112, but it should not be difficult to change. Can you do it @pablomm ?

Discussion:

Patch:
diff --git a/docs/modules/ml/classification.rst b/docs/modules/ml/classification.rst
index fb4e15f3..9524a4ae 100644
--- a/docs/modules/ml/classification.rst
+++ b/docs/modules/ml/classification.rst
@@ -11,10 +11,9 @@ Nearest Neighbors
 
 This module contains `nearest neighbors
 <https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm>`_ estimators to
-perform classification. In the examples `K-nearest neighbors classification
-<../../../auto_examples/plot_k_neighbors_classification.html>`_ and
-`Radius neighbors classification
-<../../../auto_examples/plot_radius_neighbors_classification.html>`_
+perform classification. In the examples
+:ref:`sphx_glr_auto_examples_plot_k_neighbors_classification.py`  and
+:ref:`sphx_glr_auto_examples_plot_radius_neighbors_classification.py`
 it is explained the basic usage of these estimators.
 
 .. autosummary::
diff --git a/docs/modules/ml/clustering.rst b/docs/modules/ml/clustering.rst
index 63979165..ce07f534 100644
--- a/docs/modules/ml/clustering.rst
+++ b/docs/modules/ml/clustering.rst
@@ -13,7 +13,7 @@ The following classes implement both, the K-Means and the Fuzzy K-Means
 algorithms respectively. In order to show the results in a visual way,
 the module :mod:`skfda.exploratory.visualization.clustering_plots
 <skfda.exploratory.visualization.clustering_plots>` can be used.
-See the `Clustering Example <../auto_examples/plot_clustering.html>`_ for a
+See the example :ref:`sphx_glr_auto_examples_plot_clustering.py` for a
 detailed explanation.
 
 .. autosummary::
diff --git a/docs/modules/ml/regression.rst b/docs/modules/ml/regression.rst
index 0191022b..72ba60f4 100644
--- a/docs/modules/ml/regression.rst
+++ b/docs/modules/ml/regression.rst
@@ -20,16 +20,13 @@ Nearest Neighbors
 
 This module contains `nearest neighbors
 <https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm>`_ estimators to
-perform regression. In the examples `Neighbors Scalar Regression
-<../../../auto_examples/plot_neighbors_scalar_regression.html>`_ and
-`Neighbors Functional Regression
-<../../../auto_examples/plot_neighbors_functional_regression.html>`_
+perform regression. In the examples
+:ref:`sphx_glr_auto_examples_plot_neighbors_scalar_regression.py` and
+:ref:`sphx_glr_auto_examples_plot_neighbors_functional_regression.py`
 it is explained the basic usage of these estimators.
 
 .. autosummary::
    :toctree: autosummary
 
-   skfda.ml.regression.KNeighborsScalarRegressor
-   skfda.ml.regression.RadiusNeighborsScalarRegressor
-   skfda.ml.regression.KNeighborsFunctionalRegressor
-   skfda.ml.regression.RadiusNeighborsFunctionalRegressor
+   skfda.ml.regression.KNeighborsRegressor
+   skfda.ml.regression.RadiusNeighborsRegressor
diff --git a/examples/plot_k_neighbors_classification.py b/examples/plot_k_neighbors_classification.py
index 1d9912de..1b4d2d41 100644
--- a/examples/plot_k_neighbors_classification.py
+++ b/examples/plot_k_neighbors_classification.py
@@ -185,7 +185,7 @@ ax.set_ylim((0.9, 1))
 # the integration, but the result should be similar.
 #
 
-knn = KNeighborsClassifier(metric='euclidean', sklearn_metric=True)
+knn = KNeighborsClassifier(metric='euclidean', multivariate_metric=True)
 gscv2 = GridSearchCV(knn, param_grid, cv=ss)
 gscv2.fit(X, y)
 
diff --git a/examples/plot_neighbors_functional_regression.py b/examples/plot_neighbors_functional_regression.py
index e5aac928..f3524c2c 100644
--- a/examples/plot_neighbors_functional_regression.py
+++ b/examples/plot_neighbors_functional_regression.py
@@ -11,10 +11,10 @@ Shows the usage of the nearest neighbors regressor with functional response.
 # sphinx_gallery_thumbnail_number = 4
 
 from sklearn.model_selection import train_test_split
-
 import matplotlib.pyplot as plt
+
 import skfda
-from skfda.ml.regression import KNeighborsFunctionalRegressor
+from skfda.ml.regression import KNeighborsRegressor
 from skfda.representation.basis import Fourier
 
 
@@ -22,9 +22,8 @@ from skfda.representation.basis import Fourier
 #
 # In this example we are going to show the usage of the nearest neighbors
 # regressors with functional response. There is available a K-nn version,
-# :class:`~skfda.ml.regression.KNeighborsFunctionalRegressor`, and other one
-# based in the radius,
-# :class:`~skfda.ml.regression.RadiusNeighborsFunctionalRegressor`.
+# :class:`~skfda.ml.regression.KNeighborsRegressor`, and other one based in
+# the radius, :class:`~skfda.ml.regression.RadiusNeighborsRegressor`.
 #
 #
 # As in the :ref:`scalar response example
@@ -81,7 +80,7 @@ X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.1,
 #
 
 
-knn = KNeighborsFunctionalRegressor(n_neighbors=5, weights='distance')
+knn = KNeighborsRegressor(n_neighbors=5, weights='distance')
 knn.fit(X_train, y_train)
 
 ##############################################################################
diff --git a/examples/plot_neighbors_scalar_regression.py b/examples/plot_neighbors_scalar_regression.py
index 7c232d6e..025121fb 100644
--- a/examples/plot_neighbors_scalar_regression.py
+++ b/examples/plot_neighbors_scalar_regression.py
@@ -14,17 +14,20 @@ from sklearn.model_selection import train_test_split, GridSearchCV, KFold
 
 import matplotlib.pyplot as plt
 import numpy as np
+
 import skfda
-from skfda.ml.regression import KNeighborsScalarRegressor
+from skfda.ml.regression import KNeighborsRegressor
+
 
 
 ##############################################################################
 #
 # In this example, we are going to show the usage of the nearest neighbors
 # regressors with scalar response. There is available a K-nn version,
-# :class:`~skfda.ml.regression.KNeighborsScalarRegressor`, and other one based
-# in the radius,
-# :class:`~skfda.ml.regression.RadiusNeighborsScalarRegressor`.
+# :class:`KNeighborsRegressor
+# <skfda.ml.regression.KNeighborsRegressor>`, and other one based in the
+# radius, :class:`RadiusNeighborsRegressor
+# <skfda.ml.regression.RadiusNeighborsRegressor>`.
 #
 # Firstly we will fetch a dataset to show the basic usage.
 #
@@ -80,15 +83,16 @@ X_train, X_test, y_train, y_test = train_test_split(X, log_prec,
 # Firstly we will try make a prediction with the default values of the
 # estimator, using 5 neighbors and the :math:`\mathbb{L}^2` distance.
 #
-# We can fit the :class:`~skfda.ml.regression.KNeighborsScalarRegressor` in
-# the same way than the scikit-learn estimators. This estimator is an
-# extension of the sklearn :class:`~sklearn.neighbors.KNeighborsRegressor`,
-# but accepting a :class:`~skfda.representation.grid.FDataGrid` as input
-# instead of an array with multivariate data.
+# We can fit the :class:`~skfda.ml.regression.KNeighborsRegressor` in the 
+# same way than the
+# sklearn estimators. This estimator is an extension of the sklearn
+# :class:`~sklearn.neighbors.KNeighborsRegressor`, but accepting a
+# :class:`~skfda.representation.grid.FDataGrid` as input instead of an array with
+# multivariate data.
 #
 
 
-knn = KNeighborsScalarRegressor(weights='distance')
+knn = KNeighborsRegressor(weights='distance')
 knn.fit(X_train, y_train)
 
 ##############################################################################
@@ -153,7 +157,7 @@ param_grid = {'n_neighbors': np.arange(1, 12, 2),
               'weights': ['uniform', 'distance']}
 
 
-knn = KNeighborsScalarRegressor(metric='euclidean', sklearn_metric=True)
+knn = KNeighborsRegressor(metric='euclidean', multivariate_metric=True)
 gscv = GridSearchCV(knn, param_grid, cv=KFold(n_splits=3,
                                               shuffle=True, random_state=0))
 gscv.fit(X, log_prec)
diff --git a/examples/plot_radius_neighbors_classification.py b/examples/plot_radius_neighbors_classification.py
index 91cdc20d..8591fe39 100644
--- a/examples/plot_radius_neighbors_classification.py
+++ b/examples/plot_radius_neighbors_classification.py
@@ -159,7 +159,7 @@ print(test_score)
 #
 
 radius_nn = RadiusNeighborsClassifier(radius=3, metric='euclidean',
-                                      weights='distance', sklearn_metric=True)
+                                      weights='distance', multivariate_metric=True)
 
 
 radius_nn.fit(X_train, y_train)
diff --git a/skfda/_neighbors/__init__.py b/skfda/_neighbors/__init__.py
index 9134bfb0..58316566 100644
--- a/skfda/_neighbors/__init__.py
+++ b/skfda/_neighbors/__init__.py
@@ -4,17 +4,11 @@ Includes the following classes:
  - KNeighborsClassifier
  - RadiusNeighborsClassifier
  - NearestCentroids
- - KNeighborsScalarRegressor
- - RadiusNeighborsScalarRegressor
- - KNeighborsFunctionalRegressor
- - RadiusNeighborsFunctionalRegressor
-"""
+ - KNeighborsRegressor
+ - RadiusNeighborsRegressor
 
+"""
 from .unsupervised import NearestNeighbors
-
+from .regression import KNeighborsRegressor, RadiusNeighborsRegressor
 from .classification import (KNeighborsClassifier, RadiusNeighborsClassifier,
                              NearestCentroids)
-from .regression import (KNeighborsFunctionalRegressor,
-                         KNeighborsScalarRegressor,
-                         RadiusNeighborsFunctionalRegressor,
-                         RadiusNeighborsScalarRegressor)
diff --git a/skfda/_neighbors/base.py b/skfda/_neighbors/base.py
index 2bf971d6..5e73364c 100644
--- a/skfda/_neighbors/base.py
+++ b/skfda/_neighbors/base.py
@@ -2,14 +2,13 @@
 
 from abc import ABC, abstractmethod
 
-import scipy.integrate
 from sklearn.base import BaseEstimator
-from sklearn.neighbors import NearestNeighbors as _NearestNeighbors
 from sklearn.utils.validation import check_is_fitted as sklearn_check_is_fitted
+from sklearn.base import RegressorMixin
 
 import numpy as np
 
-from .. import FDataGrid
+from .. import FDataGrid, FData
 from ..exploratory.stats import mean as l2_mean
 from ..misc.metrics import lp_distance
 
@@ -46,7 +45,7 @@ def _from_multivariate(data_matrix, sample_points, shape, **kwargs):
     return FDataGrid(data_matrix.reshape(shape), sample_points, **kwargs)
 
 
-def _to_sklearn_metric(metric, sample_points):
+def _to_multivariate_metric(metric, sample_points):
     r"""Transform a metric between FDatagrid in a sklearn compatible one.
 
     Given a metric between FDatagrids returns a compatible metric used to
@@ -66,7 +65,7 @@ def _to_sklearn_metric(metric, sample_points):
         >>> import numpy as np
         >>> from skfda import FDataGrid
         >>> from skfda.misc.metrics import lp_distance
-        >>> from skfda._neighbors.base import _to_sklearn_metric
+        >>> from skfda._neighbors.base import _to_multivariate_metric
 
         Calculate the Lp distance between fd and fd2.
 
@@ -78,7 +77,7 @@ def _to_sklearn_metric(metric, sample_points):
 
         Creation of the sklearn-style metric.
 
-        >>> sklearn_lp_distance = _to_sklearn_metric(lp_distance, [x])
+        >>> sklearn_lp_distance = _to_multivariate_metric(lp_distance, [x])
         >>> sklearn_lp_distance(np.ones(len(x)), np.zeros(len(x))).round(2)
         1.0
 
@@ -86,13 +85,13 @@ def _to_sklearn_metric(metric, sample_points):
     # Shape -> (n_samples = 1, domain_dims...., image_dimension (-1))
     shape = [1] + [len(axis) for axis in sample_points] + [-1]
 
-    def sklearn_metric(x, y, _check=False, **kwargs):
+    def multivariate_metric(x, y, _check=False, **kwargs):
 
         return metric(_from_multivariate(x, sample_points, shape),
                       _from_multivariate(y, sample_points, shape),
                       _check=_check, **kwargs)
 
-    return sklearn_metric
+    return multivariate_metric
 
 
 class NeighborsBase(ABC, BaseEstimator):
@@ -102,7 +101,7 @@ class NeighborsBase(ABC, BaseEstimator):
     def __init__(self, n_neighbors=None, radius=None,
                  weights='uniform', algorithm='auto',
                  leaf_size=30, metric='l2', metric_params=None,
-                 n_jobs=None, sklearn_metric=False):
+                 n_jobs=None, multivariate_metric=False):
 
         self.n_neighbors = n_neighbors
         self.radius = radius
@@ -112,7 +111,7 @@ class NeighborsBase(ABC, BaseEstimator):
         self.metric = metric
         self.metric_params = metric_params
         self.n_jobs = n_jobs
-        self.sklearn_metric = sklearn_metric
+        self.multivariate_metric = multivariate_metric
 
     def _check_is_fitted(self):
         """Check if the estimator is fitted.
@@ -161,17 +160,18 @@ class NeighborsMixin:
             self._sample_points = X.sample_points
             self._shape = X.data_matrix.shape[1:]
 
-            if not self.sklearn_metric:
+            if not self.multivariate_metric:
                 # Constructs sklearn metric to manage vector
                 if self.metric == 'l2':
                     metric = lp_distance
                 else:
                     metric = self.metric
-                sk_metric = _to_sklearn_metric(metric, self._sample_points)
+                sklearn_metric = _to_multivariate_metric(metric,
+                                                         self._sample_points)
             else:
-                sk_metric = self.metric
+                sklearn_metric = self.metric
 
-            self.estimator_ = self._init_estimator(sk_metric)
+            self.estimator_ = self._init_estimator(sklearn_metric)
             self.estimator_.fit(self._transform_to_multivariate(X), y)
 
         return self
@@ -436,57 +436,32 @@ class NeighborsClassifierMixin:
         return self.estimator_.predict(X)
 
 
-class NeighborsScalarRegresorMixin:
-    """Mixin class for scalar regressor based in nearest neighbors"""
+class NeighborsRegressorMixin(NeighborsMixin, RegressorMixin):
+    """Mixin class for the regressors based on neighbors"""
 
-    def predict(self, X):
-        """Predict the target for the provided data
-        Parameters
-        ----------
-        X (:class:`FDataGrid` or array-like): FDataGrid with the test
-            samples or array (n_query, n_indexed) if metric ==
-            'precomputed'.
-        Returns
-        -------
-        y : array of int, shape = [n_samples] or [n_samples, n_outputs]
-            Target values
-        Notes
-        -----
-        This method wraps the corresponding sklearn routine in the module
-        ``sklearn.neighbors``.
-
-        """
-        self._check_is_fitted()
-
-        X = self._transform_to_multivariate(X)
-
-        return self.estimator_.predict(X)
-
-
-class NearestNeighborsMixinInit:
-    def _init_estimator(self, sk_metric):
-        """Initialize the sklearn nearest neighbors estimator.
+    def fit(self, X, y):
+        """Fit the model using X as training data and y as responses.
 
         Args:
-            sk_metric: (pyfunc or 'precomputed'): Metric compatible with
-                sklearn API or matrix (n_samples, n_samples) with precomputed
-                distances.
-
+            X (:class:`FDataGrid`, array_matrix): Training data. FDataGrid
+                with the training data or array matrix with shape
+                [n_samples, n_samples] if metric='precomputed'.
+            Y (:class:`FData` or array_like): Training data. FData
+                with the training respones (functional response case)
+                or array matrix with length `n_samples` in the multivariate
+                response case.
         Returns:
-            Sklearn K Neighbors estimator initialized.
+            Estimator: self.
 
         """
-        return _NearestNeighbors(
-            n_neighbors=self.n_neighbors, radius=self.radius,
-            algorithm=self.algorithm, leaf_size=self.leaf_size,
-            metric=sk_metric, metric_params=self.metric_params,
-            n_jobs=self.n_jobs)
+        self._functional = isinstance(y, FData)
 
+        if self._functional:
+            return self._functional_fit(X, y)
+        else:
+            return super().fit(X, y)
 
-class NeighborsFunctionalRegressorMixin:
-    """Mixin class for the functional regressors based in neighbors"""
-
-    def fit(self, X, y):
+    def _functional_fit(self, X, y):
         """Fit the model using X as training data.
 
         Args:
@@ -508,18 +483,19 @@ class NeighborsFunctionalRegressorMixin:
             self._sample_points = X.sample_points
             self._shape = X.data_matrix.shape[1:]
 
-            if not self.sklearn_metric:
+            if not self.multivariate_metric:
 
                 if self.metric == 'l2':
                     metric = lp_distance
                 else:
                     metric = self.metric
 
-                sk_metric = _to_sklearn_metric(metric, self._sample_points)
+                sklearn_metric = _to_multivariate_metric(metric,
+                                                         self._sample_points)
             else:
-                sk_metric = self.metric
+                sklearn_metric = self.metric
 
-            self.estimator_ = self._init_estimator(sk_metric)
+            self.estimator_ = self._init_estimator(sklearn_metric)
             self.estimator_.fit(self._transform_to_multivariate(X))
 
         if self.regressor == 'mean':
@@ -529,11 +505,11 @@ class NeighborsFunctionalRegressorMixin:
 
         # Choose proper local regressor
         if self.weights == 'uniform':
-            self.local_regressor = self._uniform_local_regression
+            self._local_regressor = self._uniform_local_regression
         elif self.weights == 'distance':
-            self.local_regressor = self._distance_local_regression
+            self._local_regressor = self._distance_local_regression
         else:
-            self.local_regressor = self._weighted_local_regression
+            self._local_regressor = self._weighted_local_regression
 
         # Store the responses
         self._y = y
@@ -565,6 +541,71 @@ class NeighborsFunctionalRegressorMixin:
         return self._regressor(neighbors, weights)
 
     def predict(self, X):
+        """Predict the target for the provided data
+        Parameters
+        ----------
+        X (:class:`FDataGrid` or array-like): FDataGrid with the test
+            samples or array (n_query, n_indexed) if metric ==
+            'precomputed'.
+        Returns
+        -------
+        y : array of shape = [n_samples] or [n_samples, n_outputs]
+            or :class:`FData` containing as many samples as X.
+
+        """
+        self._check_is_fitted()
+
+        # Choose type of prediction
+        if self._functional:
+            return self._functional_predict(X)
+        else:
+            return self._multivariate_predict(X)
+
+    def _multivariate_predict(self, X):
+        """Predict the target for the provided data
+        Parameters
+        ----------
+        X (:class:`FDataGrid` or array-like): FDataGrid with the test
+            samples or array (n_query, n_indexed) if metric ==
+            'precomputed'.
+        Returns
+        -------
+        y : array of int, shape = [n_samples] or [n_samples, n_outputs]
+            Target values
+        Notes
+        -----
+        This method wraps the corresponding sklearn routine in the module
+        ``sklearn.neighbors``.
+
+        """
+        X = self._transform_to_multivariate(X)
+
+        return self.estimator_.predict(X)
+
+    def _init_estimator(self, sklearn_metric):
+        """Initialize the sklearn nearest neighbors estimator.
+
+        Args:
+            sklearn_metric: (pyfunc or 'precomputed'): Metric compatible with
+                sklearn API or matrix (n_samples, n_samples) with precomputed
+                distances.
+
+        Returns:
+            Sklearn K Neighbors estimator initialized.
+
+        """
+        if self._functional:
+            from sklearn.neighbors import NearestNeighbors as _NearestNeighbors
+
+            return _NearestNeighbors(
+                n_neighbors=self.n_neighbors, radius=self.radius,
+                algorithm=self.algorithm, leaf_size=self.leaf_size,
+                metric=sklearn_metric, metric_params=self.metric_params,
+                n_jobs=self.n_jobs)
+        else:
+            return self._init_multivariate_estimator(sklearn_metric)
+
+    def _functional_predict(self, X):
         """Predict functional responses.
 
         Args:
@@ -577,7 +618,6 @@ class NeighborsFunctionalRegressorMixin:
             y : :class:`FDataGrid` containing as many samples as X.
 
         """
-        self._check_is_fitted()
 
         X = self._transform_to_multivariate(X)
 
@@ -586,13 +626,14 @@ class NeighborsFunctionalRegressorMixin:
         if len(neighbors[0]) == 0:
             pred = self._outlier_response(neighbors)
         else:
-            pred = self.local_regressor(self._y[neighbors[0]], distances[0])
+            pred = self._local_regressor(self._y[neighbors[0]], distances[0])
 
         for i, idx in enumerate(neighbors[1:]):
             if len(idx) == 0:
                 new_pred = self._outlier_response(neighbors)
             else:
-                new_pred = self.local_regressor(self._y[idx], distances[i + 1])
+                new_pred = self._local_regressor(self._y[idx],
+                                                 distances[i + 1])
 
             pred = pred.concatenate(new_pred)
 
@@ -613,6 +654,51 @@ class NeighborsFunctionalRegressorMixin:
             return self.outlier_response
 
     def score(self, X, y, sample_weight=None):
+        r"""Return the coefficient of determination R^2 of the prediction.
+
+        In the multivariate response case, the coefficient :math:`R^2` is
+        defined as
+
+        .. math::
+            1 - \frac{\sum_{i=1}^{n} (y_i - \hat y_i)^2}
+            {\sum_{i=1}^{n} (y_i - \frac{1}{n}\sum_{i=1}^{n}y_i)^2}
+
+        where :math:`\hat{y}_i` is the prediction associated to the test sample
+        :math:`X_i`, and :math:`{y}_i` is the true response. See
+        :func:`sklearn.metrics.r2_score <sklearn.metrics.r2_score>` for more
+        information.
+
+
+        In the functional case it is returned an extension of the coefficient
+        of determination :math:`R^2`, defined as
+
+        .. math::
+            1 - \frac{\sum_{i=1}^{n}\int (y_i(t) - \hat{y}_i(t))^2dt}
+            {\sum_{i=1}^{n} \int (y_i(t)- \frac{1}{n}\sum_{i=1}^{n}y_i(t))^2dt}
+
+
+        The best possible score is 1.0 and it can be negative
+        (because the model can be arbitrarily worse). A constant model that
+        always predicts the expected value of y, disregarding the input
+        features, would get a R^2 score of 0.0.
+
+        Args:
+            X (FDataGrid): Test samples to be predicted.
+            y (FData or array-like): True responses of the test samples.
+            sample_weight (array_like, shape = [n_samples], optional): Sample
+                weights.
+
+        Returns:
+            (float): Coefficient of determination.
+
+        """
+        if self._functional:
+            return self._functional_score(X, y, sample_weight=sample_weight)
+        else:
+            # Default sklearn multivariate score
+            return super().score(X, y, sample_weight=sample_weight)
+
+    def _functional_score(self, X, y, sample_weight=None):
         r"""Return an extension of the coefficient of determination R^2.
 
         The coefficient is defined as
@@ -639,6 +725,11 @@ class NeighborsFunctionalRegressorMixin:
             (float): Coefficient of determination.
 
         """
+
+        # TODO: If it is created a module in ml.regression with other
+        # score metrics, move it.
+        from scipy.integrate import simps
+
         if y.dim_codomain != 1 or y.dim_domain != 1:
             raise ValueError("Score not implemented for multivariate "
                              "functional data.")
@@ -676,7 +767,7 @@ class NeighborsFunctionalRegressorMixin:
         sum_u = np.sum(data_u, axis=0)
         sum_v = np.sum(data_v, axis=0)
 
-        int_u = scipy.integrate.simps(sum_u, x=u.sample_points[0])
-        int_v = scipy.integrate.simps(sum_v, x=v.sample_points[0])
+        int_u = simps(sum_u, x=u.sample_points[0])
+        int_v = simps(sum_v, x=v.sample_points[0])
 
         return 1 - int_u / int_v
diff --git a/skfda/_neighbors/classification.py b/skfda/_neighbors/classification.py
index 1034725c..c8f63482 100644
--- a/skfda/_neighbors/classification.py
+++ b/skfda/_neighbors/classification.py
@@ -1,19 +1,15 @@
 """Neighbor models for supervised classification."""
 
-from .base import (NeighborsBase, NeighborsMixin, KNeighborsMixin,
-                   NeighborsClassifierMixin, RadiusNeighborsMixin)
 
 from sklearn.utils.multiclass import check_classification_targets
 from sklearn.preprocessing import LabelEncoder
 from sklearn.base import ClassifierMixin, BaseEstimator
 from sklearn.utils.validation import check_is_fitted as sklearn_check_is_fitted
 
-from sklearn.neighbors import KNeighborsClassifier as _KNeighborsClassifier
-from sklearn.neighbors import (RadiusNeighborsClassifier as
-                               _RadiusNeighborsClassifier)
-
 from ..misc.metrics import lp_distance, pairwise_distance
 from ..exploratory.stats import mean as l2_mean
+from .base import (NeighborsBase, NeighborsMixin, KNeighborsMixin,
+                   NeighborsClassifierMixin, RadiusNeighborsMixin)
 
 
 class KNeighborsClassifier(NeighborsBase, NeighborsMixin, KNeighborsMixin,
@@ -61,7 +57,7 @@ class KNeighborsClassifier(NeighborsBase, NeighborsMixin, KNeighborsMixin,
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors.
         Doesn't affect :meth:`fit` method.
-    sklearn_metric : boolean, optional (default = False)
+    multivariate_metric : boolean, optional (default = False)
         Indicates if the metric used is a sklearn distance between vectors (see
         :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
         the module :mod:`skfda.misc.metrics`.
@@ -95,11 +91,12 @@ class KNeighborsClassifier(NeighborsBase, NeighborsMixin, KNeighborsMixin,
 
     See also
     --------
-    RadiusNeighborsClassifier
-    KNeighborsScalarRegressor
-    RadiusNeighborsScalarRegressor
-    NearestNeighbors
-    NearestCentroids
+    :class:`~skfda.ml.classification.RadiusNeighborsClassifier`
+    :class:`~skfda.ml.classification.NearestCentroids`
+    :class:`~skfda.ml.regression.KNeighborsRegressor`
+    :class:`~skfda.ml.regression.RadiusNeighborsRegressor`
+    :class:`~skfda.ml.clustering.NearestNeighbors`
+
     Notes
     -----
     See Nearest Neighbors in the sklearn online documentation for a discussion
@@ -120,20 +117,20 @@ class KNeighborsClassifier(NeighborsBase, NeighborsMixin, KNeighborsMixin,
 
     def __init__(self, n_neighbors=5, weights='uniform', algorithm='auto',
                  leaf_size=30, metric='l2', metric_params=None,
-                 n_jobs=1, sklearn_metric=False):
+                 n_jobs=1, multivariate_metric=False):
         """Initialize the classifier."""
 
         super().__init__(n_neighbors=n_neighbors,
                          weights=weights, algorithm=algorithm,
                          leaf_size=leaf_size, metric=metric,
                          metric_params=metric_params, n_jobs=n_jobs,
-                         sklearn_metric=sklearn_metric)
+                         multivariate_metric=multivariate_metric)
 
-    def _init_estimator(self, sk_metric):
+    def _init_estimator(self, sklearn_metric):
         """Initialize the sklearn K neighbors estimator.
 
         Args:
-            sk_metric: (pyfunc or 'precomputed'): Metric compatible with
+            sklearn_metric: (pyfunc or 'precomputed'): Metric compatible with
                 sklearn API or matrix (n_samples, n_samples) with precomputed
                 distances.
 
@@ -141,10 +138,13 @@ class KNeighborsClassifier(NeighborsBase, NeighborsMixin, KNeighborsMixin,
             Sklearn K Neighbors estimator initialized.
 
         """
+        from sklearn.neighbors import (KNeighborsClassifier as
+                                       _KNeighborsClassifier)
+
         return _KNeighborsClassifier(
             n_neighbors=self.n_neighbors, weights=self.weights,
             algorithm=self.algorithm, leaf_size=self.leaf_size,
-            metric=sk_metric, metric_params=self.metric_params,
+            metric=sklearn_metric, metric_params=self.metric_params,
             n_jobs=self.n_jobs)
 
     def predict_proba(self, X):
@@ -219,7 +219,7 @@ class RadiusNeighborsClassifier(NeighborsBase, NeighborsMixin,
         The number of parallel jobs to run for neighbors search.
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors.
-    sklearn_metric : boolean, optional (default = False)
+    multivariate_metric : boolean, optional (default = False)
         Indicates if the metric used is a sklearn distance between vectors (see
         :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
         the module :mod:`skfda.misc.metrics`.
@@ -248,11 +248,11 @@ class RadiusNeighborsClassifier(NeighborsBase, NeighborsMixin,
 
     See also
     --------
-    KNeighborsClassifier
-    KNeighborsScalarRegressor
-    RadiusNeighborsScalarRegressor
-    NearestNeighbors
-    NearestCentroids
+    :class:`~skfda.ml.classification.KNeighborsClassifier`
+    :class:`~skfda.ml.classification.NearestCentroids`
+    :class:`~skfda.ml.regression.KNeighborsRegressor`
+    :class:`~skfda.ml.regression.RadiusNeighborsRegressor`
+    :class:`~skfda.ml.clustering.NearestNeighbors`
 
     Notes
     -----
@@ -268,21 +268,21 @@ class RadiusNeighborsClassifier(NeighborsBase, NeighborsMixin,
 
     def __init__(self, radius=1.0, weights='uniform', algorithm='auto',
                  leaf_size=30, metric='l2', metric_params=None,
-                 outlier_label=None, n_jobs=1, sklearn_metric=False):
+                 outlier_label=None, n_jobs=1, multivariate_metric=False):
         """Initialize the classifier."""
 
         super().__init__(radius=radius, weights=weights, algorithm=algorithm,
                          leaf_size=leaf_size, metric=metric,
                          metric_params=metric_params, n_jobs=n_jobs,
-                         sklearn_metric=sklearn_metric)
+                         multivariate_metric=multivariate_metric)
 
         self.outlier_label = outlier_label
 
-    def _init_estimator(self, sk_metric):
+    def _init_estimator(self, sklearn_metric):
         """Initialize the sklearn radius neighbors estimator.
 
         Args:
-            sk_metric: (pyfunc or 'precomputed'): Metric compatible with
+            sklearn_metric: (pyfunc or 'precomputed'): Metric compatible with
                 sklearn API or matrix (n_samples, n_samples) with precomputed
                 distances.
 
@@ -290,10 +290,13 @@ class RadiusNeighborsClassifier(NeighborsBase, NeighborsMixin,
             Sklearn Radius Neighbors estimator initialized.
 
         """
+        from sklearn.neighbors import (RadiusNeighborsClassifier as
+                                       _RadiusNeighborsClassifier)
+
         return _RadiusNeighborsClassifier(
             radius=self.radius, weights=self.weights,
             algorithm=self.algorithm, leaf_size=self.leaf_size,
-            metric=sk_metric, metric_params=self.metric_params,
+            metric=sklearn_metric, metric_params=self.metric_params,
             outlier_label=self.outlier_label, n_jobs=self.n_jobs)
 
 
@@ -349,11 +352,11 @@ class NearestCentroids(BaseEstimator, ClassifierMixin):
 
     See also
     --------
-    KNeighborsClassifier
-    RadiusNeighborsClassifier
-    KNeighborsScalarRegressor
-    RadiusNeighborsScalarRegressor
-    NearestNeighbors
+    :class:`~skfda.ml.classification.KNeighborsClassifier`
+    :class:`~skfda.ml.classification.RadiusNeighborsClassifier`
+    :class:`~skfda.ml.regression.KNeighborsRegressor`
+    :class:`~skfda.ml.regression.RadiusNeighborsRegressor`
+    :class:`~skfda.ml.clustering.NearestNeighbors`
 
     """
 
diff --git a/skfda/_neighbors/regression.py b/skfda/_neighbors/regression.py
index af16b44d..8300215e 100644
--- a/skfda/_neighbors/regression.py
+++ b/skfda/_neighbors/regression.py
@@ -1,21 +1,14 @@
 """Neighbor models for regression."""
 
-from sklearn.neighbors import KNeighborsRegressor as _KNeighborsRegressor
-from sklearn.neighbors import (RadiusNeighborsRegressor as
-                               _RadiusNeighborsRegressor)
-from sklearn.base import RegressorMixin
+from .base import (NeighborsBase, KNeighborsMixin, RadiusNeighborsMixin,
+                   NeighborsRegressorMixin)
 
-from .base import (NeighborsBase, NeighborsMixin,
-                   KNeighborsMixin, RadiusNeighborsMixin,
-                   NeighborsScalarRegresorMixin,
-                   NeighborsFunctionalRegressorMixin,
-                   NearestNeighborsMixinInit)
 
+class KNeighborsRegressor(NeighborsBase, NeighborsRegressorMixin,
+                          KNeighborsMixin):
+    """Regression based on k-nearest neighbors.
 
-class KNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
-                                KNeighborsMixin, RegressorMixin,
-                                NeighborsScalarRegresorMixin):
-    """Regression based on k-nearest neighbors with scalar response.
+    Regression with scalar, multivariate or functional response.
 
     The target is predicted by local interpolation of the targets associated of
     the nearest neighbors in the training set.
@@ -36,6 +29,11 @@ class KNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
           array of distances, and returns an array of the same shape
           containing the weights.
 
+    regressor : callable, optional ((default =
+        :func:`mean <skfda.exploratory.stats.mean>`))
+        Function to perform the local regression in the functional response
+        case. By default used the mean. Can the neighbors of a test sample,
+        and if weights != 'uniform' an array of weights as second parameter.
     algorithm : {'auto', 'ball_tree', 'brute'}, optional
         Algorithm used to compute the nearest neighbors:
 
@@ -61,7 +59,7 @@ class KNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors.
         Doesn't affect :meth:`fit` method.
-    sklearn_metric : boolean, optional (default = False)
+    multivariate_metric : boolean, optional (default = False)
         Indicates if the metric used is a sklearn distance between vectors (see
         :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
         the module :mod:`skfda.misc.metrics`.
@@ -69,31 +67,50 @@ class KNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
     --------
     Firstly, we will create a toy dataset with gaussian-like samples shifted.
 
+    >>> from skfda.ml.regression import KNeighborsRegressor
     >>> from skfda.datasets import make_multimodal_samples
     >>> from skfda.datasets import make_multimodal_landmarks
     >>> y = make_multimodal_landmarks(n_samples=30, std=.5, random_state=0)
-    >>> y = y.flatten()
-    >>> fd = make_multimodal_samples(n_samples=30, std=.5, random_state=0)
+    >>> y_train = y.flatten()
+    >>> X_train = make_multimodal_samples(n_samples=30, std=.5, random_state=0)
+    >>> X_test = make_multimodal_samples(n_samples=5, std=.05, random_state=0)
 
     We will fit a K-Nearest Neighbors regressor to regress a scalar response.
 
-    >>> from skfda.ml.regression import KNeighborsScalarRegressor
-    >>> neigh = KNeighborsScalarRegressor()
-    >>> neigh.fit(fd, y)
-    KNeighborsScalarRegressor(algorithm='auto', leaf_size=30,...)
+    >>> neigh = KNeighborsRegressor()
+    >>> neigh.fit(X_train, y_train)
+    KNeighborsRegressor(algorithm='auto', leaf_size=30,...)
 
     We can predict the modes of new samples
 
-    >>> neigh.predict(fd[:4]).round(2) # Predict first 4 locations
-    array([ 0.79,  0.27,  0.71,  0.79])
+    >>> neigh.predict(X_test).round(2) # Predict test data
+    array([ 0.38, 0.14, 0.27, 0.52, 0.38])
+
+
+    Now we will create a functional response to train the model
+
+    >>> y_train = 5 * X_train + 1
+    >>> y_train
+    FDataGrid(...)
+
+    We train the estimator with the functional response
+
+    >>> neigh.fit(X_train, y_train)
+    KNeighborsRegressor(algorithm='auto', leaf_size=30,...)
+
+    And predict the responses as in the first case.
+
+    >>> neigh.predict(X_test)
+    FDataGrid(...)
 
     See also
     --------
-    KNeighborsClassifier
-    RadiusNeighborsClassifier
-    RadiusNeighborsScalarRegressor
-    NearestNeighbors
-    NearestCentroids
+    :class:`~skfda.ml.classification.KNeighborsClassifier`
+    :class:`~skfda.ml.classification.RadiusNeighborsClassifier`
+    :class:`~skfda.ml.classification.NearestCentroids`
+    :class:`~skfda.ml.regression.RadiusNeighborsRegressor`
+    :class:`~skfda.ml.clustering.NearestNeighbors`
+
     Notes
     -----
     See Nearest Neighbors in the sklearn online documentation for a discussion
@@ -112,22 +129,23 @@ class KNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
 
     """
 
-    def __init__(self, n_neighbors=5, weights='uniform', algorithm='auto',
-                 leaf_size=30, metric='l2', metric_params=None,
-                 n_jobs=1, sklearn_metric=False):
-        """Initialize the classifier."""
+    def __init__(self, n_neighbors=5, weights='uniform', regressor='mean',
+                 algorithm='auto', leaf_size=30, metric='l2',
+                 metric_params=None, n_jobs=1, multivariate_metric=False):
+        """Initialize the regressor."""
 
         super().__init__(n_neighbors=n_neighbors,
                          weights=weights, algorithm=algorithm,
                          leaf_size=leaf_size, metric=metric,
                          metric_params=metric_params, n_jobs=n_jobs,
-                         sklearn_metric=sklearn_metric)
+                         multivariate_metric=multivariate_metric)
+        self.regressor = regressor
 
-    def _init_estimator(self, sk_metric):
+    def _init_multivariate_estimator(self, sklearn_metric):
         """Initialize the sklearn K neighbors estimator.
 
         Args:
-            sk_metric: (pyfunc or 'precomputed'): Metric compatible with
+            sklearn_metric: (pyfunc or 'precomputed'): Metric compatible with
                 sklearn API or matrix (n_samples, n_samples) with precomputed
                 distances.
 
@@ -135,17 +153,25 @@ class KNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
             Sklearn K Neighbors estimator initialized.
 
         """
+        from sklearn.neighbors import (KNeighborsRegressor as
+                                       _KNeighborsRegressor)
+
         return _KNeighborsRegressor(
             n_neighbors=self.n_neighbors, weights=self.weights,
             algorithm=self.algorithm, leaf_size=self.leaf_size,
-            metric=sk_metric, metric_params=self.metric_params,
+            metric=sklearn_metric, metric_params=self.metric_params,
             n_jobs=self.n_jobs)
 
+    def _query(self, X):
+        """Return distances and neighbors of given sample."""
+        return self.estimator_.kneighbors(X)
+
+
+class RadiusNeighborsRegressor(NeighborsBase, NeighborsRegressorMixin,
+                               RadiusNeighborsMixin):
+    """Regression based on neighbors within a fixed radius.
 
-class RadiusNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
-                                     RadiusNeighborsMixin, RegressorMixin,
-                                     NeighborsScalarRegresorMixin):
-    """Scalar regression based on neighbors within a fixed radius.
+    Regression with scalar, multivariate or functional response.
 
     The target is predicted by local interpolation of the targets associated of
     the nearest neighbors in the training set.
@@ -168,6 +194,11 @@ class RadiusNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
           containing the weights.
 
         Uniform weights are used by default.
+    regressor : callable, optional ((default =
+        :func:`mean <skfda.exploratory.stats.mean>`))
+        Function to perform the local regression in the functional response
+        case. By default used the mean. Can the neighbors of a test sample,
+        and if weights != 'uniform' an array of weights as second parameter.
     algorithm : {'auto', 'ball_tree', 'brute'}, optional
         Algorithm used to compute the nearest neighbors:
 
@@ -188,174 +219,67 @@ class RadiusNeighborsScalarRegressor(NeighborsBase, NeighborsMixin,
         for a list of available metrics.
     metric_params : dict, optional (default = None)
         Additional keyword arguments for the metric function.
+    outlier_response : :class:`FData`, optional (default = None)
+        Default response in the functional response case for test samples
+        without neighbors.
     n_jobs : int or None, optional (default=None)
         The number of parallel jobs to run for neighbors search.
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors.
-    sklearn_metric : boolean, optional (default = False)
+    multivariate_metric : boolean, optional (default = False)
         Indicates if the metric used is a sklearn distance between vectors (see
         :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
         the module :mod:`skfda.misc.metrics`.
+
     Examples
     --------
     Firstly, we will create a toy dataset with gaussian-like samples shifted.
 
+    >>> from skfda.ml.regression import RadiusNeighborsRegressor
     >>> from skfda.datasets import make_multimodal_samples
     >>> from skfda.datasets import make_multimodal_landmarks
     >>> y = make_multimodal_landmarks(n_samples=30, std=.5, random_state=0)
-    >>> y = y.flatten()
-    >>> fd = make_multimodal_samples(n_samples=30, std=.5, random_state=0)
-
-
-    We will fit a K-Nearest Neighbors regressor to regress a scalar response.
-
-    >>> from skfda.ml.regression import RadiusNeighborsScalarRegressor
-    >>> neigh = RadiusNeighborsScalarRegressor(radius=.2)
-    >>> neigh.fit(fd, y)
-    RadiusNeighborsScalarRegressor(algorithm='auto', leaf_size=30,...)
-
-    We can predict the modes of new samples.
-
-    >>> neigh.predict(fd[:4]).round(2) # Predict first 4 locations
-    array([ 0.84,  0.27,  0.66,  0.79])
-
-    See also
-    --------
-    KNeighborsClassifier
-    RadiusNeighborsClassifier
-    KNeighborsScalarRegressor
-    NearestNeighbors
-    NearestCentroids
-    Notes
-    -----
-    See Nearest Neighbors in the sklearn online documentation for a discussion
-    of the choice of ``algorithm`` and ``leaf_size``.
-
-    This class wraps the sklearn classifier
-    `sklearn.neighbors.RadiusNeighborsClassifier`.
-
-    https://en.wikipedia.org/wiki/K-nearest_neighbor_algorithm
-
-    """
-
-    def __init__(self, radius=1.0, weights='uniform', algorithm='auto',
-                 leaf_size=30, metric='l2', metric_params=None,
-                 n_jobs=1, sklearn_metric=False):
-        """Initialize the classifier."""
-
-        super().__init__(radius=radius, weights=weights, algorithm=algorithm,
-                         leaf_size=leaf_size, metric=metric,
-                         metric_params=metric_params, n_jobs=n_jobs,
-                         sklearn_metric=sklearn_metric)
-
-    def _init_estimator(self, sk_metric):
-        """Initialize the sklearn radius neighbors estimator.
-
-        Args:
-            sk_metric: (pyfunc or 'precomputed'): Metric compatible with
-                sklearn API or matrix (n_samples, n_samples) with precomputed
-                distances.
-
-        Returns:
-            Sklearn Radius Neighbors estimator initialized.
-
-        """
-        return _RadiusNeighborsRegressor(
-            radius=self.radius, weights=self.weights,
-            algorithm=self.algorithm, leaf_size=self.leaf_size,
-            metric=sk_metric, metric_params=self.metric_params,
-            n_jobs=self.n_jobs)
-
-
-class KNeighborsFunctionalRegressor(NearestNeighborsMixinInit,
-                                    NeighborsBase, KNeighborsMixin,
-                                    NeighborsFunctionalRegressorMixin):
-    """Functional regression based on neighbors within a fixed radius.
+    >>> y_train = y.flatten()
+    >>> X_train = make_multimodal_samples(n_samples=30, std=.5, random_state=0)
+    >>> X_test = make_multimodal_samples(n_samples=5, std=.05, random_state=0)
 
-    The target is predicted by local interpolation of the targets
-    associated of the nearest neighbors in the training set.
+    We will fit a Radius-Nearest Neighbors regressor to regress a scalar
+    response.
 
-    Parameters
-    ----------
-    n_neighbors : int, optional (default = 5)
-        Number of neighbors to use by default for :meth:`kneighbors` queries.
-    weights : str or callable
-        weight function used in prediction.  Possible values:
+    >>> neigh = RadiusNeighborsRegressor(radius=0.2)
+    >>> neigh.fit(X_train, y_train)
+    RadiusNeighborsRegressor(algorithm='auto', leaf_size=30,...)
 
-        - 'uniform' : uniform weights.  All points in each neighborhood
-          are weighted equally.
-        - 'distance' : weight points by the inverse of their distance.
-          in this case, closer neighbors of a query point will have a
-          greater influence than neighbors which are further away.
-        - [callable] : a user-defined function which accepts an
-          array of distances, and returns an array of the same shape
-          containing the weights.
+    We can predict the modes of new samples
 
-        Uniform weights are used by default.
-    regressor : callable, optional ((default =
-        :func:`mean <skfda.exploratory.stats.mean>`))
-        Function to perform the local regression. By default used the mean. Can
-        accept a user-defined function wich accepts a :class:`FDataGrid` with
-        the neighbors of a test sample, and if weights != 'uniform' an array
-        of weights as second parameter.
-    algorithm : {'auto', 'ball_tree', 'brute'}, optional
-        Algorithm used to compute the nearest neighbors:
+    >>> neigh.predict(X_test).round(2) # Predict test data
+    array([ 0.39, 0.07, 0.26, 0.5 , 0.46])
 
-        - 'ball_tree' will use :class:`sklearn.neighbors.BallTree`.
-        - 'brute' will use a brute-force search.
-        - 'auto' will attempt to decide the most appropriate algorithm
-          based on the values passed to :meth:`fit` method.
 
-    leaf_size : int, optional (default = 30)
-        Leaf size passed to BallTree. This can affect the
-        speed of the construction and query, as well as the memory
-        required to store the tree. The optimal value depends on the
-        nature of the problem.
-    metric : string or callable, (default
-        :func:`lp_distance <skfda.metrics.lp_distance>`)
-        the distance metric to use for the tree.  The default metric is
-        the L2 distance. See the documentation of the metrics module
-        for a list of available metrics.
-    metric_params : dict, optional (default = None)
-        Additional keyword arguments for the metric function.
-    n_jobs : int or None, optional (default=None)
-        The number of parallel jobs to run for neighbors search.
-        ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
-        ``-1`` means using all processors.
-    sklearn_metric : boolean, optional (default = False)
-        Indicates if the metric used is a sklearn distance between vectors (see
-        :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
-        the module :mod:`skfda.misc.metrics`.
-    Examples
-    --------
-    Firstly, we will create a toy dataset with gaussian-like samples shifted,
-    and we will try to predict 5 X +1.
+    Now we will create a functional response to train the model
 
-    >>> from skfda.datasets import make_multimodal_samples
-    >>> X_train = make_multimodal_samples(n_samples=30, std=.05,
-    ...                                   random_state=0)
     >>> y_train = 5 * X_train + 1
-    >>> X_test = make_multimodal_samples(n_samples=5, std=.05, random_state=0)
+    >>> y_train
+    FDataGrid(...)
 
-    We will fit a K-Nearest Neighbors functional regressor.
+    We train the estimator with the functional response
 
-    >>> from skfda.ml.regression import KNeighborsFunctionalRegressor
-    >>> neigh = KNeighborsFunctionalRegressor()
     >>> neigh.fit(X_train, y_train)
-    KNeighborsFunctionalRegressor(algorithm='auto', leaf_size=30,...)
+    RadiusNeighborsRegressor(algorithm='auto', leaf_size=30,...)
 
-    We can predict the response of new samples.
+    And predict the responses as in the first case.
 
     >>> neigh.predict(X_test)
     FDataGrid(...)
 
     See also
     --------
-    KNeighborsClassifier
-    RadiusNeighborsClassifier
-    KNeighborsScalarRegressor
-    NearestNeighbors
-    NearestCentroids
+    :class:`~skfda.ml.classification.KNeighborsClassifier`
+    :class:`~skfda.ml.classification.RadiusNeighborsClassifier`
+    :class:`~skfda.ml.classification.NearestCentroids`
+    :class:`~skfda.ml.regression.KNeighborsRegressor`
+    :class:`~skfda.ml.clustering.NearestNeighbors`
+
     Notes
     -----
     See Nearest Neighbors in the sklearn online documentation for a discussion
@@ -368,140 +292,39 @@ class KNeighborsFunctionalRegressor(NearestNeighborsMixinInit,
 
     """
 
-    def __init__(self, n_neighbors=5, weights='uniform', regressor='mean',
+    def __init__(self, radius=1.0, weights='uniform', regressor='mean',
                  algorithm='auto', leaf_size=30, metric='l2',
-                 metric_params=None, n_jobs=1, sklearn_metric=False):
+                 metric_params=None, outlier_response=None, n_jobs=1,
+                 multivariate_metric=False):
         """Initialize the classifier."""
 
-        super().__init__(n_neighbors=n_neighbors, radius=1.,
-                         weights=weights, algorithm=algorithm,
+        super().__init__(radius=radius, weights=weights, algorithm=algorithm,
                          leaf_size=leaf_size, metric=metric,
                          metric_params=metric_params, n_jobs=n_jobs,
-                         sklearn_metric=sklearn_metric)
+                         multivariate_metric=multivariate_metric)
         self.regressor = regressor
+        self.outlier_response = outlier_response
 
-    def _query(self, X):
-        """Return distances and neighbors of given sample"""
-        return self.estimator_.kneighbors(X)
-
-
-class RadiusNeighborsFunctionalRegressor(NearestNeighborsMixinInit,
-                                         NeighborsBase, RadiusNeighborsMixin,
-                                         NeighborsFunctionalRegressorMixin):
-    """Functional regression based on neighbors within a fixed radius.
-
-    The target is predicted by local interpolation of the targets
-    associated of the nearest neighbors in the training set.
-
-    Parameters
-    ----------
-    radius : float, optional (default = 1.0)
-        Range of parameter space to use by default for :meth:`radius_neighbors`
-        queries.
-    weights : str or callable
-        weight function used in prediction.  Possible values:
-
-        - 'uniform' : uniform weights.  All points in each neighborhood
-          are weighted equally.
-        - 'distance' : weight points by the inverse of their distance.
-          in this case, closer neighbors of a query point will have a
-          greater influence than neighbors which are further away.
-        - [callable] : a user-defined function which accepts an
-          array of distances, and returns an array of the same shape
-          containing the weights.
-
-        Uniform weights are used by default.
-    regressor : callable, optional ((default =
-        :func:`mean <skfda.exploratory.stats.mean>`))
-        Function to perform the local regression. By default used the mean. Can
-        accept a user-defined function wich accepts a :class:`FDataGrid` with
-        the neighbors of a test sample, and if weights != 'uniform' an array
-        of weights as second parameter.
-    algorithm : {'auto', 'ball_tree', 'brute'}, optional
-        Algorithm used to compute the nearest neighbors:
-
-        - 'ball_tree' will use :class:`sklearn.neighbors.BallTree`.
-        - 'brute' will use a brute-force search.
-        - 'auto' will attempt to decide the most appropriate algorithm
-          based on the values passed to :meth:`fit` method.
-
-    leaf_size : int, optional (default = 30)
-        Leaf size passed to BallTree. This can affect the
-        speed of the construction and query, as well as the memory
-        required to store the tree. The optimal value depends on the
-        nature of the problem.
-    metric : string or callable, (default
-        :func:`lp_distance <skfda.metrics.lp_distance>`)
-        the distance metric to use for the tree.  The default metric is
-        the L2 distance. See the documentation of the metrics module
-        for a list of available metrics.
-    metric_params : dict, optional (default = None)
-        Additional keyword arguments for the metric function.
-    outlier_response : :class:`FDataGrid`, optional (default = None)
-        Default response for test samples without neighbors.
-    n_jobs : int or None, optional (default=None)
-        The number of parallel jobs to run for neighbors search.
-        ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
-        ``-1`` means using all processors.
-    sklearn_metric : boolean, optional (default = False)
-        Indicates if the metric used is a sklearn distance between vectors (see
-        :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
-        the module :mod:`skfda.misc.metrics`.
-    Examples
-    --------
-    Firstly, we will create a toy dataset with gaussian-like samples shifted,
-    and we will try to predict the response 5 X +1.
-
-    >>> from skfda.datasets import make_multimodal_samples
-    >>> X_train = make_multimodal_samples(n_samples=30, std=.05,
-    ...                                   random_state=0)
-    >>> y_train = 5 * X_train + 1
-    >>> X_test = make_multimodal_samples(n_samples=5, std=.05, random_state=0)
-
-    We will fit a Radius Nearest Neighbors functional regressor.
-
-    >>> from skfda.ml.regression import RadiusNeighborsFunctionalRegressor
-    >>> neigh = RadiusNeighborsFunctionalRegressor(radius=.03)
-    >>> neigh.fit(X_train, y_train)
-    RadiusNeighborsFunctionalRegressor(algorithm='auto', leaf_size=30,...)
-
-    We can predict the response of new samples.
-
-    >>> neigh.predict(X_test)
-    FDataGrid(...)
-
-    See also
-    --------
-    KNeighborsClassifier
-    RadiusNeighborsClassifier
-    KNeighborsScalarRegressor
-    NearestNeighbors
-    NearestCentroids
-    Notes
-    -----
-    See Nearest Neighbors in the sklearn online documentation for a discussion
-    of the choice of ``algorithm`` and ``leaf_size``.
-
-    This class wraps the sklearn classifier
-    `sklearn.neighbors.RadiusNeighborsClassifier`.
+    def _init_multivariate_estimator(self, sklearn_metric):
+        """Initialize the sklearn radius neighbors estimator.
 
-    https://en.wikipedia.org/wiki/K-nearest_neighbor_algorithm
+        Args:
+            sklearn_metric: (pyfunc or 'precomputed'): Metric compatible with
+                sklearn API or matrix (n_samples, n_samples) with precomputed
+                distances.
 
-    """
+        Returns:
+            Sklearn Radius Neighbors estimator initialized.
 
-    def __init__(self, radius=1., weights='uniform', regressor='mean',
-                 algorithm='auto', leaf_size=30, metric='l2',
-                 metric_params=None, outlier_response=None, n_jobs=1,
-                 sklearn_metric=False):
-        """Initialize the classifier."""
+        """
+        from sklearn.neighbors import (RadiusNeighborsRegressor as
+                                       _RadiusNeighborsRegressor)
 
-        super().__init__(n_neighbors=5, radius=radius,
-                         weights=weights, algorithm=algorithm,
-                         leaf_size=leaf_size, metric=metric,
-                         metric_params=metric_params, n_jobs=n_jobs,
-                         sklearn_metric=sklearn_metric)
-        self.regressor = regressor
-        self.outlier_response = outlier_response
+        return _RadiusNeighborsRegressor(
+            radius=self.radius, weights=self.weights,
+            algorithm=self.algorithm, leaf_size=self.leaf_size,
+            metric=sklearn_metric, metric_params=self.metric_params,
+            n_jobs=self.n_jobs)
 
     def _query(self, X):
         """Return distances and neighbors of given sample"""
diff --git a/skfda/_neighbors/unsupervised.py b/skfda/_neighbors/unsupervised.py
index 1a103cb7..9e2fbee1 100644
--- a/skfda/_neighbors/unsupervised.py
+++ b/skfda/_neighbors/unsupervised.py
@@ -1,11 +1,11 @@
 """Unsupervised learner for implementing neighbor searches."""
 
-from .base import (NearestNeighborsMixinInit, NeighborsBase, NeighborsMixin,
-                   KNeighborsMixin, RadiusNeighborsMixin, _to_sklearn_metric)
+from .base import (NeighborsBase, NeighborsMixin, KNeighborsMixin,
+                   RadiusNeighborsMixin)
 
 
-class NearestNeighbors(NearestNeighborsMixinInit, NeighborsBase,
-                       NeighborsMixin, KNeighborsMixin, RadiusNeighborsMixin):
+class NearestNeighbors(NeighborsBase, NeighborsMixin, KNeighborsMixin,
+                       RadiusNeighborsMixin):
     """Unsupervised learner for implementing neighbor searches.
 
     Parameters
@@ -40,7 +40,7 @@ class NearestNeighbors(NearestNeighborsMixinInit, NeighborsBase,
         ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
         ``-1`` means using all processors.
         Doesn't affect :meth:`fit` method.
-    sklearn_metric : boolean, optional (default = False)
+    multivariate_metric : boolean, optional (default = False)
         Indicates if the metric used is a sklearn distance between vectors (see
         :class:`sklearn.neighbors.DistanceMetric`) or a functional metric of
         the module :mod:`skfda.misc.metrics`.
@@ -82,11 +82,12 @@ class NearestNeighbors(NearestNeighborsMixinInit, NeighborsBase,
 
     See also
     --------
-    KNeighborsClassifier
-    RadiusNeighborsClassifier
-    KNeighborsScalarRegressor
-    RadiusNeighborsScalarRegressor
-    NearestCentroids
+    :class:`~skfda.ml.classification.KNeighborsClassifier`
+    :class:`~skfda.ml.classification.RadiusNeighborsClassifier`
+    :class:`~skfda.ml.classification.NearestCentroids`
+    :class:`~skfda.ml.regression.KNeighborsRegressor`
+    :class:`~skfda.ml.regression.RadiusNeighborsRegressor`
+
     Notes
     -----
     See Nearest Neighbors in the sklearn online documentation for a discussion
@@ -101,10 +102,31 @@ class NearestNeighbors(NearestNeighborsMixinInit, NeighborsBase,
 
     def __init__(self, n_neighbors=5, radius=1.0, algorithm='auto',
                  leaf_size=30, metric='l2', metric_params=None,
-                 n_jobs=1, sklearn_metric=False):
+                 n_jobs=1, multivariate_metric=False):
         """Initialize the nearest neighbors searcher."""
 
         super().__init__(n_neighbors=n_neighbors, radius=radius,
                          algorithm=algorithm, leaf_size=leaf_size,
                          metric=metric, metric_params=metric_params,
-                         n_jobs=n_jobs, sklearn_metric=sklearn_metric)
+                         n_jobs=n_jobs,
+                         multivariate_metric=multivariate_metric)
+
+    def _init_estimator(self, sklearn_metric):
+        """Initialize the sklearn nearest neighbors estimator.
+
+        Args:
+            sklearn_metric: (pyfunc or 'precomputed'): Metric compatible with
+                sklearn API or matrix (n_samples, n_samples) with precomputed
+                distances.
+
+        Returns:
+            Sklearn K Neighbors estimator initialized.
+
+        """
+        from sklearn.neighbors import NearestNeighbors as _NearestNeighbors
+
+        return _NearestNeighbors(
+            n_neighbors=self.n_neighbors, radius=self.radius,
+            algorithm=self.algorithm, leaf_size=self.leaf_size,
+            metric=sklearn_metric, metric_params=self.metric_params,
+            n_jobs=self.n_jobs)
diff --git a/skfda/ml/regression/__init__.py b/skfda/ml/regression/__init__.py
index 00f57c77..c2a67127 100644
--- a/skfda/ml/regression/__init__.py
+++ b/skfda/ml/regression/__init__.py
@@ -1,8 +1,4 @@
 
 
-from ..._neighbors import (KNeighborsScalarRegressor,
-                           RadiusNeighborsScalarRegressor,
-                           KNeighborsFunctionalRegressor,
-                           RadiusNeighborsFunctionalRegressor)
-
+from ..._neighbors import KNeighborsRegressor, RadiusNeighborsRegressor
 from .linear_model import LinearScalarRegression


Test patch:
diff --git a/tests/test_neighbors.py b/tests/test_neighbors.py
index 354251eb..d4df75fd 100644
--- a/tests/test_neighbors.py
+++ b/tests/test_neighbors.py
@@ -10,10 +10,7 @@ from skfda.ml.classification import (KNeighborsClassifier,
                                      RadiusNeighborsClassifier,
                                      NearestCentroids)
 from skfda.ml.clustering import NearestNeighbors
-from skfda.ml.regression import (KNeighborsScalarRegressor,
-                                 RadiusNeighborsScalarRegressor,
-                                 KNeighborsFunctionalRegressor,
-                                 RadiusNeighborsFunctionalRegressor)
+from skfda.ml.regression import KNeighborsRegressor, RadiusNeighborsRegressor
 from skfda.representation.basis import Fourier
 
 
@@ -22,8 +19,9 @@ class TestNeighbors(unittest.TestCase):
     def setUp(self):
         """Creates test data"""
         random_state = np.random.RandomState(0)
-        modes_location = np.concatenate((random_state.normal(-.3, .04, size=15),
-                                         random_state.normal(.3, .04, size=15)))
+        modes_location = np.concatenate(
+            (random_state.normal(-.3, .04, size=15),
+             random_state.normal(.3, .04, size=15)))
 
         idx = np.arange(30)
         random_state.shuffle(idx)
@@ -71,8 +69,8 @@ class TestNeighbors(unittest.TestCase):
 
         # Dummy test, with weight = distance, only the sample with distance 0
         # will be returned, obtaining the exact location
-        knnr = KNeighborsScalarRegressor(weights='distance')
-        rnnr = RadiusNeighborsScalarRegressor(weights='distance', radius=.1)
+        knnr = KNeighborsRegressor(weights='distance')
+        rnnr = RadiusNeighborsRegressor(weights='distance', radius=.1)
 
         knnr.fit(self.X, self.modes_location)
         rnnr.fit(self.X, self.modes_location)
@@ -83,6 +81,7 @@ class TestNeighbors(unittest.TestCase):
                                              self.modes_location)
 
     def test_kneighbors(self):
+        """Test k neighbor searches for all k-neighbors estimators"""
 
         nn = NearestNeighbors()
         nn.fit(self.X)
@@ -90,7 +89,7 @@ class TestNeighbors(unittest.TestCase):
         knn = KNeighborsClassifier()
         knn.fit(self.X, self.y)
 
-        knnr = KNeighborsScalarRegressor()
+        knnr = KNeighborsRegressor()
         knnr.fit(self.X, self.modes_location)
 
         for neigh in [nn, knn, knnr]:
@@ -120,7 +119,7 @@ class TestNeighbors(unittest.TestCase):
         knn = RadiusNeighborsClassifier(radius=.1)
         knn.fit(self.X, self.y)
 
-        knnr = RadiusNeighborsScalarRegressor(radius=.1)
+        knnr = RadiusNeighborsRegressor(radius=.1)
         knnr.fit(self.X, self.modes_location)
 
         for neigh in [nn, knn, knnr]:
@@ -143,7 +142,7 @@ class TestNeighbors(unittest.TestCase):
                 self.assertEqual(graph[0, i] == 0.0, i not in links[0])
 
     def test_knn_functional_response(self):
-        knnr = KNeighborsFunctionalRegressor(n_neighbors=1)
+        knnr = KNeighborsRegressor(n_neighbors=1)
 
         knnr.fit(self.X, self.X)
 
@@ -153,8 +152,8 @@ class TestNeighbors(unittest.TestCase):
 
     def test_knn_functional_response_sklearn(self):
         # Check sklearn metric
-        knnr = KNeighborsFunctionalRegressor(n_neighbors=1, metric='euclidean',
-                                             sklearn_metric=True)
+        knnr = KNeighborsRegressor(n_neighbors=1, metric='euclidean',
+                                             multivariate_metric=True)
         knnr.fit(self.X, self.X)
 
         res = knnr.predict(self.X)
@@ -162,7 +161,7 @@ class TestNeighbors(unittest.TestCase):
                                              self.X.data_matrix)
 
     def test_knn_functional_response_precomputed(self):
-        knnr = KNeighborsFunctionalRegressor(n_neighbors=4, weights='distance',
+        knnr = KNeighborsRegressor(n_neighbors=4, weights='distance',
                                              metric='precomputed')
         d = pairwise_distance(lp_distance)
         distances = d(self.X[:4], self.X[:4])
@@ -174,7 +173,7 @@ class TestNeighbors(unittest.TestCase):
                                              self.X[:4].data_matrix)
 
     def test_radius_functional_response(self):
-        knnr = RadiusNeighborsFunctionalRegressor(metric=lp_distance,
+        knnr = RadiusNeighborsRegressor(metric=lp_distance,
                                                   weights='distance',
                                                   regressor=l2_mean)
 
@@ -190,7 +189,7 @@ class TestNeighbors(unittest.TestCase):
 
             return np.array([w == 0 for w in weights], dtype=float)
 
-        knnr = KNeighborsFunctionalRegressor(weights=weights, n_neighbors=5)
+        knnr = KNeighborsRegressor(weights=weights, n_neighbors=5)
         response = self.X.to_basis(Fourier(domain_range=(-1, 1), nbasis=10))
         knnr.fit(self.X, response)
 
@@ -200,7 +199,7 @@ class TestNeighbors(unittest.TestCase):
 
     def test_functional_regression_distance_weights(self):
 
-        knnr = KNeighborsFunctionalRegressor(
+        knnr = KNeighborsRegressor(
             weights='distance', n_neighbors=10)
         knnr.fit(self.X[:10], self.X[:10])
         res = knnr.predict(self.X[11])
@@ -216,7 +215,7 @@ class TestNeighbors(unittest.TestCase):
                                              response.data_matrix)
 
     def test_functional_response_basis(self):
-        knnr = KNeighborsFunctionalRegressor(weights='distance', n_neighbors=5)
+        knnr = KNeighborsRegressor(weights='distance', n_neighbors=5)
         response = self.X.to_basis(Fourier(domain_range=(-1, 1), nbasis=10))
         knnr.fit(self.X, response)
 
@@ -225,7 +224,7 @@ class TestNeighbors(unittest.TestCase):
                                              response.coefficients)
 
     def test_radius_outlier_functional_response(self):
-        knnr = RadiusNeighborsFunctionalRegressor(radius=0.001)
+        knnr = RadiusNeighborsRegressor(radius=0.001)
         knnr.fit(self.X[3:6], self.X[3:6])
 
         # No value given
@@ -233,7 +232,7 @@ class TestNeighbors(unittest.TestCase):
             knnr.predict(self.X[:10])
 
         # Test response
-        knnr = RadiusNeighborsFunctionalRegressor(radius=0.001,
+        knnr = RadiusNeighborsRegressor(radius=0.001,
                                                   outlier_response=self.X[0])
         knnr.fit(self.X[:6], self.X[:6])
 
@@ -255,7 +254,7 @@ class TestNeighbors(unittest.TestCase):
 
     def test_functional_regressor_exceptions(self):
 
-        knnr = RadiusNeighborsFunctionalRegressor()
+        knnr = RadiusNeighborsRegressor()
 
         with np.testing.assert_raises(ValueError):
             knnr.fit(self.X[:3], self.X[:4])
@@ -274,7 +273,7 @@ class TestNeighbors(unittest.TestCase):
 
     def test_search_neighbors_sklearn(self):
 
-        nn = NearestNeighbors(metric='euclidean', sklearn_metric=True,
+        nn = NearestNeighbors(metric='euclidean', multivariate_metric=True,
                               n_neighbors=2)
         nn.fit(self.X[:4], self.y[:4])
 
@@ -283,9 +282,18 @@ class TestNeighbors(unittest.TestCase):
         result = np.array([[0, 3], [1, 2], [2, 1], [3, 0]])
         np.testing.assert_array_almost_equal(neighbors, result)
 
+    def test_score_scalar_response(self):
+
+        neigh = KNeighborsRegressor()
+
+        neigh.fit(self.X, self.modes_location)
+        r = neigh.score(self.X, self.modes_location)
+        np.testing.assert_almost_equal(r, 0.9975889963743335)
+
+
     def test_score_functional_response(self):
 
-        neigh = KNeighborsFunctionalRegressor()
+        neigh = KNeighborsRegressor()
 
         y = 5 * self.X + 1
         neigh.fit(self.X, y)
@@ -301,7 +309,7 @@ class TestNeighbors(unittest.TestCase):
         np.testing.assert_almost_equal(r, 0.9982527586114364)
 
     def test_score_functional_response_exceptions(self):
-        neigh = RadiusNeighborsFunctionalRegressor()
+        neigh = RadiusNeighborsRegressor()
         neigh.fit(self.X, self.X)
 
         with np.testing.assert_raises(ValueError):
@@ -309,7 +317,7 @@ class TestNeighbors(unittest.TestCase):
 
     def test_multivariate_response_score(self):
 
-        neigh = RadiusNeighborsFunctionalRegressor()
+        neigh = RadiusNeighborsRegressor()
         y = make_multimodal_samples(n_samples=5, dim_domain=2, random_state=0)
         neigh.fit(self.X[:5], y)