GH-45733: [C++][Python] Add biased/unbiased toggle to skew and kurtosis functions (#45762)

### Rationale for this change

Allow unbiased computations of skew and kurtosis. Toggle variable to compute biased/unbiased.

### What changes are included in this PR?

Add `biased` option to compute biased/unbiased versions when computing Skew and kurtosis.

### Are these changes tested?

Yes, tests added.

### Are there any user-facing changes?

There is a new `biased` option on `SkewOptions`

* GitHub Issue: #45733

Authored-by: Raúl Cumplido <[email protected]>
Signed-off-by: Antoine Pitrou <[email protected]>

Author: raulcd

cpp/src/arrow/acero/hash_aggregate_test.cc

@@ -1504,11 +1504,15 @@ TEST_P(GroupBy, VarianceOptionsAndSkewOptions) {
       /*ddof=*/0, /*skip_nulls=*/false, /*min_count=*/3);
 
   auto skew_keep_nulls = std::make_shared<SkewOptions>(/*skip_nulls=*/false,
+                                                       /*biased=*/true,
                                                        /*min_count=*/0);
-  auto skew_min_count =
-      std::make_shared<SkewOptions>(/*skip_nulls=*/true, /*min_count=*/3);
+  auto skew_min_count = std::make_shared<SkewOptions>(/*skip_nulls=*/true,
+                                                      /*biased=*/true, /*min_count=*/3);
   auto skew_keep_nulls_min_count = std::make_shared<SkewOptions>(
-      /*skip_nulls=*/false, /*min_count=*/3);
+      /*skip_nulls=*/false, /*biased=*/true, /*min_count=*/3);
+
+  auto skew_unbiased = std::make_shared<SkewOptions>(
+      /*skip_nulls=*/false, /*biased=*/false, /*min_count=*/0);
 
   for (std::string value_column : {"argument", "argument1"}) {
     for (bool use_threads : {false}) {
@@ -1553,26 +1557,30 @@ TEST_P(GroupBy, VarianceOptionsAndSkewOptions) {
                   {"hash_skew", skew_keep_nulls, value_column, "hash_skew"},
                   {"hash_skew", skew_min_count, value_column, "hash_skew"},
                   {"hash_skew", skew_keep_nulls_min_count, value_column, "hash_skew"},
+                  {"hash_skew", skew_unbiased, value_column, "hash_skew"},
                   {"hash_kurtosis", skew_keep_nulls, value_column, "hash_kurtosis"},
                   {"hash_kurtosis", skew_min_count, value_column, "hash_kurtosis"},
                   {"hash_kurtosis", skew_keep_nulls_min_count, value_column,
                    "hash_kurtosis"},
+                  {"hash_kurtosis", skew_unbiased, value_column, "hash_kurtosis"},
               },
               use_threads));
       expected = ArrayFromJSON(struct_({
                                    field("key", int64()),
                                    field("hash_skew", float64()),
                                    field("hash_skew", float64()),
                                    field("hash_skew", float64()),
+                                   field("hash_skew", float64()),
+                                   field("hash_kurtosis", float64()),
                                    field("hash_kurtosis", float64()),
                                    field("hash_kurtosis", float64()),
                                    field("hash_kurtosis", float64()),
                                }),
                                R"([
-         [1, null,      0.707106,  null,     null,      -1.5,      null     ],
-         [2, 0.213833,  0.213833,  0.213833, -1.720164, -1.720164, -1.720164],
-         [3, 0.0,       null,      null,     -2.0,       null,     null     ],
-         [4, null,      0.707106,  null,     null,      -1.5,      null     ]
+         [1, null,      0.707106,  null,     null,    null,      -1.5,      null,      null    ],
+         [2, 0.213833,  0.213833,  0.213833, 0.37037, -1.720164, -1.720164, -1.720164, -3.90123],
+         [3, 0.0,       null,      null,     null,    -2.0,       null,     null,      null    ],
+         [4, null,      0.707106,  null,     null,    null,      -1.5,      null,      null    ]
          ])");
       ValidateOutput(actual);
       AssertDatumsApproxEqual(expected, actual, /*verbose=*/true);

