scikit-learn · boricles · Sep 13, 2020 · Sep 13, 2020 · Sep 14, 2020 · Sep 14, 2020
diff --git a/sklearn/gaussian_process/_gpr.py b/sklearn/gaussian_process/_gpr.py
@@ -196,10 +196,15 @@ def fit(self, X, y):
        if self.normalize_y:
            self._y_train_mean = np.mean(y, axis=0)
            self._y_train_std = np.std(y, axis=0)
-
+            # Assign a standard deviation of one to constant
+            # targets for avoiding division by zero errors
+            if np.iterable(self._y_train_std):
+                self._y_train_std[self._y_train_std == 0] = 1
+            else:
+                self._y_train_std = (
+                    self._y_train_std if self._y_train_std != 0 else 1)
            # Remove mean and make unit variance
            y = (y - self._y_train_mean) / self._y_train_std
-
        else:
            self._y_train_mean = np.zeros(1)
            self._y_train_std = 1

diff --git a/sklearn/gaussian_process/tests/test_gpr.py b/sklearn/gaussian_process/tests/test_gpr.py
@@ -536,3 +536,17 @@ def test_bound_check_fixed_hyperparameter():
                        periodicity_bounds="fixed")  # seasonal component
    kernel = k1 + k2
    GaussianProcessRegressor(kernel=kernel).fit(X, y)
+
+
+@pytest.mark.parametrize('kernel', kernels)
+def test_constant_target(kernel):
+    # Test for constant target(s), the private _y_train_std attribute has
+    # a standard deviation of 1 instead of 0
+    gpr = GaussianProcessRegressor(kernel=kernel, normalize_y=True)
+    # create a constant targets of 2
+    y = np.full((1, X.shape[0]), 2).ravel()
+    expected_std = np.full((1, X.shape[0]), 1).ravel()
+    gpr.fit(X, y)
+    _y_train_std = gpr._y_train_std
+
+    assert_array_equal(_y_train_std, expected_std)