From 23b4165a6b603bf98062fd32b5055f9b701971f6 Mon Sep 17 00:00:00 2001
From: Jiaming Yuan <jm.yuan@outlook.com>
Date: Sat, 20 Mar 2021 01:56:17 +0800
Subject: [PATCH] Fix gamma deviance (#6761)

---
 src/metric/elementwise_metric.cu  | 25 +++++++++++++++++--------
 tests/python/test_eval_metrics.py | 20 ++++++++++++++++++++
 2 files changed, 37 insertions(+), 8 deletions(-)

diff --git a/src/metric/elementwise_metric.cu b/src/metric/elementwise_metric.cu
index f106fc547..b8322d647 100644
--- a/src/metric/elementwise_metric.cu
+++ b/src/metric/elementwise_metric.cu
@@ -274,18 +274,27 @@ struct EvalPoissonNegLogLik {
   }
 };
 
+/**
+ * Gamma deviance
+ *
+ *   Expected input:
+ *   label >= 0
+ *   predt >= 0
+ */
 struct EvalGammaDeviance {
-  const char *Name() const {
-    return "gamma-deviance";
+  const char *Name() const { return "gamma-deviance"; }
+
+  XGBOOST_DEVICE bst_float EvalRow(bst_float label, bst_float predt) const {
+    predt += kRtEps;
+    label += kRtEps;
+    return std::log(predt / label) + label / predt - 1;
   }
 
-  XGBOOST_DEVICE bst_float EvalRow(bst_float label, bst_float pred) const {
-    bst_float epsilon = 1.0e-9;
-    bst_float tmp = label / (pred + epsilon);
-    return tmp - std::log(tmp) - 1;
-  }
   static bst_float GetFinal(bst_float esum, bst_float wsum) {
-    return 2 * esum;
+    if (wsum <= 0) {
+      wsum = kRtEps;
+    }
+    return 2 * esum / wsum;
   }
 };
 
diff --git a/tests/python/test_eval_metrics.py b/tests/python/test_eval_metrics.py
index 53d42b1a0..411749124 100644
--- a/tests/python/test_eval_metrics.py
+++ b/tests/python/test_eval_metrics.py
@@ -103,3 +103,23 @@ class TestEvalMetrics:
         assert gbdt_01.predict(dvalid)[0] == gbdt_02.predict(dvalid)[0]
         assert gbdt_01.predict(dvalid)[0] == gbdt_03.predict(dvalid)[0]
         assert gbdt_03.predict(dvalid)[0] != gbdt_04.predict(dvalid)[0]
+
+    @pytest.mark.skipif(**tm.no_sklearn())
+    def test_gamma_deviance(self):
+        from sklearn.metrics import mean_gamma_deviance
+        rng = np.random.RandomState(1994)
+        n_samples = 100
+        n_features = 30
+
+        X = rng.randn(n_samples, n_features)
+        y = rng.randn(n_samples)
+        y = y - y.min() * 100
+
+        reg = xgb.XGBRegressor(tree_method="hist", objective="reg:gamma", n_estimators=10)
+        reg.fit(X, y, eval_metric="gamma-deviance")
+
+        booster = reg.get_booster()
+        score = reg.predict(X)
+        gamma_dev = float(booster.eval(xgb.DMatrix(X, y)).split(":")[1].split(":")[0])
+        skl_gamma_dev = mean_gamma_deviance(y, score)
+        np.testing.assert_allclose(gamma_dev, skl_gamma_dev, rtol=1e-6)