cpp/src/arrow/compute/api_aggregate.cc

@@ -111,6 +111,7 @@ static auto kVarianceOptionsType = GetFunctionOptionsType<VarianceOptions>(
     DataMember("min_count", &VarianceOptions::min_count));
 static auto kSkewOptionsType = GetFunctionOptionsType<SkewOptions>(
     DataMember("skip_nulls", &SkewOptions::skip_nulls),
+    DataMember("biased", &SkewOptions::biased),
     DataMember("min_count", &SkewOptions::min_count));
 static auto kQuantileOptionsType = GetFunctionOptionsType<QuantileOptions>(
     DataMember("q", &QuantileOptions::q),
@@ -154,9 +155,10 @@ VarianceOptions::VarianceOptions(int ddof, bool skip_nulls, uint32_t min_count)
       min_count(min_count) {}
 constexpr char VarianceOptions::kTypeName[];
 
-SkewOptions::SkewOptions(bool skip_nulls, uint32_t min_count)
+SkewOptions::SkewOptions(bool skip_nulls, bool biased, uint32_t min_count)
     : FunctionOptions(internal::kSkewOptionsType),
       skip_nulls(skip_nulls),
+      biased(biased),
       min_count(min_count) {}
 
 QuantileOptions::QuantileOptions(double q, enum Interpolation interpolation,

cpp/src/arrow/compute/api_aggregate.h

@@ -117,13 +117,18 @@ class ARROW_EXPORT VarianceOptions : public FunctionOptions {
 /// \brief Control Skew and Kurtosis kernel behavior
 class ARROW_EXPORT SkewOptions : public FunctionOptions {
  public:
-  explicit SkewOptions(bool skip_nulls = true, uint32_t min_count = 0);
+  explicit SkewOptions(bool skip_nulls = true, bool biased = true,
+                       uint32_t min_count = 0);
   static constexpr char const kTypeName[] = "SkewOptions";
   static SkewOptions Defaults() { return SkewOptions{}; }
 
   /// If true (the default), null values are ignored. Otherwise, if any value is null,
   /// emit null.
   bool skip_nulls;
+  /// If true (the default), the calculated value is biased. If false, the calculated
+  /// value includes a correction factor to reduce bias, making it more accurate for
+  /// small sample sizes.
+  bool biased;
   /// If less than this many non-null values are observed, emit null.
   uint32_t min_count;
 };

cpp/src/arrow/compute/kernels/aggregate_test.cc

@@ -3693,6 +3693,11 @@ TEST_F(TestSkewKurtosis, Options) {
     AssertSkewKurtosisInvalid(type, {"[]", "[]", "[]"}, options);
     AssertSkewKurtosisAre(type, "[0, 1, null, 2]", options, 0.0, -1.5);
     AssertSkewKurtosisAre(type, {"[0, 1]", "[]", "[null, 2]"}, options, 0.0, -1.5);
+    options.biased = false;
+    AssertSkewKurtosisInvalid(type, "[0, 1]", options);
+    AssertSkewKurtosisAre(type, {"[1, 2, 3]", "[40, null]"}, options, 1.9889477403978211,
+                          3.9631931024230695);
+    options.biased = true;
     options.min_count = 3;
     AssertSkewKurtosisAre(type, "[0, 1, null, 2]", options, 0.0, -1.5);
     AssertSkewKurtosisAre(type, {"[0, 1]", "[]", "[null, 2]"}, options, 0.0, -1.5);

cpp/src/arrow/compute/kernels/aggregate_var_std.cc

@@ -176,6 +176,7 @@ struct StatisticImpl : public ScalarAggregator {
       : out_type(out_type),
         stat_type(stat_type),
         skip_nulls(options.skip_nulls),
+        biased(options.biased),
         min_count(options.min_count),
         ddof(0),
         state(moments_level_for_statistic(stat_type), decimal_scale, skip_nulls) {}
@@ -197,7 +198,9 @@ struct StatisticImpl : public ScalarAggregator {
 
   Status Finalize(KernelContext*, Datum* out) override {
     if (state.count() <= ddof || state.count() < min_count ||
-        (!state.all_valid && !skip_nulls)) {
+        (!state.all_valid && !skip_nulls) ||
+        (stat_type == StatisticType::Skew && !biased && state.count() < 3) ||
+        (stat_type == StatisticType::Kurtosis && !biased && state.count() < 4)) {
       out->value = std::make_shared<DoubleScalar>();
     } else {
       switch (stat_type) {
@@ -208,10 +211,10 @@ struct StatisticImpl : public ScalarAggregator {
           out->value = std::make_shared<DoubleScalar>(state.moments.Variance(ddof));
           break;
         case StatisticType::Skew:
-          out->value = std::make_shared<DoubleScalar>(state.moments.Skew());
+          out->value = std::make_shared<DoubleScalar>(state.moments.Skew(biased));
           break;
         case StatisticType::Kurtosis:
-          out->value = std::make_shared<DoubleScalar>(state.moments.Kurtosis());
+          out->value = std::make_shared<DoubleScalar>(state.moments.Kurtosis(biased));
           break;
         default:
           return Status::NotImplemented("Unsupported statistic type ",
@@ -224,6 +227,7 @@ struct StatisticImpl : public ScalarAggregator {
   std::shared_ptr<DataType> out_type;
   StatisticType stat_type;
   bool skip_nulls;
+  bool biased;
   uint32_t min_count;
   int ddof = 0;
   MomentsState<ArrowType> state;

cpp/src/arrow/compute/kernels/aggregate_var_std_internal.h

@@ -84,14 +84,31 @@ struct Moments {
 
   double Stddev(int ddof) const { return sqrt(Variance(ddof)); }
 
-  double Skew() const {
-    // This may return NaN for m2 == 0 and m3 == 0, which is expected
-    return sqrt(count) * m3 / sqrt(m2 * m2 * m2);
+  double Skew(bool biased = true) const {
+    double result;
+    // This may return NaN for m2 == 0 and m3 == 0, which is expected.
+    if (biased) {
+      result = sqrt(count) * m3 / sqrt(m2 * m2 * m2);
+    } else {
+      auto m2_avg = m2 / count;
+      result = sqrt(count * (count - 1)) / (count - 2) * (m3 / count) /
+               sqrt(m2_avg * m2_avg * m2_avg);
+    }
+    return result;
   }
 
-  double Kurtosis() const {
-    // This may return NaN for m2 == 0 and m4 == 0, which is expected
-    return count * m4 / (m2 * m2) - 3;
+  double Kurtosis(bool biased = true) const {
+    double result;
+    // This may return NaN for m2 == 0 and m4 == 0, which is expected.
+    if (biased) {
+      result = count * m4 / (m2 * m2) - 3;
+    } else {
+      auto m2_avg = m2 / count;
+      result = 1.0 / ((count - 2) * (count - 3)) *
+               (((count * count) - 1.0) * (m4 / count) / (m2_avg * m2_avg) -
+                3 * ((count - 1) * (count - 1)));
+    }
+    return result;
   }
 
   void MergeFrom(int level, const Moments& other) { *this = Merge(level, *this, other); }

cpp/src/arrow/compute/kernels/hash_aggregate_numeric.cc

@@ -452,35 +452,37 @@ struct GroupedStatisticImpl : public GroupedAggregator {
   Status InitInternal(ExecContext* ctx, const KernelInitArgs& args,
                       StatisticType stat_type, const VarianceOptions& options) {
     return InitInternal(ctx, args, stat_type, options.ddof, options.skip_nulls,
-                        options.min_count);
+                        /*biased=*/false, options.min_count);
   }
 
   // Init helper for hash_skew and hash_kurtosis
   Status InitInternal(ExecContext* ctx, const KernelInitArgs& args,
                       StatisticType stat_type, const SkewOptions& options) {
     return InitInternal(ctx, args, stat_type, /*ddof=*/0, options.skip_nulls,
-                        options.min_count);
+                        options.biased, options.min_count);
   }
 
   Status InitInternal(ExecContext* ctx, const KernelInitArgs& args,
-                      StatisticType stat_type, int ddof, bool skip_nulls,
+                      StatisticType stat_type, int ddof, bool skip_nulls, bool biased,
                       uint32_t min_count) {
     if constexpr (is_decimal_type<Type>::value) {
       int32_t decimal_scale =
           checked_cast<const DecimalType&>(*args.inputs[0].type).scale();
-      return InitInternal(ctx, stat_type, decimal_scale, ddof, skip_nulls, min_count);
+      return InitInternal(ctx, stat_type, decimal_scale, ddof, skip_nulls, biased,
+                          min_count);
     } else {
-      return InitInternal(ctx, stat_type, /*decimal_scale=*/0, ddof, skip_nulls,
+      return InitInternal(ctx, stat_type, /*decimal_scale=*/0, ddof, skip_nulls, biased,
                           min_count);
     }
   }
 
   Status InitInternal(ExecContext* ctx, StatisticType stat_type, int32_t decimal_scale,
-                      int ddof, bool skip_nulls, uint32_t min_count) {
+                      int ddof, bool skip_nulls, bool biased, uint32_t min_count) {
     stat_type_ = stat_type;
     moments_level_ = moments_level_for_statistic(stat_type_);
     decimal_scale_ = decimal_scale;
     skip_nulls_ = skip_nulls;
+    biased_ = biased;
     min_count_ = min_count;
     ddof_ = ddof;
     ctx_ = ctx;
@@ -539,7 +541,7 @@ struct GroupedStatisticImpl : public GroupedAggregator {
   Status ConsumeGeneric(const ExecSpan& batch) {
     GroupedStatisticImpl<Type> state;
     RETURN_NOT_OK(state.InitInternal(ctx_, stat_type_, decimal_scale_, ddof_, skip_nulls_,
-                                     min_count_));
+                                     biased_, min_count_));
     RETURN_NOT_OK(state.Resize(num_groups_));
     int64_t* counts = state.counts_.mutable_data();
     double* means = state.means_.mutable_data();
@@ -612,7 +614,7 @@ struct GroupedStatisticImpl : public GroupedAggregator {
       var_std.resize(num_groups_);
       GroupedStatisticImpl<Type> state;
       RETURN_NOT_OK(state.InitInternal(ctx_, stat_type_, decimal_scale_, ddof_,
-                                       skip_nulls_, min_count_));
+                                       skip_nulls_, biased_, min_count_));
       RETURN_NOT_OK(state.Resize(num_groups_));
       int64_t* other_counts = state.counts_.mutable_data();
       double* other_means = state.means_.mutable_data();
@@ -739,7 +741,9 @@ struct GroupedStatisticImpl : public GroupedAggregator {
     const double* m3s = m3s_data();
     const double* m4s = m4s_data();
     for (int64_t i = 0; i < num_groups_; ++i) {
-      if (counts[i] > ddof_ && counts[i] >= min_count_) {
+      if (counts[i] > ddof_ && counts[i] >= min_count_ &&
+          (stat_type_ != StatisticType::Skew || biased_ || counts[i] > 2) &&
+          (stat_type_ != StatisticType::Kurtosis || biased_ || counts[i] > 3)) {
         const auto moments = Moments(counts[i], means[i], m2s[i], m3s[i], m4s[i]);
         switch (stat_type_) {
           case StatisticType::Var:
@@ -749,10 +753,10 @@ struct GroupedStatisticImpl : public GroupedAggregator {
             results[i] = moments.Stddev(ddof_);
             break;
           case StatisticType::Skew:
-            results[i] = moments.Skew();
+            results[i] = moments.Skew(biased_);
             break;
           case StatisticType::Kurtosis:
-            results[i] = moments.Kurtosis();
+            results[i] = moments.Kurtosis(biased_);
             break;
           default:
             return Status::NotImplemented("Statistic type ",
@@ -809,6 +813,7 @@ struct GroupedStatisticImpl : public GroupedAggregator {
   int moments_level_;
   int32_t decimal_scale_;
   bool skip_nulls_;
+  bool biased_;
   uint32_t min_count_;
   int ddof_;
   int64_t num_groups_ = 0;

python/pyarrow/_compute.pyx

@@ -1909,8 +1909,8 @@ class VarianceOptions(_VarianceOptions):
 
 
 cdef class _SkewOptions(FunctionOptions):
-    def _set_options(self, skip_nulls, min_count):
-        self.wrapped.reset(new CSkewOptions(skip_nulls, min_count))
+    def _set_options(self, skip_nulls, biased, min_count):
+        self.wrapped.reset(new CSkewOptions(skip_nulls, biased, min_count))
 
 
 class SkewOptions(_SkewOptions):
@@ -1920,11 +1920,14 @@ class SkewOptions(_SkewOptions):
     Parameters
     ----------
     {_skip_nulls_doc()}
+    biased : bool, default True
+        Whether the calculated value is biased.
+        If False, the value computed includes a correction factor to reduce bias.
     {_min_count_doc(default=0)}
     """
 
-    def __init__(self, *, skip_nulls=True, min_count=0):
-        self._set_options(skip_nulls, min_count)
+    def __init__(self, *, skip_nulls=True, biased=True, min_count=0):
+        self._set_options(skip_nulls, biased, min_count)
 
 
 cdef class _SplitOptions(FunctionOptions):

python/pyarrow/includes/libarrow.pxd

@@ -2628,8 +2628,9 @@ cdef extern from "arrow/compute/api.h" namespace "arrow::compute" nogil:
 
     cdef cppclass CSkewOptions \
             "arrow::compute::SkewOptions"(CFunctionOptions):
-        CSkewOptions(c_bool skip_nulls, uint32_t min_count)
+        CSkewOptions(c_bool skip_nulls, c_bool biased, uint32_t min_count)
         c_bool skip_nulls
+        c_bool biased
         uint32_t min_count
 
     cdef cppclass CScalarAggregateOptions \

python/pyarrow/tests/test_compute.py

@@ -457,6 +457,21 @@ def test_kurtosis():
     assert pc.kurtosis(data, min_count=4).as_py() is None
 
 
+@pytest.mark.parametrize("input, expected", (
+    (
+        [1.0, 2.0, 3.0, 40.0, None],
+        {'skew': 1.988947740397821, 'kurtosis': 3.9631931024230695}
+    ),
+    ([1, 2, 40], {'skew': 1.7281098503730385, 'kurtosis': None}),
+    ([1, 40], {'skew': None, 'kurtosis': None}),
+))
+def test_unbiased_skew_and_kurtosis(input, expected):
+    arrow_skew = pc.skew(input, skip_nulls=True, biased=False)
+    arrow_kurtosis = pc.kurtosis(input, skip_nulls=True, biased=False)
+    assert arrow_skew == pa.scalar(expected['skew'], type=pa.float64())
+    assert arrow_kurtosis == pa.scalar(expected['kurtosis'], type=pa.float64())
+
+
 def test_count_substring():
     for (ty, offset) in [(pa.string(), pa.int32()),
                          (pa.large_string(), pa.int64())]: