You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@arrow.apache.org by ap...@apache.org on 2025/03/07 20:17:56 UTC
(arrow) branch main updated: GH-45676: [C++][Python][Compute] Add skew and kurtosis functions (#45677)
This is an automated email from the ASF dual-hosted git repository.
apitrou pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/arrow.git
The following commit(s) were added to refs/heads/main by this push:
new 13c898cbe9 GH-45676: [C++][Python][Compute] Add skew and kurtosis functions (#45677)
13c898cbe9 is described below
commit 13c898cbe9ef08daf1bc7bdf08280c5bfa16086d
Author: Antoine Pitrou <an...@python.org>
AuthorDate: Fri Mar 7 21:17:50 2025 +0100
GH-45676: [C++][Python][Compute] Add skew and kurtosis functions (#45677)
### What changes are included in this PR?
* Add "skew" and "kurtosis" scalar aggregate functions
* Add "hash_skew" and "hash_kurtosis" grouped aggregated functions
* Add support for decimal32/decimal64 to "variance", "stddev", "hash_variance" and "hash_stddev"
* Small optimization for "hash_variance" and "hash_stddev": avoid allocating the identity group id mapping when consuming an ExecBatch
### Are these changes tested?
Yes, by included unit tests.
### Are there any user-facing changes?
No, except new APIs.
* GitHub Issue: #45676
Authored-by: Antoine Pitrou <an...@python.org>
Signed-off-by: Antoine Pitrou <an...@python.org>
---
cpp/src/arrow/acero/hash_aggregate_test.cc | 379 ++++++++++---------
cpp/src/arrow/compute/api_aggregate.cc | 17 +
cpp/src/arrow/compute/api_aggregate.h | 42 +++
cpp/src/arrow/compute/kernels/aggregate_internal.h | 11 +-
cpp/src/arrow/compute/kernels/aggregate_test.cc | 174 ++++++++-
cpp/src/arrow/compute/kernels/aggregate_var_std.cc | 374 ++++++++++---------
.../compute/kernels/aggregate_var_std_internal.h | 110 +++++-
cpp/src/arrow/compute/kernels/hash_aggregate.cc | 409 ++++++++++++++-------
cpp/src/arrow/testing/gtest_util_test.cc | 10 +
cpp/src/arrow/testing/math.cc | 6 +
cpp/src/arrow/util/math_internal.h | 26 ++
cpp/src/arrow/util/math_test.cc | 33 ++
docs/source/cpp/compute.rst | 8 +
python/pyarrow/_compute.pyx | 19 +
python/pyarrow/compute.py | 1 +
python/pyarrow/includes/libarrow.pxd | 6 +
python/pyarrow/tests/test_compute.py | 15 +
17 files changed, 1120 insertions(+), 520 deletions(-)
diff --git a/cpp/src/arrow/acero/hash_aggregate_test.cc b/cpp/src/arrow/acero/hash_aggregate_test.cc
index e50113968e..347bb96269 100644
--- a/cpp/src/arrow/acero/hash_aggregate_test.cc
+++ b/cpp/src/arrow/acero/hash_aggregate_test.cc
@@ -77,6 +77,7 @@ using compute::PivotWiderOptions;
using compute::RowSegmenter;
using compute::ScalarAggregateOptions;
using compute::Segment;
+using compute::SkewOptions;
using compute::SortIndices;
using compute::SortKey;
using compute::SortOrder;
@@ -1080,9 +1081,67 @@ TEST_P(GroupBy, MeanOverflow) {
}
}
-TEST_P(GroupBy, VarianceAndStddev) {
+TEST_P(GroupBy, VarianceStddevSkewKurtosis) {
+ for (auto value_type : {int32(), float64()}) {
+ ARROW_SCOPED_TRACE("value_type = ", *value_type);
+ auto batch = RecordBatchFromJSON(
+ schema({field("argument", value_type), field("key", int64())}), R"([
+ [1, 1],
+ [null, 1],
+ [0, 2],
+ [null, 3],
+ [4, null],
+ [3, 1],
+ [0, 2],
+ [-1, 2],
+ [1, null],
+ [null, 3]
+ ])");
+
+ ASSERT_OK_AND_ASSIGN(Datum aggregated_and_grouped,
+ GroupByTest(
+ {
+ batch->GetColumnByName("argument"),
+ batch->GetColumnByName("argument"),
+ batch->GetColumnByName("argument"),
+ batch->GetColumnByName("argument"),
+ },
+ {
+ batch->GetColumnByName("key"),
+ },
+ {},
+ {
+ {"hash_variance", nullptr},
+ {"hash_stddev", nullptr},
+ {"hash_skew", nullptr},
+ {"hash_kurtosis", nullptr},
+ },
+ false));
+
+ auto expected = ArrayFromJSON(struct_({
+ field("key_0", int64()),
+ field("hash_variance", float64()),
+ field("hash_stddev", float64()),
+ field("hash_skew", float64()),
+ field("hash_kurtosis", float64()),
+ }),
+ R"([
+ [1, 1.0, 1.0, 0.0, -2.0],
+ [2, 0.22222222222222224, 0.4714045207910317, -0.7071067811865478, -1.5],
+ [3, null, null, null, null],
+ [null, 2.25, 1.5, 0.0, -2.0]
+ ])");
+ AssertDatumsApproxEqual(expected, aggregated_and_grouped,
+ /*verbose=*/true);
+ }
+}
+
+TEST_P(GroupBy, VarianceAndStddevDdof) {
+ // Test ddof
+ auto variance_options = std::make_shared<VarianceOptions>(/*ddof=*/2);
+
auto batch = RecordBatchFromJSON(
- schema({field("argument", int32()), field("key", int64())}), R"([
+ schema({field("argument", float64()), field("key", int64())}), R"([
[1, 1],
[null, 1],
[0, 2],
@@ -1094,83 +1153,7 @@ TEST_P(GroupBy, VarianceAndStddev) {
[1, null],
[null, 3]
])");
-
ASSERT_OK_AND_ASSIGN(Datum aggregated_and_grouped,
- GroupByTest(
- {
- batch->GetColumnByName("argument"),
- batch->GetColumnByName("argument"),
- },
- {
- batch->GetColumnByName("key"),
- },
- {},
- {
- {"hash_variance", nullptr},
- {"hash_stddev", nullptr},
- },
- false));
-
- AssertDatumsApproxEqual(ArrayFromJSON(struct_({
- field("key_0", int64()),
- field("hash_variance", float64()),
- field("hash_stddev", float64()),
- }),
- R"([
- [1, 1.0, 1.0 ],
- [2, 0.22222222222222224, 0.4714045207910317],
- [3, null, null ],
- [null, 2.25, 1.5 ]
- ])"),
- aggregated_and_grouped,
- /*verbose=*/true);
-
- batch = RecordBatchFromJSON(
- schema({field("argument", float64()), field("key", int64())}), R"([
- [1.0, 1],
- [null, 1],
- [0.0, 2],
- [null, 3],
- [4.0, null],
- [3.0, 1],
- [0.0, 2],
- [-1.0, 2],
- [1.0, null],
- [null, 3]
- ])");
-
- ASSERT_OK_AND_ASSIGN(aggregated_and_grouped, GroupByTest(
- {
- batch->GetColumnByName("argument"),
- batch->GetColumnByName("argument"),
- },
- {
- batch->GetColumnByName("key"),
- },
- {},
- {
- {"hash_variance", nullptr},
- {"hash_stddev", nullptr},
- },
- false));
-
- AssertDatumsApproxEqual(ArrayFromJSON(struct_({
- field("key_0", int64()),
- field("hash_variance", float64()),
- field("hash_stddev", float64()),
- }),
- R"([
- [1, 1.0, 1.0 ],
- [2, 0.22222222222222224, 0.4714045207910317],
- [3, null, null ],
- [null, 2.25, 1.5 ]
- ])"),
- aggregated_and_grouped,
- /*verbose=*/true);
-
- // Test ddof
- auto variance_options = std::make_shared<VarianceOptions>(/*ddof=*/2);
- ASSERT_OK_AND_ASSIGN(aggregated_and_grouped,
GroupByTest(
{
batch->GetColumnByName("argument"),
@@ -1201,55 +1184,59 @@ TEST_P(GroupBy, VarianceAndStddev) {
/*verbose=*/true);
}
-TEST_P(GroupBy, VarianceAndStddevDecimal) {
- auto batch = RecordBatchFromJSON(
- schema({field("argument0", decimal128(3, 2)), field("argument1", decimal128(3, 2)),
- field("key", int64())}),
- R"([
- ["1.00", "1.00", 1],
- [null, null, 1],
- ["0.00", "0.00", 2],
- ["4.00", "4.00", null],
- ["3.00", "3.00", 1],
- ["0.00", "0.00", 2],
- ["-1.00", "-1.00", 2],
- ["1.00", "1.00", null]
- ])");
+TEST_P(GroupBy, VarianceStddevSkewKurtosisDecimal) {
+ for (auto value_type :
+ {decimal32(3, 2), decimal64(3, 2), decimal128(3, 2), decimal256(3, 2)}) {
+ ARROW_SCOPED_TRACE("value_type = ", *value_type);
+ auto batch = RecordBatchFromJSON(
+ schema({field("argument", value_type), field("key", int64())}),
+ R"([
+ ["1.00", 1],
+ [null, 1],
+ ["0.00", 2],
+ ["4.00", null],
+ ["3.00", 1],
+ ["0.00", 2],
+ ["-1.00", 2],
+ ["1.00", null]
+ ])");
- ASSERT_OK_AND_ASSIGN(Datum aggregated_and_grouped,
- GroupByTest(
- {
- batch->GetColumnByName("argument0"),
- batch->GetColumnByName("argument0"),
- batch->GetColumnByName("argument1"),
- batch->GetColumnByName("argument1"),
- },
- {
- batch->GetColumnByName("key"),
- },
- {},
- {
- {"hash_variance", nullptr},
- {"hash_stddev", nullptr},
- {"hash_variance", nullptr},
- {"hash_stddev", nullptr},
- },
- false));
+ ASSERT_OK_AND_ASSIGN(Datum aggregated_and_grouped,
+ GroupByTest(
+ {
+ batch->GetColumnByName("argument"),
+ batch->GetColumnByName("argument"),
+ batch->GetColumnByName("argument"),
+ batch->GetColumnByName("argument"),
+ },
+ {
+ batch->GetColumnByName("key"),
+ },
+ {},
+ {
+ {"hash_variance", nullptr},
+ {"hash_stddev", nullptr},
+ {"hash_skew", nullptr},
+ {"hash_kurtosis", nullptr},
+ },
+ false));
- AssertDatumsApproxEqual(ArrayFromJSON(struct_({
- field("key_0", int64()),
- field("hash_variance", float64()),
- field("hash_stddev", float64()),
- field("hash_variance", float64()),
- field("hash_stddev", float64()),
- }),
- R"([
- [1, 1.0, 1.0, 1.0, 1.0 ],
- [2, 0.22222222222222224, 0.4714045207910317, 0.22222222222222224, 0.4714045207910317],
- [null, 2.25, 1.5, 2.25, 1.5 ]
- ])"),
- aggregated_and_grouped,
- /*verbose=*/true);
+ auto expected = ArrayFromJSON(struct_({
+ field("key_0", int64()),
+ field("hash_variance", float64()),
+ field("hash_stddev", float64()),
+ field("hash_skew", float64()),
+ field("hash_kurtosis", float64()),
+ }),
+ R"([
+ [1, 1.0, 1.0, 0.0, -2.0],
+ [2, 0.22222222222222224, 0.4714045207910317, -0.7071067811865478, -1.5],
+ [null, 2.25, 1.5, 0.0, -2.0]
+ ])");
+
+ AssertDatumsApproxEqual(expected, aggregated_and_grouped,
+ /*verbose=*/true);
+ }
}
TEST_P(GroupBy, TDigest) {
@@ -1487,7 +1474,7 @@ TEST_P(GroupBy, StddevVarianceTDigestScalar) {
}
}
-TEST_P(GroupBy, VarianceOptions) {
+TEST_P(GroupBy, VarianceOptionsAndSkewOptions) {
BatchesWithSchema input;
input.batches = {
ExecBatchFromJSON(
@@ -1503,81 +1490,93 @@ TEST_P(GroupBy, VarianceOptions) {
"[[null, null, 1]]"),
ExecBatchFromJSON({int32(), float32(), int64()}, "[[2, 2.0, 1], [3, 3.0, 2]]"),
ExecBatchFromJSON({int32(), float32(), int64()}, "[[4, 4.0, 2], [2, 2.0, 4]]"),
- ExecBatchFromJSON({int32(), float32(), int64()}, "[[null, null, 4]]"),
+ ExecBatchFromJSON({int32(), float32(), int64()}, "[[null, null, 4], [6, 6.0, 3]]"),
};
input.schema = schema(
{field("argument", int32()), field("argument1", float32()), field("key", int64())});
- auto keep_nulls = std::make_shared<VarianceOptions>(/*ddof=*/0, /*skip_nulls=*/false,
- /*min_count=*/0);
- auto min_count =
+ auto var_keep_nulls =
+ std::make_shared<VarianceOptions>(/*ddof=*/0, /*skip_nulls=*/false,
+ /*min_count=*/0);
+ auto var_min_count =
std::make_shared<VarianceOptions>(/*ddof=*/0, /*skip_nulls=*/true, /*min_count=*/3);
- auto keep_nulls_min_count = std::make_shared<VarianceOptions>(
+ auto var_keep_nulls_min_count = std::make_shared<VarianceOptions>(
/*ddof=*/0, /*skip_nulls=*/false, /*min_count=*/3);
- for (bool use_threads : {false}) {
- SCOPED_TRACE(use_threads ? "parallel/merged" : "serial");
- ASSERT_OK_AND_ASSIGN(
- Datum actual,
- RunGroupBy(
- input, {"key"},
- {
- {"hash_stddev", keep_nulls, "argument", "hash_stddev"},
- {"hash_stddev", min_count, "argument", "hash_stddev"},
- {"hash_stddev", keep_nulls_min_count, "argument", "hash_stddev"},
- {"hash_variance", keep_nulls, "argument", "hash_variance"},
- {"hash_variance", min_count, "argument", "hash_variance"},
- {"hash_variance", keep_nulls_min_count, "argument", "hash_variance"},
- },
- use_threads));
- Datum expected = ArrayFromJSON(struct_({
- field("key", int64()),
- field("hash_stddev", float64()),
- field("hash_stddev", float64()),
- field("hash_stddev", float64()),
- field("hash_variance", float64()),
- field("hash_variance", float64()),
- field("hash_variance", float64()),
- }),
- R"([
- [1, null, 0.471405, null, null, 0.222222, null ],
- [2, 1.29904, 1.29904, 1.29904, 1.6875, 1.6875, 1.6875],
- [3, 0.0, null, null, 0.0, null, null ],
- [4, null, 0.471405, null, null, 0.222222, null ]
- ])");
- ValidateOutput(expected);
- AssertDatumsApproxEqual(expected, actual, /*verbose=*/true);
+ auto skew_keep_nulls = std::make_shared<SkewOptions>(/*skip_nulls=*/false,
+ /*min_count=*/0);
+ auto skew_min_count =
+ std::make_shared<SkewOptions>(/*skip_nulls=*/true, /*min_count=*/3);
+ auto skew_keep_nulls_min_count = std::make_shared<SkewOptions>(
+ /*skip_nulls=*/false, /*min_count=*/3);
- ASSERT_OK_AND_ASSIGN(
- actual,
- RunGroupBy(
- input, {"key"},
- {
- {"hash_stddev", keep_nulls, "argument1", "hash_stddev"},
- {"hash_stddev", min_count, "argument1", "hash_stddev"},
- {"hash_stddev", keep_nulls_min_count, "argument1", "hash_stddev"},
- {"hash_variance", keep_nulls, "argument1", "hash_variance"},
- {"hash_variance", min_count, "argument1", "hash_variance"},
- {"hash_variance", keep_nulls_min_count, "argument1", "hash_variance"},
- },
- use_threads));
- expected = ArrayFromJSON(struct_({
- field("key", int64()),
- field("hash_stddev", float64()),
- field("hash_stddev", float64()),
- field("hash_stddev", float64()),
- field("hash_variance", float64()),
- field("hash_variance", float64()),
- field("hash_variance", float64()),
- }),
- R"([
+ for (std::string value_column : {"argument", "argument1"}) {
+ for (bool use_threads : {false}) {
+ SCOPED_TRACE(use_threads ? "parallel/merged" : "serial");
+ ASSERT_OK_AND_ASSIGN(
+ Datum actual,
+ RunGroupBy(
+ input, {"key"},
+ {
+ {"hash_stddev", var_keep_nulls, value_column, "hash_stddev"},
+ {"hash_stddev", var_min_count, value_column, "hash_stddev"},
+ {"hash_stddev", var_keep_nulls_min_count, value_column, "hash_stddev"},
+ {"hash_variance", var_keep_nulls, value_column, "hash_variance"},
+ {"hash_variance", var_min_count, value_column, "hash_variance"},
+ {"hash_variance", var_keep_nulls_min_count, value_column,
+ "hash_variance"},
+ },
+ use_threads));
+ Datum expected = ArrayFromJSON(struct_({
+ field("key", int64()),
+ field("hash_stddev", float64()),
+ field("hash_stddev", float64()),
+ field("hash_stddev", float64()),
+ field("hash_variance", float64()),
+ field("hash_variance", float64()),
+ field("hash_variance", float64()),
+ }),
+ R"([
[1, null, 0.471405, null, null, 0.222222, null ],
[2, 1.29904, 1.29904, 1.29904, 1.6875, 1.6875, 1.6875],
- [3, 0.0, null, null, 0.0, null, null ],
+ [3, 2.5, null, null, 6.25, null, null ],
[4, null, 0.471405, null, null, 0.222222, null ]
- ])");
- ValidateOutput(expected);
- AssertDatumsApproxEqual(expected, actual, /*verbose=*/true);
+ ])");
+ ValidateOutput(actual);
+ AssertDatumsApproxEqual(expected, actual, /*verbose=*/true);
+
+ ASSERT_OK_AND_ASSIGN(
+ actual,
+ RunGroupBy(
+ input, {"key"},
+ {
+ {"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_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"},
+ },
+ use_threads));
+ expected = ArrayFromJSON(struct_({
+ field("key", int64()),
+ field("hash_skew", float64()),
+ field("hash_skew", float64()),
+ field("hash_skew", 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 ]
+ ])");
+ ValidateOutput(actual);
+ AssertDatumsApproxEqual(expected, actual, /*verbose=*/true);
+ }
}
}
diff --git a/cpp/src/arrow/compute/api_aggregate.cc b/cpp/src/arrow/compute/api_aggregate.cc
index 37119d2b67..20d3ce2faf 100644
--- a/cpp/src/arrow/compute/api_aggregate.cc
+++ b/cpp/src/arrow/compute/api_aggregate.cc
@@ -109,6 +109,9 @@ static auto kVarianceOptionsType = GetFunctionOptionsType<VarianceOptions>(
DataMember("ddof", &VarianceOptions::ddof),
DataMember("skip_nulls", &VarianceOptions::skip_nulls),
DataMember("min_count", &VarianceOptions::min_count));
+static auto kSkewOptionsType = GetFunctionOptionsType<SkewOptions>(
+ DataMember("skip_nulls", &SkewOptions::skip_nulls),
+ DataMember("min_count", &SkewOptions::min_count));
static auto kQuantileOptionsType = GetFunctionOptionsType<QuantileOptions>(
DataMember("q", &QuantileOptions::q),
DataMember("interpolation", &QuantileOptions::interpolation),
@@ -151,6 +154,11 @@ 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)
+ : FunctionOptions(internal::kSkewOptionsType),
+ skip_nulls(skip_nulls),
+ min_count(min_count) {}
+
QuantileOptions::QuantileOptions(double q, enum Interpolation interpolation,
bool skip_nulls, uint32_t min_count)
: FunctionOptions(internal::kQuantileOptionsType),
@@ -203,6 +211,7 @@ void RegisterAggregateOptions(FunctionRegistry* registry) {
DCHECK_OK(registry->AddFunctionOptionsType(kCountOptionsType));
DCHECK_OK(registry->AddFunctionOptionsType(kModeOptionsType));
DCHECK_OK(registry->AddFunctionOptionsType(kVarianceOptionsType));
+ DCHECK_OK(registry->AddFunctionOptionsType(kSkewOptionsType));
DCHECK_OK(registry->AddFunctionOptionsType(kQuantileOptionsType));
DCHECK_OK(registry->AddFunctionOptionsType(kTDigestOptionsType));
DCHECK_OK(registry->AddFunctionOptionsType(kPivotOptionsType));
@@ -271,6 +280,14 @@ Result<Datum> Variance(const Datum& value, const VarianceOptions& options,
return CallFunction("variance", {value}, &options, ctx);
}
+Result<Datum> Skew(const Datum& value, const SkewOptions& options, ExecContext* ctx) {
+ return CallFunction("skew", {value}, &options, ctx);
+}
+
+Result<Datum> Kurtosis(const Datum& value, const SkewOptions& options, ExecContext* ctx) {
+ return CallFunction("kurtosis", {value}, &options, ctx);
+}
+
Result<Datum> Quantile(const Datum& value, const QuantileOptions& options,
ExecContext* ctx) {
return CallFunction("quantile", {value}, &options, ctx);
diff --git a/cpp/src/arrow/compute/api_aggregate.h b/cpp/src/arrow/compute/api_aggregate.h
index 66876ee66e..61bab4cdb8 100644
--- a/cpp/src/arrow/compute/api_aggregate.h
+++ b/cpp/src/arrow/compute/api_aggregate.h
@@ -114,6 +114,20 @@ class ARROW_EXPORT VarianceOptions : public FunctionOptions {
uint32_t min_count;
};
+/// \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);
+ 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 less than this many non-null values are observed, emit null.
+ uint32_t min_count;
+};
+
/// \brief Control Quantile kernel behavior
///
/// By default, returns the median value.
@@ -503,6 +517,34 @@ Result<Datum> Variance(const Datum& value,
const VarianceOptions& options = VarianceOptions::Defaults(),
ExecContext* ctx = NULLPTR);
+/// \brief Calculate the skewness of a numeric array
+///
+/// \param[in] value input datum, expecting Array or ChunkedArray
+/// \param[in] options see SkewOptions for more information
+/// \param[in] ctx the function execution context, optional
+/// \return datum of the computed skewness as a DoubleScalar
+///
+/// \since 20.0.0
+/// \note API not yet finalized
+ARROW_EXPORT
+Result<Datum> Skew(const Datum& value,
+ const SkewOptions& options = SkewOptions::Defaults(),
+ ExecContext* ctx = NULLPTR);
+
+/// \brief Calculate the kurtosis of a numeric array
+///
+/// \param[in] value input datum, expecting Array or ChunkedArray
+/// \param[in] options see SkewOptions for more information
+/// \param[in] ctx the function execution context, optional
+/// \return datum of the computed kurtosis as a DoubleScalar
+///
+/// \since 20.0.0
+/// \note API not yet finalized
+ARROW_EXPORT
+Result<Datum> Kurtosis(const Datum& value,
+ const SkewOptions& options = SkewOptions::Defaults(),
+ ExecContext* ctx = NULLPTR);
+
/// \brief Calculate the quantiles of a numeric array
///
/// \param[in] value input datum, expecting Array or ChunkedArray
diff --git a/cpp/src/arrow/compute/kernels/aggregate_internal.h b/cpp/src/arrow/compute/kernels/aggregate_internal.h
index 9dab049821..23aa20eddc 100644
--- a/cpp/src/arrow/compute/kernels/aggregate_internal.h
+++ b/cpp/src/arrow/compute/kernels/aggregate_internal.h
@@ -17,6 +17,9 @@
#pragma once
+#include <cmath>
+#include <initializer_list>
+
#include "arrow/compute/kernels/util_internal.h"
#include "arrow/type.h"
#include "arrow/type_traits.h"
@@ -24,9 +27,7 @@
#include "arrow/util/int128_internal.h"
#include "arrow/util/logging.h"
-namespace arrow {
-namespace compute {
-namespace internal {
+namespace arrow::compute::internal {
// Find the largest compatible primitive type for a primitive type.
template <typename I, typename Enable = void>
@@ -254,6 +255,4 @@ SumType SumArray(const ArraySpan& data) {
data, [](ValueType v) { return static_cast<SumType>(v); });
}
-} // namespace internal
-} // namespace compute
-} // namespace arrow
+} // namespace arrow::compute::internal
diff --git a/cpp/src/arrow/compute/kernels/aggregate_test.cc b/cpp/src/arrow/compute/kernels/aggregate_test.cc
index 36bb783d35..d64c740a8e 100644
--- a/cpp/src/arrow/compute/kernels/aggregate_test.cc
+++ b/cpp/src/arrow/compute/kernels/aggregate_test.cc
@@ -42,6 +42,7 @@
#include "arrow/testing/gtest_util.h"
#include "arrow/testing/matchers.h"
+#include "arrow/testing/math.h"
#include "arrow/testing/random.h"
#include "arrow/util/logging.h"
@@ -3386,6 +3387,9 @@ TEST_F(TestVarStdKernelMergeStability, Basics) {
#ifndef __MINGW32__ // MinGW has precision issues
// XXX: The reference value from numpy is actually wrong due to floating
// point limits. The correct result should equals variance(90, 0) = 4050.
+ // The problem is that the mean is not exactly representable as floating-point,
+ // and that small inaccuracy produces a large deviation when plugged into the M2
+ // calculation.
std::vector<std::string> chunks = {"[40000008000000490]", "[40000008000000400]"};
this->AssertVarStdIs(chunks, options, 3904.0);
#endif
@@ -3430,12 +3434,21 @@ TEST_F(TestVarStdKernelUInt32, Basics) {
this->AssertVarStdIs("[0, 0, 4294967295]", options, 6.148914688373205e+18);
}
-// https://en.wikipedia.org/wiki/Kahan_summation_algorithm
void KahanSum(double& sum, double& adjust, double addend) {
- double y = addend - adjust;
- double t = sum + y;
- adjust = (t - sum) - y;
- sum = t;
+ // Backported enhancement from Neumaier's algorithm: consider case where
+ // sum is small compared to addend.
+ // https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements
+ if (abs(sum) >= abs(addend)) {
+ double y = addend - adjust;
+ double t = sum + y;
+ adjust = (t - sum) - y;
+ sum = t;
+ } else {
+ double y = sum - adjust;
+ double t = addend + y;
+ adjust = (t - addend) - y;
+ sum = t;
+ }
}
// Calculate reference variance with Welford's online algorithm + Kahan summation
@@ -3534,7 +3547,8 @@ TEST_F(TestVarStdKernelIntegerLength, Basics) {
TEST(TestVarStdKernel, Decimal) {
// Effectively treated as double, sanity check results here
- for (const auto& ty : {decimal128(3, 2), decimal256(3, 2)}) {
+ for (const auto& ty :
+ {decimal32(3, 2), decimal64(3, 2), decimal128(3, 2), decimal256(3, 2)}) {
CheckVarStd(ArrayFromJSON(ty, R"(["1.00"])"), VarianceOptions(), 0);
CheckVarStd(ArrayFromJSON(ty, R"([null, "1.00", "2.00", "3.00"])"), VarianceOptions(),
0.6666666666666666);
@@ -3544,6 +3558,154 @@ TEST(TestVarStdKernel, Decimal) {
}
}
+//
+// Skew and Kurtosis
+//
+
+constexpr int kSkewUlps = 3;
+constexpr int kKurtosisUlps = 6;
+
+void CheckSkewKurtosis(const Datum& array, const SkewOptions& options,
+ double expected_skew, double expected_kurtosis, int n_ulps = -1) {
+ ARROW_SCOPED_TRACE("type = ", *array.type());
+ ASSERT_OK_AND_ASSIGN(Datum out_skew, Skew(array, options));
+ ASSERT_OK_AND_ASSIGN(Datum out_kurtosis, Kurtosis(array, options));
+ const auto& skew = checked_cast<const DoubleScalar&>(*out_skew.scalar());
+ const auto& kurtosis = checked_cast<const DoubleScalar&>(*out_kurtosis.scalar());
+ ASSERT_TRUE(skew.is_valid && kurtosis.is_valid);
+ AssertWithinUlp(expected_skew, skew.value, n_ulps >= 0 ? n_ulps : kSkewUlps);
+ AssertWithinUlp(expected_kurtosis, kurtosis.value,
+ n_ulps >= 0 ? n_ulps : kKurtosisUlps);
+}
+
+class TestSkewKurtosis : public ::testing::Test {
+ public:
+ void AssertSkewKurtosisAre(const Array& array, const SkewOptions& options,
+ double expected_skew, double expected_kurtosis,
+ int n_ulps = -1) {
+ CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+ }
+
+ void AssertSkewKurtosisAre(const std::shared_ptr<ChunkedArray>& array,
+ const SkewOptions& options, double expected_skew,
+ double expected_kurtosis, int n_ulps = -1) {
+ CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+ }
+
+ void AssertSkewKurtosisAre(const std::shared_ptr<DataType>& type, std::string_view json,
+ const SkewOptions& options, double expected_skew,
+ double expected_kurtosis, int n_ulps = -1) {
+ auto array = ArrayFromJSON(type, json);
+ CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+ }
+
+ void AssertSkewKurtosisAre(const std::shared_ptr<DataType>& type,
+ const std::vector<std::string>& json,
+ const SkewOptions& options, double expected_skew,
+ double expected_kurtosis, int n_ulps = -1) {
+ auto array = ChunkedArrayFromJSON(type, json);
+ CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+ }
+
+ void AssertSkewKurtosisInvalid(const Array& array, const SkewOptions& options) {
+ AssertSkewKurtosisInvalidInternal(array, options);
+ }
+
+ void AssertSkewKurtosisInvalid(const std::shared_ptr<ChunkedArray>& array,
+ const SkewOptions& options) {
+ AssertSkewKurtosisInvalidInternal(array, options);
+ }
+
+ void AssertSkewKurtosisInvalid(const std::shared_ptr<DataType>& type,
+ std::string_view json, const SkewOptions& options) {
+ auto array = ArrayFromJSON(type, json);
+ AssertSkewKurtosisInvalidInternal(array, options);
+ }
+
+ void AssertSkewKurtosisInvalid(const std::shared_ptr<DataType>& type,
+ const std::vector<std::string>& json,
+ const SkewOptions& options) {
+ auto array = ChunkedArrayFromJSON(type, json);
+ AssertSkewKurtosisInvalidInternal(array, options);
+ }
+
+ private:
+ void AssertSkewKurtosisInvalidInternal(const Datum& array, const SkewOptions& options) {
+ ASSERT_OK_AND_ASSIGN(Datum out_skew, Skew(array, options));
+ ASSERT_OK_AND_ASSIGN(Datum out_kurtosis, Kurtosis(array, options));
+ const auto& skew = checked_cast<const DoubleScalar&>(*out_skew.scalar());
+ const auto& kurtosis = checked_cast<const DoubleScalar&>(*out_kurtosis.scalar());
+ ASSERT_FALSE(skew.is_valid || kurtosis.is_valid);
+ }
+};
+
+TEST_F(TestSkewKurtosis, Basics) {
+ // Test sample from SciPy, with results obtained using numpy.float128
+ auto options = SkewOptions::Defaults();
+ AssertSkewKurtosisAre(float64(), "[1.165, 0.6268, 0.0751, 0.3516, -0.6965]", options,
+ -0.29322304336607355496, -0.83411431970273759);
+ // Results are slightly different because the input doesn't losslessly convert
+ // to float32.
+ AssertSkewKurtosisAre(float32(), "[1.165, 0.6268, 0.0751, 0.3516, -0.6965]", options,
+ -0.2932230870440958164, -0.8341143229437093939);
+}
+
+TEST_F(TestSkewKurtosis, Chunked) {
+ auto options = SkewOptions::Defaults();
+ AssertSkewKurtosisAre(float64(), {"[1.165, 0.6268]", "[]", "[0.0751, 0.3516, -0.6965]"},
+ options, -0.29322304336607355496, -0.83411431970273759);
+ AssertSkewKurtosisAre(float32(), {"[1.165, 0.6268]", "[]", "[0.0751, 0.3516, -0.6965]"},
+ options, -0.2932230870440958164, -0.8341143229437093939);
+}
+
+TEST_F(TestSkewKurtosis, Decimal) {
+ auto options = SkewOptions::Defaults();
+ for (auto type :
+ {decimal32(5, 4), decimal64(5, 4), decimal128(5, 4), decimal256(5, 4)}) {
+ AssertSkewKurtosisAre(type, R"(["1.1650", "0.6268", "0.0751", "0.3516", "-0.6965"])",
+ options, -0.29322304336607355496, -0.83411431970273759);
+ }
+}
+
+TEST_F(TestSkewKurtosis, Integral) {
+ auto options = SkewOptions::Defaults();
+ for (auto type : IntTypes()) {
+ AssertSkewKurtosisAre(type, "[1, 2, 3, 5]", options, 0.4346507595746657,
+ -1.1542857142857144);
+ }
+}
+
+TEST_F(TestSkewKurtosis, SpecialCases) {
+ auto options = SkewOptions::Defaults();
+ for (auto type : {float64(), float32()}) {
+ AssertSkewKurtosisAre(type, "[0, 1, 2]", options, 0.0, -1.5, /*n_ulps=*/0);
+ AssertSkewKurtosisAre(type, "[1]", options, std::nan(""), std::nan(""));
+ AssertSkewKurtosisAre(type, "[1, 1, 1, 1, 1, 1]", options, std::nan(""),
+ std::nan(""));
+ }
+}
+
+TEST_F(TestSkewKurtosis, Options) {
+ for (auto type : {float64(), float32()}) {
+ auto options = SkewOptions::Defaults();
+ AssertSkewKurtosisInvalid(type, "[]", options);
+ AssertSkewKurtosisInvalid(type, std::vector<std::string>{}, 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.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);
+ options.skip_nulls = false;
+ AssertSkewKurtosisInvalid(type, "[0, 1, null, 2]", options);
+ AssertSkewKurtosisInvalid(type, {"[0, 1]", "[]", "[null, 2]"}, options);
+ options.skip_nulls = true;
+ options.min_count = 4;
+ AssertSkewKurtosisInvalid(type, "[0, 1, null, 2]", options);
+ AssertSkewKurtosisInvalid(type, {"[0, 1]", "[]", "[null, 2]"}, options);
+ }
+}
+
//
// Quantile
//
diff --git a/cpp/src/arrow/compute/kernels/aggregate_var_std.cc b/cpp/src/arrow/compute/kernels/aggregate_var_std.cc
index e4189f9b62..8d2da195b0 100644
--- a/cpp/src/arrow/compute/kernels/aggregate_var_std.cc
+++ b/cpp/src/arrow/compute/kernels/aggregate_var_std.cc
@@ -15,18 +15,21 @@
// specific language governing permissions and limitations
// under the License.
+#include <array>
#include <cmath>
+#include <type_traits>
#include "arrow/compute/api_aggregate.h"
#include "arrow/compute/kernels/aggregate_internal.h"
#include "arrow/compute/kernels/aggregate_var_std_internal.h"
#include "arrow/compute/kernels/common_internal.h"
#include "arrow/util/bit_run_reader.h"
+#include "arrow/util/checked_cast.h"
#include "arrow/util/int128_internal.h"
-namespace arrow {
-namespace compute {
-namespace internal {
+namespace arrow::compute::internal {
+
+using ::arrow::internal::checked_cast;
namespace {
@@ -34,13 +37,14 @@ using arrow::internal::int128_t;
using arrow::internal::VisitSetBitRunsVoid;
template <typename ArrowType>
-struct VarStdState {
+struct MomentsState {
using ArrayType = typename TypeTraits<ArrowType>::ArrayType;
using CType = typename TypeTraits<ArrowType>::CType;
- using ThisType = VarStdState<ArrowType>;
+ using SumType = typename internal::GetSumType<ArrowType>::SumType;
+ using ThisType = MomentsState<ArrowType>;
- explicit VarStdState(int32_t decimal_scale, VarianceOptions options)
- : decimal_scale(decimal_scale), options(options) {}
+ MomentsState(int level, int32_t decimal_scale, bool skip_nulls)
+ : level(level), decimal_scale(decimal_scale), skip_nulls(skip_nulls) {}
template <typename T>
double ToDouble(T value) const {
@@ -51,89 +55,90 @@ struct VarStdState {
double ToDouble(const Decimal128& value) const { return value.ToDouble(decimal_scale); }
double ToDouble(const Decimal256& value) const { return value.ToDouble(decimal_scale); }
- // float/double/int64/decimal: calculate `m2` (sum((X-mean)^2)) with `two pass
- // algorithm`
- // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass_algorithm
- template <typename T = ArrowType>
- enable_if_t<is_floating_type<T>::value || (sizeof(CType) > 4) ||
- (!is_integer_type<T>::value && sizeof(CType) == 4)>
- Consume(const ArraySpan& array) {
+ int64_t count() const { return moments.count; }
+
+ void Consume(const ArraySpan& array) {
+ constexpr bool kCanUseIntArithmetic = std::is_integral_v<CType> && sizeof(CType) <= 4;
+
this->all_valid = array.GetNullCount() == 0;
- int64_t count = array.length - array.GetNullCount();
- if (count == 0 || (!this->all_valid && !options.skip_nulls)) {
+ int64_t valid_count = array.length - array.GetNullCount();
+ if (valid_count == 0 || (!this->all_valid && !this->skip_nulls)) {
return;
}
- using SumType = typename internal::GetSumType<T>::SumType;
+ if constexpr (kCanUseIntArithmetic) {
+ if (level == 2) {
+ // int32/16/8: textbook one pass algorithm for M2 with integer arithmetic
+
+ // max number of elements that sum will not overflow int64 (2Gi int32 elements)
+ // for uint32: 0 <= sum < 2^63 (int64 >= 0)
+ // for int32: -2^62 <= sum < 2^62
+ constexpr int64_t kMaxChunkLength = 1ULL << (63 - sizeof(CType) * 8);
+ int64_t start_index = 0;
+
+ ArraySpan slice = array;
+ while (valid_count > 0) {
+ // process in chunks that overflow will never happen
+ slice.SetSlice(start_index + array.offset,
+ std::min(kMaxChunkLength, array.length - start_index));
+ const int64_t count = slice.length - slice.GetNullCount();
+ start_index += slice.length;
+ valid_count -= count;
+
+ if (count > 0) {
+ IntegerVarStd var_std;
+ const CType* values = slice.GetValues<CType>(1);
+ VisitSetBitRunsVoid(slice.buffers[0].data, slice.offset, slice.length,
+ [&](int64_t pos, int64_t len) {
+ for (int64_t i = 0; i < len; ++i) {
+ const auto value = values[pos + i];
+ var_std.ConsumeOne(value);
+ }
+ });
+
+ // merge variance
+ auto slice_moments = Moments(var_std.count, var_std.mean(), var_std.m2());
+ this->moments.MergeFrom(level, slice_moments);
+ }
+ }
+ return;
+ }
+ }
+
+ // float/double/int64/decimal: calculate each moment in a separate pass.
+ // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass_algorithm
SumType sum = internal::SumArray<CType, SumType, SimdLevel::NONE>(array);
- const double mean = ToDouble(sum) / count;
+ const double mean = ToDouble(sum) / valid_count;
const double m2 = internal::SumArray<CType, double, SimdLevel::NONE>(
array, [this, mean](CType value) {
const double v = ToDouble(value);
return (v - mean) * (v - mean);
});
-
- ThisType state(decimal_scale, options);
- state.count = count;
- state.mean = mean;
- state.m2 = m2;
- this->MergeFrom(state);
- }
-
- // int32/16/8: textbook one pass algorithm with integer arithmetic
- template <typename T = ArrowType>
- enable_if_t<is_integer_type<T>::value && (sizeof(CType) <= 4)> Consume(
- const ArraySpan& array) {
- // max number of elements that sum will not overflow int64 (2Gi int32 elements)
- // for uint32: 0 <= sum < 2^63 (int64 >= 0)
- // for int32: -2^62 <= sum < 2^62
- constexpr int64_t max_length = 1ULL << (63 - sizeof(CType) * 8);
-
- this->all_valid = array.GetNullCount() == 0;
- if (!this->all_valid && !options.skip_nulls) return;
- int64_t start_index = 0;
- int64_t valid_count = array.length - array.GetNullCount();
-
- ArraySpan slice = array;
- while (valid_count > 0) {
- // process in chunks that overflow will never happen
- slice.SetSlice(start_index + array.offset,
- std::min(max_length, array.length - start_index));
- const int64_t count = slice.length - slice.GetNullCount();
- start_index += slice.length;
- valid_count -= count;
-
- if (count > 0) {
- IntegerVarStd<ArrowType> var_std;
- const CType* values = slice.GetValues<CType>(1);
- VisitSetBitRunsVoid(slice.buffers[0].data, slice.offset, slice.length,
- [&](int64_t pos, int64_t len) {
- for (int64_t i = 0; i < len; ++i) {
- const auto value = values[pos + i];
- var_std.ConsumeOne(value);
- }
- });
-
- // merge variance
- ThisType state(decimal_scale, options);
- state.count = var_std.count;
- state.mean = var_std.mean();
- state.m2 = var_std.m2();
- this->MergeFrom(state);
+ double m3 = 0, m4 = 0;
+ if (level >= 3) {
+ m3 = internal::SumArray<CType, double, SimdLevel::NONE>(
+ array, [this, mean](CType value) {
+ const double v = ToDouble(value);
+ return (v - mean) * (v - mean) * (v - mean);
+ });
+ if (level >= 4) {
+ m4 = internal::SumArray<CType, double, SimdLevel::NONE>(
+ array, [this, mean](CType value) {
+ const double v = ToDouble(value);
+ return (v - mean) * (v - mean) * (v - mean) * (v - mean);
+ });
}
}
+ this->moments.MergeFrom(level, Moments(valid_count, mean, m2, m3, m4));
}
- // Scalar: textbook algorithm
void Consume(const Scalar& scalar, const int64_t count) {
- this->m2 = 0;
if (scalar.is_valid) {
- this->count = count;
- this->mean = ToDouble(UnboxScalar<ArrowType>::Unbox(scalar));
+ double value = ToDouble(UnboxScalar<ArrowType>::Unbox(scalar));
+ this->moments = Moments::FromScalar(level, value, count);
} else {
- this->count = 0;
- this->mean = 0;
+ this->moments = Moments();
this->all_valid = false;
}
}
@@ -142,35 +147,38 @@ struct VarStdState {
// https://www.emathzone.com/tutorials/basic-statistics/combined-variance.html
void MergeFrom(const ThisType& state) {
this->all_valid = this->all_valid && state.all_valid;
- if (state.count == 0) {
- return;
- }
- if (this->count == 0) {
- this->count = state.count;
- this->mean = state.mean;
- this->m2 = state.m2;
- return;
- }
- MergeVarStd(this->count, this->mean, state.count, state.mean, state.m2, &this->count,
- &this->mean, &this->m2);
+ this->moments.MergeFrom(level, state.moments);
}
+ const int level;
const int32_t decimal_scale;
- const VarianceOptions options;
- int64_t count = 0;
- double mean = 0;
- double m2 = 0; // m2 = count*s2 = sum((X-mean)^2)
+ const bool skip_nulls;
+ Moments moments;
bool all_valid = true;
};
template <typename ArrowType>
-struct VarStdImpl : public ScalarAggregator {
- using ThisType = VarStdImpl<ArrowType>;
+struct StatisticImpl : public ScalarAggregator {
+ using ThisType = StatisticImpl<ArrowType>;
using ArrayType = typename TypeTraits<ArrowType>::ArrayType;
- explicit VarStdImpl(int32_t decimal_scale, const std::shared_ptr<DataType>& out_type,
- const VarianceOptions& options, VarOrStd return_type)
- : out_type(out_type), state(decimal_scale, options), return_type(return_type) {}
+ StatisticImpl(StatisticType stat_type, int32_t decimal_scale,
+ const std::shared_ptr<DataType>& out_type, const VarianceOptions& options)
+ : out_type(out_type),
+ stat_type(stat_type),
+ skip_nulls(options.skip_nulls),
+ min_count(options.min_count),
+ ddof(options.ddof),
+ state(moments_level_for_statistic(stat_type), decimal_scale, skip_nulls) {}
+
+ StatisticImpl(StatisticType stat_type, int32_t decimal_scale,
+ const std::shared_ptr<DataType>& out_type, const SkewOptions& options)
+ : out_type(out_type),
+ stat_type(stat_type),
+ skip_nulls(options.skip_nulls),
+ min_count(options.min_count),
+ ddof(0),
+ state(moments_level_for_statistic(stat_type), decimal_scale, skip_nulls) {}
Status Consume(KernelContext*, const ExecSpan& batch) override {
if (batch[0].is_array()) {
@@ -188,92 +196,90 @@ struct VarStdImpl : public ScalarAggregator {
}
Status Finalize(KernelContext*, Datum* out) override {
- if (state.count <= state.options.ddof || state.count < state.options.min_count ||
- (!state.all_valid && !state.options.skip_nulls)) {
+ if (state.count() <= ddof || state.count() < min_count ||
+ (!state.all_valid && !skip_nulls)) {
out->value = std::make_shared<DoubleScalar>();
} else {
- double var = state.m2 / (state.count - state.options.ddof);
- out->value =
- std::make_shared<DoubleScalar>(return_type == VarOrStd::Var ? var : sqrt(var));
+ switch (stat_type) {
+ case StatisticType::Std:
+ out->value = std::make_shared<DoubleScalar>(state.moments.Stddev(ddof));
+ break;
+ case StatisticType::Var:
+ out->value = std::make_shared<DoubleScalar>(state.moments.Variance(ddof));
+ break;
+ case StatisticType::Skew:
+ out->value = std::make_shared<DoubleScalar>(state.moments.Skew());
+ break;
+ case StatisticType::Kurtosis:
+ out->value = std::make_shared<DoubleScalar>(state.moments.Kurtosis());
+ break;
+ default:
+ return Status::NotImplemented("Unsupported statistic type ",
+ static_cast<int>(stat_type));
+ }
}
return Status::OK();
}
std::shared_ptr<DataType> out_type;
- VarStdState<ArrowType> state;
- VarOrStd return_type;
-};
-
-struct VarStdInitState {
- std::unique_ptr<KernelState> state;
- KernelContext* ctx;
- const DataType& in_type;
- const std::shared_ptr<DataType>& out_type;
- const VarianceOptions& options;
- VarOrStd return_type;
-
- VarStdInitState(KernelContext* ctx, const DataType& in_type,
- const std::shared_ptr<DataType>& out_type,
- const VarianceOptions& options, VarOrStd return_type)
- : ctx(ctx),
- in_type(in_type),
- out_type(out_type),
- options(options),
- return_type(return_type) {}
-
- Status Visit(const DataType&) {
- return Status::NotImplemented("No variance/stddev implemented");
- }
-
- Status Visit(const HalfFloatType&) {
- return Status::NotImplemented("No variance/stddev implemented");
- }
-
- template <typename Type>
- enable_if_number<Type, Status> Visit(const Type&) {
- state.reset(
- new VarStdImpl<Type>(/*decimal_scale=*/0, out_type, options, return_type));
- return Status::OK();
- }
-
- template <typename Type>
- enable_if_decimal<Type, Status> Visit(const Type&) {
- state.reset(new VarStdImpl<Type>(checked_cast<const DecimalType&>(in_type).scale(),
- out_type, options, return_type));
- return Status::OK();
- }
-
- Result<std::unique_ptr<KernelState>> Create() {
- RETURN_NOT_OK(VisitTypeInline(in_type, this));
- return std::move(state);
- }
+ StatisticType stat_type;
+ bool skip_nulls;
+ uint32_t min_count;
+ int ddof = 0;
+ MomentsState<ArrowType> state;
};
-Result<std::unique_ptr<KernelState>> StddevInit(KernelContext* ctx,
- const KernelInitArgs& args) {
- VarStdInitState visitor(
- ctx, *args.inputs[0].type, args.kernel->signature->out_type().type(),
- static_cast<const VarianceOptions&>(*args.options), VarOrStd::Std);
- return visitor.Create();
+template <typename OptionsType>
+Result<std::unique_ptr<KernelState>> StatisticInit(
+ StatisticType stat_type, const DataType& in_type,
+ const std::shared_ptr<DataType>& out_type, const OptionsType& options) {
+ auto make_kernel_state = [&](auto&& type, int32_t decimal_scale = 0) {
+ using Type = std::decay_t<decltype(type)>;
+ return std::unique_ptr<KernelState>(
+ new StatisticImpl<Type>(stat_type, decimal_scale, out_type, options));
+ };
+
+ auto visit = [&](auto&& type) -> Result<std::unique_ptr<KernelState>> {
+ using Type = std::decay_t<decltype(type)>;
+ // Decimals
+ if constexpr (is_decimal_type<Type>::value) {
+ return make_kernel_state(type, type.scale());
+ }
+ // Numbers (except half-float)
+ if constexpr (is_number_type<Type>::value && !is_half_float_type<Type>::value) {
+ return make_kernel_state(type);
+ }
+ return Status::NotImplemented("No variance/stddev implemented for ",
+ in_type.ToString());
+ };
+ return VisitType(in_type, visit);
}
-Result<std::unique_ptr<KernelState>> VarianceInit(KernelContext* ctx,
- const KernelInitArgs& args) {
- VarStdInitState visitor(
- ctx, *args.inputs[0].type, args.kernel->signature->out_type().type(),
- static_cast<const VarianceOptions&>(*args.options), VarOrStd::Var);
- return visitor.Create();
+template <typename OptionsType, StatisticType kStatType>
+Result<std::unique_ptr<KernelState>> StatisticInit(KernelContext* ctx,
+ const KernelInitArgs& args) {
+ const DataType& in_type = *args.inputs[0].type;
+ const std::shared_ptr<DataType>& out_type = args.kernel->signature->out_type().type();
+ const OptionsType& options = checked_cast<const OptionsType&>(*args.options);
+
+ return StatisticInit(kStatType, in_type, out_type, options);
}
-void AddVarStdKernels(KernelInit init,
- const std::vector<std::shared_ptr<DataType>>& types,
- ScalarAggregateFunction* func) {
+void AddStatisticKernels(KernelInit init,
+ const std::vector<std::shared_ptr<DataType>>& types,
+ ScalarAggregateFunction* func) {
for (const auto& ty : types) {
auto sig = KernelSignature::Make({InputType(ty->id())}, float64());
AddAggKernel(std::move(sig), init, func);
}
}
+void AddStatisticKernels(KernelInit init, ScalarAggregateFunction* func) {
+ AddStatisticKernels(init, NumericTypes(), func);
+ AddStatisticKernels(
+ init, {decimal32(1, 1), decimal64(1, 1), decimal128(1, 1), decimal256(1, 1)}, func);
+}
+
const FunctionDoc stddev_doc{
"Calculate the standard deviation of a numeric array",
("The number of degrees of freedom can be controlled using VarianceOptions.\n"
@@ -292,21 +298,53 @@ const FunctionDoc variance_doc{
{"array"},
"VarianceOptions"};
+const FunctionDoc skew_doc{
+ "Calculate the skewness of a numeric array",
+ ("Nulls are ignored by default. If there are not enough non-null values\n"
+ "in the array to satisfy `min_count`, null is returned.\n"
+ "The behavior of nulls and the `min_count` parameter can be changed\n"
+ "in SkewOptions."),
+ {"array"},
+ "SkewOptions"};
+
+const FunctionDoc kurtosis_doc{
+ "Calculate the kurtosis of a numeric array",
+ ("Nulls are ignored by default. If there are not enough non-null values\n"
+ "in the array to satisfy `min_count`, null is returned.\n"
+ "The behavior of nulls and the `min_count` parameter can be changed\n"
+ "in SkewOptions."),
+ {"array"},
+ "SkewOptions"};
+
std::shared_ptr<ScalarAggregateFunction> AddStddevAggKernels() {
- static auto default_std_options = VarianceOptions::Defaults();
+ static const auto default_std_options = VarianceOptions::Defaults();
auto func = std::make_shared<ScalarAggregateFunction>("stddev", Arity::Unary(),
stddev_doc, &default_std_options);
- AddVarStdKernels(StddevInit, NumericTypes(), func.get());
- AddVarStdKernels(StddevInit, {decimal128(1, 1), decimal256(1, 1)}, func.get());
+ AddStatisticKernels(StatisticInit<VarianceOptions, StatisticType::Std>, func.get());
return func;
}
std::shared_ptr<ScalarAggregateFunction> AddVarianceAggKernels() {
- static auto default_var_options = VarianceOptions::Defaults();
+ static const auto default_var_options = VarianceOptions::Defaults();
auto func = std::make_shared<ScalarAggregateFunction>(
"variance", Arity::Unary(), variance_doc, &default_var_options);
- AddVarStdKernels(VarianceInit, NumericTypes(), func.get());
- AddVarStdKernels(VarianceInit, {decimal128(1, 1), decimal256(1, 1)}, func.get());
+ AddStatisticKernels(StatisticInit<VarianceOptions, StatisticType::Var>, func.get());
+ return func;
+}
+
+std::shared_ptr<ScalarAggregateFunction> AddSkewAggKernels() {
+ static const auto default_options = SkewOptions::Defaults();
+ auto func = std::make_shared<ScalarAggregateFunction>("skew", Arity::Unary(), skew_doc,
+ &default_options);
+ AddStatisticKernels(StatisticInit<SkewOptions, StatisticType::Skew>, func.get());
+ return func;
+}
+
+std::shared_ptr<ScalarAggregateFunction> AddKurtosisAggKernels() {
+ static const auto default_options = SkewOptions::Defaults();
+ auto func = std::make_shared<ScalarAggregateFunction>("kurtosis", Arity::Unary(),
+ kurtosis_doc, &default_options);
+ AddStatisticKernels(StatisticInit<SkewOptions, StatisticType::Kurtosis>, func.get());
return func;
}
@@ -315,8 +353,8 @@ std::shared_ptr<ScalarAggregateFunction> AddVarianceAggKernels() {
void RegisterScalarAggregateVariance(FunctionRegistry* registry) {
DCHECK_OK(registry->AddFunction(AddVarianceAggKernels()));
DCHECK_OK(registry->AddFunction(AddStddevAggKernels()));
+ DCHECK_OK(registry->AddFunction(AddSkewAggKernels()));
+ DCHECK_OK(registry->AddFunction(AddKurtosisAggKernels()));
}
-} // namespace internal
-} // namespace compute
-} // namespace arrow
+} // namespace arrow::compute::internal
diff --git a/cpp/src/arrow/compute/kernels/aggregate_var_std_internal.h b/cpp/src/arrow/compute/kernels/aggregate_var_std_internal.h
index 675ebfd91d..f7c35bea96 100644
--- a/cpp/src/arrow/compute/kernels/aggregate_var_std_internal.h
+++ b/cpp/src/arrow/compute/kernels/aggregate_var_std_internal.h
@@ -18,24 +18,27 @@
#pragma once
#include "arrow/util/int128_internal.h"
+#include "arrow/util/logging.h"
+#include "arrow/util/math_internal.h"
-namespace arrow {
-namespace compute {
-namespace internal {
+#include <cmath>
+#include <type_traits>
+
+namespace arrow::compute::internal {
using arrow::internal::int128_t;
// Accumulate sum/squared sum (using naive summation)
// Shared implementation between scalar/hash aggregate variance/stddev kernels
-template <typename ArrowType>
struct IntegerVarStd {
- using c_type = typename ArrowType::c_type;
-
int64_t count = 0;
int64_t sum = 0;
int128_t square_sum = 0;
- void ConsumeOne(const c_type value) {
+ template <typename Integer>
+ void ConsumeOne(Integer value) {
+ static_assert(std::is_integral_v<Integer>);
+ static_assert(sizeof(Integer) <= 4);
sum += value;
square_sum += static_cast<uint64_t>(value) * value;
count++;
@@ -53,16 +56,87 @@ struct IntegerVarStd {
}
};
-static inline void MergeVarStd(int64_t count1, double mean1, int64_t count2, double mean2,
- double m22, int64_t* out_count, double* out_mean,
- double* out_m2) {
- double mean = (mean1 * count1 + mean2 * count2) / (count1 + count2);
- *out_m2 += m22 + count1 * (mean1 - mean) * (mean1 - mean) +
- count2 * (mean2 - mean) * (mean2 - mean);
- *out_count += count2;
- *out_mean = mean;
+enum class StatisticType { Var, Std, Skew, Kurtosis };
+
+constexpr int moments_level_for_statistic(StatisticType s) {
+ switch (s) {
+ case StatisticType::Skew:
+ return 3;
+ case StatisticType::Kurtosis:
+ return 4;
+ default:
+ return 2;
+ }
}
-} // namespace internal
-} // namespace compute
-} // namespace arrow
+struct Moments {
+ int64_t count = 0;
+ double mean = 0;
+ double m2 = 0; // m2 = sum((X-mean)^2)
+ double m3 = 0; // m3 = sum((X-mean)^3)
+ double m4 = 0; // m4 = sum((X-mean)^4)
+
+ Moments() = default;
+ Moments(int64_t count, double mean, double m2, double m3 = 0, double m4 = 0)
+ : count(count), mean(mean), m2(m2), m3(m3), m4(m4) {}
+
+ double Variance(int ddof) const { return m2 / (count - ddof); }
+
+ 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 Kurtosis() const {
+ // This may return NaN for m2 == 0 and m4 == 0, which is expected
+ return count * m4 / (m2 * m2) - 3;
+ }
+
+ void MergeFrom(int level, const Moments& other) { *this = Merge(level, *this, other); }
+
+ static Moments Merge(int level, const Moments& a, const Moments& b) {
+ using ::arrow::internal::NeumaierSum;
+
+ if (a.count == 0) {
+ return b;
+ } else if (b.count == 0) {
+ return a;
+ }
+
+ // Shorter aliases for readability
+ const int64_t na = a.count, nb = b.count;
+ const int64_t n = na + nb;
+ const double mean = (a.mean * na + b.mean * nb) / n;
+ // NOTE: there is a more common formula:
+ // double delta = b.mean - a.mean;
+ // double m2 = a.m2 + b.m2 + delta * delta * na * nb / n;
+ // but it gives worse results in TestVarStdKernelMergeStability.
+ const double m2 = NeumaierSum({a.m2, b.m2, na * (a.mean - mean) * (a.mean - mean),
+ nb * (b.mean - mean) * (b.mean - mean)});
+ double m3 = 0;
+ double m4 = 0;
+ // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
+ if (level >= 3) {
+ double delta = b.mean - a.mean;
+ double delta2 = delta * delta;
+ m3 = NeumaierSum({a.m3, b.m3, delta2 * delta * na * nb * (na - nb) / (n * n),
+ 3 * delta * (na * b.m2 - nb * a.m2) / n});
+ if (level >= 4) {
+ m4 = NeumaierSum(
+ {a.m4, b.m4,
+ (delta2 * delta2) * na * nb * (na * na - na * nb + nb * nb) / (n * n * n),
+ 6 * delta2 * (na * na * b.m2 + nb * nb * a.m2) / (n * n),
+ 4 * delta * (na * b.m3 - nb * a.m3) / n});
+ }
+ }
+ return Moments(n, mean, m2, m3, m4);
+ }
+
+ static Moments FromScalar(int level, double value, int64_t count) {
+ return Moments(count, /*mean=*/value, /*m2=*/0, /*m3=*/0, /*m4=*/0);
+ }
+};
+
+} // namespace arrow::compute::internal
diff --git a/cpp/src/arrow/compute/kernels/hash_aggregate.cc b/cpp/src/arrow/compute/kernels/hash_aggregate.cc
index b52f1c8286..e84c3c2dc1 100644
--- a/cpp/src/arrow/compute/kernels/hash_aggregate.cc
+++ b/cpp/src/arrow/compute/kernels/hash_aggregate.cc
@@ -134,10 +134,12 @@ HashAggregateKernel MakeUnaryKernel(KernelInit init) {
std::move(init));
}
-Status AddHashAggKernels(
- const std::vector<std::shared_ptr<DataType>>& types,
- Result<HashAggregateKernel> make_kernel(const std::shared_ptr<DataType>&),
- HashAggregateFunction* function) {
+using HashAggregateKernelFactory =
+ std::function<Result<HashAggregateKernel>(const std::shared_ptr<DataType>&)>;
+
+Status AddHashAggKernels(const std::vector<std::shared_ptr<DataType>>& types,
+ HashAggregateKernelFactory make_kernel,
+ HashAggregateFunction* function) {
for (const auto& ty : types) {
ARROW_ASSIGN_OR_RAISE(auto kernel, make_kernel(ty));
RETURN_NOT_OK(function->AddKernel(std::move(kernel)));
@@ -845,28 +847,58 @@ using GroupedMeanFactory =
using arrow::internal::int128_t;
template <typename Type>
-struct GroupedVarStdImpl : public GroupedAggregator {
+struct GroupedStatisticImpl : public GroupedAggregator {
using CType = typename TypeTraits<Type>::CType;
+ using SumType = typename internal::GetSumType<Type>::SumType;
+ // This method is defined solely to make GroupedStatisticImpl instantiable
+ // in ConsumeImpl below. It will be redefined in subclasses.
Status Init(ExecContext* ctx, const KernelInitArgs& args) override {
- options_ = *checked_cast<const VarianceOptions*>(args.options);
- if (is_decimal_type<Type>::value) {
- const int32_t scale =
+ return Status::NotImplemented("");
+ }
+
+ // Init helper for hash_variance and hash_stddev
+ 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);
+ }
+
+ // 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);
+ }
+
+ Status InitInternal(ExecContext* ctx, const KernelInitArgs& args,
+ StatisticType stat_type, int ddof, bool skip_nulls,
+ 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, scale, args.options);
+ return InitInternal(ctx, stat_type, decimal_scale, ddof, skip_nulls, min_count);
+ } else {
+ return InitInternal(ctx, stat_type, /*decimal_scale=*/0, ddof, skip_nulls,
+ min_count);
}
- return InitInternal(ctx, 0, args.options);
}
- Status InitInternal(ExecContext* ctx, int32_t decimal_scale,
- const FunctionOptions* options) {
- options_ = *checked_cast<const VarianceOptions*>(options);
+ Status InitInternal(ExecContext* ctx, StatisticType stat_type, int32_t decimal_scale,
+ int ddof, bool skip_nulls, uint32_t min_count) {
+ stat_type_ = stat_type;
+ moments_level_ = moments_level_for_statistic(stat_type_);
decimal_scale_ = decimal_scale;
+ skip_nulls_ = skip_nulls;
+ min_count_ = min_count;
+ ddof_ = ddof;
ctx_ = ctx;
pool_ = ctx->memory_pool();
counts_ = TypedBufferBuilder<int64_t>(pool_);
means_ = TypedBufferBuilder<double>(pool_);
m2s_ = TypedBufferBuilder<double>(pool_);
+ m3s_ = TypedBufferBuilder<double>(pool_);
+ m4s_ = TypedBufferBuilder<double>(pool_);
no_nulls_ = TypedBufferBuilder<bool>(pool_);
return Status::OK();
}
@@ -877,6 +909,12 @@ struct GroupedVarStdImpl : public GroupedAggregator {
RETURN_NOT_OK(counts_.Append(added_groups, 0));
RETURN_NOT_OK(means_.Append(added_groups, 0));
RETURN_NOT_OK(m2s_.Append(added_groups, 0));
+ if (moments_level_ >= 3) {
+ RETURN_NOT_OK(m3s_.Append(added_groups, 0));
+ if (moments_level_ >= 4) {
+ RETURN_NOT_OK(m4s_.Append(added_groups, 0));
+ }
+ }
RETURN_NOT_OK(no_nulls_.Append(added_groups, true));
return Status::OK();
}
@@ -894,27 +932,30 @@ struct GroupedVarStdImpl : public GroupedAggregator {
return value.ToDouble(decimal_scale_);
}
- Status Consume(const ExecSpan& batch) override { return ConsumeImpl(batch); }
+ Status Consume(const ExecSpan& batch) override {
+ constexpr bool kCanUseIntArithmetic = std::is_integral_v<CType> && sizeof(CType) <= 4;
+
+ if constexpr (kCanUseIntArithmetic) {
+ if (moments_level_ == 2) {
+ return ConsumeIntegral(batch);
+ }
+ }
+ return ConsumeGeneric(batch);
+ }
// float/double/int64/decimal: calculate `m2` (sum((X-mean)^2)) with
- // `two pass algorithm` (see aggregate_var_std.cc)
- template <typename T = Type>
- enable_if_t<is_floating_type<T>::value || (sizeof(CType) > 4) ||
- std::is_same_v<CType, Decimal32>,
- Status>
- ConsumeImpl(const ExecSpan& batch) {
- using SumType = typename internal::GetSumType<T>::SumType;
-
- GroupedVarStdImpl<Type> state;
- RETURN_NOT_OK(state.InitInternal(ctx_, decimal_scale_, &options_));
+ // two pass algorithm (see aggregate_var_std.cc)
+ Status ConsumeGeneric(const ExecSpan& batch) {
+ GroupedStatisticImpl<Type> state;
+ RETURN_NOT_OK(state.InitInternal(ctx_, stat_type_, decimal_scale_, ddof_, skip_nulls_,
+ min_count_));
RETURN_NOT_OK(state.Resize(num_groups_));
int64_t* counts = state.counts_.mutable_data();
double* means = state.means_.mutable_data();
- double* m2s = state.m2s_.mutable_data();
uint8_t* no_nulls = state.no_nulls_.mutable_data();
- // XXX this uses naive summation; we should switch to pairwise summation as was
- // done for the scalar aggregate kernel in ARROW-11567
+ // XXX this uses naive summation; we should switch to pairwise summation
+ // (as the scalar aggregate kernel does) or Kahan summation.
std::vector<SumType> sums(num_groups_);
VisitGroupedValues<Type>(
batch,
@@ -928,27 +969,34 @@ struct GroupedVarStdImpl : public GroupedAggregator {
means[i] = ToDouble(sums[i]) / counts[i];
}
+ double* m2s = state.m2s_mutable_data();
+ double* m3s = state.m3s_mutable_data();
+ double* m4s = state.m4s_mutable_data();
+ // Having distinct VisitGroupedValuesNonNull calls based on moments_level_
+ // would increase code generation for relatively little benefit.
VisitGroupedValuesNonNull<Type>(
batch, [&](uint32_t g, typename TypeTraits<Type>::CType value) {
- const double v = ToDouble(value);
- m2s[g] += (v - means[g]) * (v - means[g]);
+ const double d = ToDouble(value) - means[g];
+ const double d2 = d * d;
+ switch (moments_level_) {
+ case 4:
+ m4s[g] += d2 * d2;
+ [[fallthrough]];
+ case 3:
+ m3s[g] += d2 * d;
+ [[fallthrough]];
+ default:
+ m2s[g] += d2;
+ break;
+ }
});
- ARROW_ASSIGN_OR_RAISE(auto mapping,
- AllocateBuffer(num_groups_ * sizeof(uint32_t), pool_));
- for (uint32_t i = 0; static_cast<int64_t>(i) < num_groups_; i++) {
- mapping->template mutable_data_as<uint32_t>()[i] = i;
- }
- ArrayData group_id_mapping(uint32(), num_groups_, {nullptr, std::move(mapping)},
- /*null_count=*/0);
- return this->Merge(std::move(state), group_id_mapping);
+ return MergeSameGroups(std::move(state));
}
- // int32/16/8: textbook one pass algorithm with integer arithmetic (see
- // aggregate_var_std.cc)
- template <typename T = Type>
- enable_if_t<is_integer_type<T>::value && (sizeof(CType) <= 4), Status> ConsumeImpl(
- const ExecSpan& batch) {
+ // int32/16/8: textbook one pass algorithm to compute `m2` with integer arithmetic
+ // (see aggregate_var_std.cc)
+ Status ConsumeIntegral(const ExecSpan& batch) {
// max number of elements that sum will not overflow int64 (2Gi int32 elements)
// for uint32: 0 <= sum < 2^63 (int64 >= 0)
// for int32: -2^62 <= sum < 2^62
@@ -963,15 +1011,7 @@ struct GroupedVarStdImpl : public GroupedAggregator {
return Status::OK();
}
- std::vector<IntegerVarStd<Type>> var_std(num_groups_);
-
- ARROW_ASSIGN_OR_RAISE(auto mapping,
- AllocateBuffer(num_groups_ * sizeof(uint32_t), pool_));
- for (uint32_t i = 0; static_cast<int64_t>(i) < num_groups_; i++) {
- mapping->template mutable_data_as<uint32_t>()[i] = i;
- }
- ArrayData group_id_mapping(uint32(), num_groups_, {nullptr, std::move(mapping)},
- /*null_count=*/0);
+ std::vector<IntegerVarStd> var_std(num_groups_);
for (int64_t start_index = 0; start_index < batch.length; start_index += max_length) {
// process in chunks that overflow will never happen
@@ -979,12 +1019,13 @@ struct GroupedVarStdImpl : public GroupedAggregator {
// reset state
var_std.clear();
var_std.resize(num_groups_);
- GroupedVarStdImpl<Type> state;
- RETURN_NOT_OK(state.InitInternal(ctx_, decimal_scale_, &options_));
+ GroupedStatisticImpl<Type> state;
+ RETURN_NOT_OK(state.InitInternal(ctx_, stat_type_, decimal_scale_, ddof_,
+ skip_nulls_, min_count_));
RETURN_NOT_OK(state.Resize(num_groups_));
int64_t* other_counts = state.counts_.mutable_data();
double* other_means = state.means_.mutable_data();
- double* other_m2s = state.m2s_.mutable_data();
+ double* other_m2s = state.m2s_mutable_data();
uint8_t* other_no_nulls = state.no_nulls_.mutable_data();
if (batch[0].is_array()) {
@@ -1033,34 +1074,63 @@ struct GroupedVarStdImpl : public GroupedAggregator {
other_means[i] = var_std[i].mean();
other_m2s[i] = var_std[i].m2();
}
- RETURN_NOT_OK(this->Merge(std::move(state), group_id_mapping));
+ RETURN_NOT_OK(MergeSameGroups(std::move(state)));
}
return Status::OK();
}
Status Merge(GroupedAggregator&& raw_other,
const ArrayData& group_id_mapping) override {
- // Combine m2 from two chunks (see aggregate_var_std.cc)
- auto other = checked_cast<GroupedVarStdImpl*>(&raw_other);
+ DCHECK_EQ(group_id_mapping.length,
+ checked_cast<GroupedStatisticImpl*>(&raw_other)->num_groups_);
+ const uint32_t* g = group_id_mapping.GetValues<uint32_t>(1);
+ return MergeInternal(std::move(raw_other),
+ [g](int64_t other_g) { return g[other_g]; });
+ }
+
+ Status MergeSameGroups(GroupedAggregator&& raw_other) {
+ return MergeInternal(std::move(raw_other), [](int64_t other_g) { return other_g; });
+ }
+
+ template <typename GroupIdMapper>
+ Status MergeInternal(GroupedAggregator&& raw_other, GroupIdMapper&& group_id_mapper) {
+ // Combine moments from two chunks
+ auto other = checked_cast<GroupedStatisticImpl*>(&raw_other);
+ DCHECK_EQ(moments_level_, other->moments_level_);
int64_t* counts = counts_.mutable_data();
double* means = means_.mutable_data();
- double* m2s = m2s_.mutable_data();
+ double* m2s = m2s_mutable_data();
+ // Moments above the current level will just be ignored.
+ double* m3s = m3s_mutable_data();
+ double* m4s = m4s_mutable_data();
uint8_t* no_nulls = no_nulls_.mutable_data();
const int64_t* other_counts = other->counts_.data();
const double* other_means = other->means_.data();
- const double* other_m2s = other->m2s_.data();
+ const double* other_m2s = other->m2s_data();
+ const double* other_m3s = other->m3s_data();
+ const double* other_m4s = other->m4s_data();
const uint8_t* other_no_nulls = other->no_nulls_.data();
- auto g = group_id_mapping.GetValues<uint32_t>(1);
- for (int64_t other_g = 0; other_g < group_id_mapping.length; ++other_g, ++g) {
+ const int64_t num_other_groups = other->num_groups_;
+
+ for (int64_t other_g = 0; other_g < num_other_groups; ++other_g) {
+ const auto g = group_id_mapper(other_g);
if (!bit_util::GetBit(other_no_nulls, other_g)) {
- bit_util::ClearBit(no_nulls, *g);
+ bit_util::ClearBit(no_nulls, g);
}
if (other_counts[other_g] == 0) continue;
- MergeVarStd(counts[*g], means[*g], other_counts[other_g], other_means[other_g],
- other_m2s[other_g], &counts[*g], &means[*g], &m2s[*g]);
+ auto moments = Moments::Merge(
+ moments_level_, Moments(counts[g], means[g], m2s[g], m3s[g], m4s[g]),
+ Moments(other_counts[other_g], other_means[other_g], other_m2s[other_g],
+ other_m3s[other_g], other_m4s[other_g]));
+ counts[g] = moments.count;
+ means[g] = moments.mean;
+ // Fill moments in reverse order, in case m3s or m4s is the same as m2s.
+ m4s[g] = moments.m4;
+ m3s[g] = moments.m3;
+ m2s[g] = moments.m2;
}
return Status::OK();
}
@@ -1073,11 +1143,30 @@ struct GroupedVarStdImpl : public GroupedAggregator {
auto* results = values->mutable_data_as<double>();
const int64_t* counts = counts_.data();
- const double* m2s = m2s_.data();
+ const double* means = means_.data();
+ const double* m2s = m2s_data();
+ const double* m3s = m3s_data();
+ const double* m4s = m4s_data();
for (int64_t i = 0; i < num_groups_; ++i) {
- if (counts[i] > options_.ddof && counts[i] >= options_.min_count) {
- const double variance = m2s[i] / (counts[i] - options_.ddof);
- results[i] = result_type_ == VarOrStd::Var ? variance : std::sqrt(variance);
+ if (counts[i] > ddof_ && counts[i] >= min_count_) {
+ const auto moments = Moments(counts[i], means[i], m2s[i], m3s[i], m4s[i]);
+ switch (stat_type_) {
+ case StatisticType::Var:
+ results[i] = moments.Variance(ddof_);
+ break;
+ case StatisticType::Std:
+ results[i] = moments.Stddev(ddof_);
+ break;
+ case StatisticType::Skew:
+ results[i] = moments.Skew();
+ break;
+ case StatisticType::Kurtosis:
+ results[i] = moments.Kurtosis();
+ break;
+ default:
+ return Status::NotImplemented("Statistic type ",
+ static_cast<int>(stat_type_));
+ }
continue;
}
@@ -1090,7 +1179,7 @@ struct GroupedVarStdImpl : public GroupedAggregator {
null_count += 1;
bit_util::SetBitTo(null_bitmap->mutable_data(), i, false);
}
- if (!options_.skip_nulls) {
+ if (!skip_nulls_) {
if (null_bitmap) {
arrow::internal::BitmapAnd(null_bitmap->data(), 0, no_nulls_.data(), 0,
num_groups_, 0, null_bitmap->mutable_data());
@@ -1106,56 +1195,91 @@ struct GroupedVarStdImpl : public GroupedAggregator {
std::shared_ptr<DataType> out_type() const override { return float64(); }
- VarOrStd result_type_;
+ const double* m2s_data() const { return m2s_.data(); }
+ // If moments_level_ < 3, the values read from m3s_data() will be ignored,
+ // but we still need to point to a valid buffer of the appropriate size.
+ // The trick is to reuse m2s_, which simplifies the code.
+ const double* m3s_data() const {
+ return (moments_level_ >= 3) ? m3s_.data() : m2s_.data();
+ }
+ const double* m4s_data() const {
+ return (moments_level_ >= 4) ? m4s_.data() : m2s_.data();
+ }
+
+ double* m2s_mutable_data() { return m2s_.mutable_data(); }
+ double* m3s_mutable_data() {
+ return (moments_level_ >= 3) ? m3s_.mutable_data() : m2s_.mutable_data();
+ }
+ double* m4s_mutable_data() {
+ return (moments_level_ >= 4) ? m4s_.mutable_data() : m2s_.mutable_data();
+ }
+
+ StatisticType stat_type_;
+ int moments_level_;
int32_t decimal_scale_;
- VarianceOptions options_;
+ bool skip_nulls_;
+ uint32_t min_count_;
+ int ddof_;
int64_t num_groups_ = 0;
// m2 = count * s2 = sum((X-mean)^2)
TypedBufferBuilder<int64_t> counts_;
- TypedBufferBuilder<double> means_, m2s_;
+ TypedBufferBuilder<double> means_, m2s_, m3s_, m4s_;
TypedBufferBuilder<bool> no_nulls_;
ExecContext* ctx_;
MemoryPool* pool_;
};
-template <typename T, VarOrStd result_type>
-Result<std::unique_ptr<KernelState>> VarStdInit(KernelContext* ctx,
- const KernelInitArgs& args) {
- auto impl = std::make_unique<GroupedVarStdImpl<T>>();
- impl->result_type_ = result_type;
- RETURN_NOT_OK(impl->Init(ctx->exec_context(), args));
- // R build with openSUSE155 requires an explicit unique_ptr construction
- return std::unique_ptr<KernelState>(std::move(impl));
-}
-
-template <VarOrStd result_type>
-struct GroupedVarStdFactory {
- template <typename T, typename Enable = enable_if_t<is_integer_type<T>::value ||
- is_floating_type<T>::value ||
- is_decimal_type<T>::value>>
- Status Visit(const T&) {
- kernel = MakeKernel(std::move(argument_type), VarStdInit<T, result_type>);
- return Status::OK();
- }
+template <typename Type, typename OptionsType, StatisticType kStatType>
+struct ConcreteGroupedStatisticImpl : public GroupedStatisticImpl<Type> {
+ using GroupedStatisticImpl<Type>::InitInternal;
- Status Visit(const HalfFloatType& type) {
- return Status::NotImplemented("Computing variance/stddev of data of type ", type);
+ Status Init(ExecContext* ctx, const KernelInitArgs& args) override {
+ const auto& options = checked_cast<const OptionsType&>(*args.options);
+ return InitInternal(ctx, args, kStatType, options);
}
+};
- Status Visit(const DataType& type) {
- return Status::NotImplemented("Computing variance/stddev of data of type ", type);
- }
+template <typename Type>
+using GroupedVarianceImpl =
+ ConcreteGroupedStatisticImpl<Type, VarianceOptions, StatisticType::Var>;
+template <typename Type>
+using GroupedStddevImpl =
+ ConcreteGroupedStatisticImpl<Type, VarianceOptions, StatisticType::Std>;
+template <typename Type>
+using GroupedSkewImpl =
+ ConcreteGroupedStatisticImpl<Type, SkewOptions, StatisticType::Skew>;
+template <typename Type>
+using GroupedKurtosisImpl =
+ ConcreteGroupedStatisticImpl<Type, SkewOptions, StatisticType::Kurtosis>;
+
+template <template <typename Type> typename GroupedImpl>
+Result<HashAggregateKernel> MakeGroupedStatisticKernel(
+ const std::shared_ptr<DataType>& type) {
+ auto make_kernel = [&](auto&& type) -> Result<HashAggregateKernel> {
+ using T = std::decay_t<decltype(type)>;
+ // Supporting all number types except float16
+ if constexpr (is_integer_type<T>::value ||
+ (is_floating_type<T>::value && !is_half_float_type<T>::value) ||
+ is_decimal_type<T>::value) {
+ return MakeKernel(InputType(T::type_id), HashAggregateInit<GroupedImpl<T>>);
+ }
+ return Status::NotImplemented("Computing higher-order statistic of data of type ",
+ type);
+ };
- static Result<HashAggregateKernel> Make(const std::shared_ptr<DataType>& type) {
- GroupedVarStdFactory factory;
- factory.argument_type = type->id();
- RETURN_NOT_OK(VisitTypeInline(*type, &factory));
- return std::move(factory.kernel);
- }
+ return VisitType(*type, make_kernel);
+}
- HashAggregateKernel kernel;
- InputType argument_type;
-};
+Status AddHashAggregateStatisticKernels(HashAggregateFunction* func,
+ HashAggregateKernelFactory make_kernel) {
+ RETURN_NOT_OK(AddHashAggKernels(SignedIntTypes(), make_kernel, func));
+ RETURN_NOT_OK(AddHashAggKernels(UnsignedIntTypes(), make_kernel, func));
+ RETURN_NOT_OK(AddHashAggKernels(FloatingPointTypes(), make_kernel, func));
+ RETURN_NOT_OK(AddHashAggKernels(
+ {decimal32(1, 1), decimal64(1, 1), decimal128(1, 1), decimal256(1, 1)}, make_kernel,
+ func));
+ return Status::OK();
+}
// ----------------------------------------------------------------------
// TDigest implementation
@@ -3749,8 +3873,8 @@ const FunctionDoc hash_product_doc{
const FunctionDoc hash_mean_doc{
"Compute the mean of values in each group",
("Null values are ignored.\n"
- "For integers and floats, NaN is returned if min_count = 0 and\n"
- "there are no values. For decimals, null is returned instead."),
+ "For integers and floats, NaN is emitted if min_count = 0 and\n"
+ "there are no values in a group. For decimals, null is emitted instead."),
{"array", "group_id_array"},
"ScalarAggregateOptions"};
@@ -3758,22 +3882,38 @@ const FunctionDoc hash_stddev_doc{
"Compute the standard deviation of values in each group",
("The number of degrees of freedom can be controlled using VarianceOptions.\n"
"By default (`ddof` = 0), the population standard deviation is calculated.\n"
- "Nulls are ignored. If there are not enough non-null values in the array\n"
- "to satisfy `ddof`, null is returned."),
+ "Nulls are ignored. If there are not enough non-null values in a group\n"
+ "to satisfy `ddof`, null is emitted."),
{"array", "group_id_array"}};
const FunctionDoc hash_variance_doc{
"Compute the variance of values in each group",
("The number of degrees of freedom can be controlled using VarianceOptions.\n"
"By default (`ddof` = 0), the population variance is calculated.\n"
- "Nulls are ignored. If there are not enough non-null values in the array\n"
- "to satisfy `ddof`, null is returned."),
+ "Nulls are ignored. If there are not enough non-null values in a group\n"
+ "to satisfy `ddof`, null is emitted."),
+ {"array", "group_id_array"}};
+
+const FunctionDoc hash_skew_doc{
+ "Compute the skewness of values in each group",
+ ("Nulls are ignored by default. If there are not enough non-null values\n"
+ "in a group to satisfy `min_count`, null is emitted.\n"
+ "The behavior of nulls and the `min_count` parameter can be changed\n"
+ "in SkewOptions."),
+ {"array", "group_id_array"}};
+
+const FunctionDoc hash_kurtosis_doc{
+ "Compute the kurtosis of values in each group",
+ ("Nulls are ignored by default. If there are not enough non-null values\n"
+ "in a group to satisfy `min_count`, null is emitted.\n"
+ "The behavior of nulls and the `min_count` parameter can be changed\n"
+ "in SkewOptions."),
{"array", "group_id_array"}};
const FunctionDoc hash_tdigest_doc{
"Compute approximate quantiles of values in each group",
("The T-Digest algorithm is used for a fast approximation.\n"
- "By default, the 0.5 quantile (i.e. median) is returned.\n"
+ "By default, the 0.5 quantile (i.e. median) is emitted.\n"
"Nulls and NaNs are ignored.\n"
"Nulls are returned if there are no valid data points."),
{"array", "group_id_array"},
@@ -3873,10 +4013,11 @@ const FunctionDoc hash_pivot_doc{
} // namespace
void RegisterHashAggregateBasic(FunctionRegistry* registry) {
- static auto default_count_options = CountOptions::Defaults();
- static auto default_scalar_aggregate_options = ScalarAggregateOptions::Defaults();
- static auto default_tdigest_options = TDigestOptions::Defaults();
- static auto default_variance_options = VarianceOptions::Defaults();
+ static const auto default_count_options = CountOptions::Defaults();
+ static const auto default_scalar_aggregate_options = ScalarAggregateOptions::Defaults();
+ static const auto default_tdigest_options = TDigestOptions::Defaults();
+ static const auto default_variance_options = VarianceOptions::Defaults();
+ static const auto default_skew_options = SkewOptions::Defaults();
{
auto func = std::make_shared<HashAggregateFunction>(
@@ -3948,28 +4089,32 @@ void RegisterHashAggregateBasic(FunctionRegistry* registry) {
{
auto func = std::make_shared<HashAggregateFunction>(
"hash_stddev", Arity::Binary(), hash_stddev_doc, &default_variance_options);
- DCHECK_OK(AddHashAggKernels(SignedIntTypes(),
- GroupedVarStdFactory<VarOrStd::Std>::Make, func.get()));
- DCHECK_OK(AddHashAggKernels(UnsignedIntTypes(),
- GroupedVarStdFactory<VarOrStd::Std>::Make, func.get()));
- DCHECK_OK(AddHashAggKernels(FloatingPointTypes(),
- GroupedVarStdFactory<VarOrStd::Std>::Make, func.get()));
- DCHECK_OK(AddHashAggKernels({decimal128(1, 1), decimal256(1, 1)},
- GroupedVarStdFactory<VarOrStd::Std>::Make, func.get()));
+ DCHECK_OK(AddHashAggregateStatisticKernels(
+ func.get(), MakeGroupedStatisticKernel<GroupedStddevImpl>));
DCHECK_OK(registry->AddFunction(std::move(func)));
}
{
auto func = std::make_shared<HashAggregateFunction>(
"hash_variance", Arity::Binary(), hash_variance_doc, &default_variance_options);
- DCHECK_OK(AddHashAggKernels(SignedIntTypes(),
- GroupedVarStdFactory<VarOrStd::Var>::Make, func.get()));
- DCHECK_OK(AddHashAggKernels(UnsignedIntTypes(),
- GroupedVarStdFactory<VarOrStd::Var>::Make, func.get()));
- DCHECK_OK(AddHashAggKernels(FloatingPointTypes(),
- GroupedVarStdFactory<VarOrStd::Var>::Make, func.get()));
- DCHECK_OK(AddHashAggKernels({decimal128(1, 1), decimal256(1, 1)},
- GroupedVarStdFactory<VarOrStd::Var>::Make, func.get()));
+ DCHECK_OK(AddHashAggregateStatisticKernels(
+ func.get(), MakeGroupedStatisticKernel<GroupedVarianceImpl>));
+ DCHECK_OK(registry->AddFunction(std::move(func)));
+ }
+
+ {
+ auto func = std::make_shared<HashAggregateFunction>(
+ "hash_skew", Arity::Binary(), hash_skew_doc, &default_skew_options);
+ DCHECK_OK(AddHashAggregateStatisticKernels(
+ func.get(), MakeGroupedStatisticKernel<GroupedSkewImpl>));
+ DCHECK_OK(registry->AddFunction(std::move(func)));
+ }
+
+ {
+ auto func = std::make_shared<HashAggregateFunction>(
+ "hash_kurtosis", Arity::Binary(), hash_kurtosis_doc, &default_skew_options);
+ DCHECK_OK(AddHashAggregateStatisticKernels(
+ func.get(), MakeGroupedStatisticKernel<GroupedKurtosisImpl>));
DCHECK_OK(registry->AddFunction(std::move(func)));
}
diff --git a/cpp/src/arrow/testing/gtest_util_test.cc b/cpp/src/arrow/testing/gtest_util_test.cc
index daf071c2b3..94c8156b8f 100644
--- a/cpp/src/arrow/testing/gtest_util_test.cc
+++ b/cpp/src/arrow/testing/gtest_util_test.cc
@@ -227,12 +227,14 @@ void CheckNotWithinUlp(Float x, Float y, int n_ulp) {
TEST(TestWithinUlp, Double) {
for (double f : {0.0, 1e-20, 1.0, 2345678.9}) {
+ CheckWithinUlp(f, f, 0);
CheckWithinUlp(f, f, 1);
CheckWithinUlp(f, f, 42);
}
CheckWithinUlp(-0.0, 0.0, 1);
CheckWithinUlp(1.0, 1.0000000000000002, 1);
CheckWithinUlp(1.0, 1.0000000000000007, 3);
+ CheckNotWithinUlp(1.0, 1.0000000000000002, 0);
CheckNotWithinUlp(1.0, 1.0000000000000007, 2);
CheckNotWithinUlp(1.0, 1.0000000000000007, 1);
// left and right have a different exponent but are still very close
@@ -243,6 +245,9 @@ TEST(TestWithinUlp, Double) {
CheckWithinUlp(123.4567, 123.45670000000015, 11);
CheckNotWithinUlp(123.4567, 123.45670000000015, 10);
+ CheckWithinUlp(HUGE_VAL, HUGE_VAL, 10);
+ CheckWithinUlp(-HUGE_VAL, -HUGE_VAL, 10);
+ CheckWithinUlp(std::nan(""), std::nan(""), 10);
CheckNotWithinUlp(HUGE_VAL, -HUGE_VAL, 10);
CheckNotWithinUlp(12.34, -HUGE_VAL, 10);
CheckNotWithinUlp(12.34, std::nan(""), 10);
@@ -252,12 +257,14 @@ TEST(TestWithinUlp, Double) {
TEST(TestWithinUlp, Float) {
for (float f : {0.0f, 1e-8f, 1.0f, 123.456f}) {
+ CheckWithinUlp(f, f, 0);
CheckWithinUlp(f, f, 1);
CheckWithinUlp(f, f, 42);
}
CheckWithinUlp(-0.0f, 0.0f, 1);
CheckWithinUlp(1.0f, 1.0000001f, 1);
CheckWithinUlp(1.0f, 1.0000013f, 11);
+ CheckNotWithinUlp(1.0f, 1.0000001f, 0);
CheckNotWithinUlp(1.0f, 1.0000013f, 10);
// left and right have a different exponent but are still very close
CheckWithinUlp(1.0f, 0.99999994f, 1);
@@ -267,6 +274,9 @@ TEST(TestWithinUlp, Float) {
CheckWithinUlp(123.456f, 123.456085f, 11);
CheckNotWithinUlp(123.456f, 123.456085f, 10);
+ CheckWithinUlp(HUGE_VALF, HUGE_VALF, 10);
+ CheckWithinUlp(-HUGE_VALF, -HUGE_VALF, 10);
+ CheckWithinUlp(std::nanf(""), std::nanf(""), 10);
CheckNotWithinUlp(HUGE_VALF, -HUGE_VALF, 10);
CheckNotWithinUlp(12.34f, -HUGE_VALF, 10);
CheckNotWithinUlp(12.34f, std::nanf(""), 10);
diff --git a/cpp/src/arrow/testing/math.cc b/cpp/src/arrow/testing/math.cc
index 2cb2fcb2a9..cc22e5a4de 100644
--- a/cpp/src/arrow/testing/math.cc
+++ b/cpp/src/arrow/testing/math.cc
@@ -53,9 +53,15 @@ bool WithinUlpOneWay(Float left, Float right, int n_ulps) {
template <typename Float>
bool WithinUlpGeneric(Float left, Float right, int n_ulps) {
+ if (std::isnan(left) || std::isnan(right)) {
+ return std::isnan(left) == std::isnan(right);
+ }
if (!std::isfinite(left) || !std::isfinite(right)) {
return left == right;
}
+ if (n_ulps == 0) {
+ return left == right;
+ }
return (std::abs(left) <= std::abs(right)) ? WithinUlpOneWay(left, right, n_ulps)
: WithinUlpOneWay(right, left, n_ulps);
}
diff --git a/cpp/src/arrow/util/math_internal.h b/cpp/src/arrow/util/math_internal.h
index db8856e708..3ff30cabf2 100644
--- a/cpp/src/arrow/util/math_internal.h
+++ b/cpp/src/arrow/util/math_internal.h
@@ -17,6 +17,10 @@
#pragma once
+#include <cassert>
+#include <cmath>
+#include <initializer_list>
+
#include "arrow/util/macros.h"
#include "arrow/util/visibility.h"
@@ -33,4 +37,26 @@ namespace arrow::internal {
ARROW_EXPORT
double NormalPPF(double p);
+/// "Improved Kahan–Babuška algorithm" by Neumaier
+///
+/// https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements
+template <typename Range = std::initializer_list<double>>
+double NeumaierSum(Range&& inputs) {
+ double sum = 0, c = 0;
+ for (const double v : inputs) {
+ double t = sum + v;
+ if (std::isfinite(t)) {
+ if (std::abs(sum) >= std::abs(v)) {
+ // If sum is bigger, low-order digits of v are lost...
+ c += (sum - t) + v;
+ } else {
+ // ... else low-order digits of sum are lost.
+ c += (v - t) + sum;
+ }
+ }
+ sum = t;
+ }
+ return sum + c;
+}
+
} // namespace arrow::internal
diff --git a/cpp/src/arrow/util/math_test.cc b/cpp/src/arrow/util/math_test.cc
index 2a801a4aa1..df36202d4c 100644
--- a/cpp/src/arrow/util/math_test.cc
+++ b/cpp/src/arrow/util/math_test.cc
@@ -16,6 +16,7 @@
// under the License.
#include <cmath>
+#include <vector>
#include <gtest/gtest.h>
@@ -79,4 +80,36 @@ TEST(NormalPPF, Basics) {
}
}
+TEST(NeumaierSum, Basics) {
+ ASSERT_EQ(0.0, NeumaierSum(std::vector<double>{}));
+ ASSERT_EQ(2.0, NeumaierSum({1., -2., 3.}));
+ ASSERT_EQ(HUGE_VAL, NeumaierSum({1., HUGE_VAL, 1.}));
+ ASSERT_EQ(-HUGE_VAL, NeumaierSum({1., -HUGE_VAL, 1.}));
+ ASSERT_TRUE(std::isnan(NeumaierSum({1., -HUGE_VAL, 1., HUGE_VAL})));
+ ASSERT_TRUE(std::isnan(NeumaierSum({1., NAN, 1.})));
+}
+
+TEST(NeumaierSum, Precision) {
+ // Test cases that would fail with a naive sum
+ ASSERT_EQ(2.0, NeumaierSum({1., 1e100, 1., -1e100}));
+ ASSERT_EQ(2.0, NeumaierSum({1., 1., 1e100, -1e100}));
+ ASSERT_EQ(2.0, NeumaierSum({1e100, -1e100, 1., 1.}));
+ ASSERT_EQ(2.0, NeumaierSum({-1e100, 1e100, 1., 1.}));
+ ASSERT_EQ(2.0, NeumaierSum({1., -1e100, 1e100, 1.}));
+ ASSERT_EQ(2.0, NeumaierSum({1., 1e100, -1e100, 1.}));
+
+ // NumPy code:
+ // array = np.arange(321000, dtype='float64')
+ // array -= np.mean(array)
+ // array *= array
+ constexpr int N = 321000;
+ std::vector<double> values(N);
+ for (int i = 0; i < N; ++i) {
+ values[i] = (i - 160499.5) * (i - 160499.5);
+ }
+ // This is the exact result (as computed using np.sum,
+ // math.fsum or fractions.Fraction).
+ ASSERT_EQ(2756346749973250, NeumaierSum(values));
+}
+
} // namespace arrow::internal
diff --git a/docs/source/cpp/compute.rst b/docs/source/cpp/compute.rst
index 750572713c..8825ffebf2 100644
--- a/docs/source/cpp/compute.rst
+++ b/docs/source/cpp/compute.rst
@@ -220,6 +220,8 @@ the input to a single output value.
+--------------------+---------+------------------+------------------------+----------------------------------+-------+
| index | Unary | Any | Scalar Int64 | :struct:`IndexOptions` | \(4) |
+--------------------+---------+------------------+------------------------+----------------------------------+-------+
+| kurtosis | Unary | Numeric | Scalar Float64 | :struct:`SkewOptions` | \(11) |
++--------------------+---------+------------------+------------------------+----------------------------------+-------+
| last | Unary | Numeric, Binary | Scalar Input type | :struct:`ScalarAggregateOptions` | \(3) |
+--------------------+---------+------------------+------------------------+----------------------------------+-------+
| max | Unary | Non-nested types | Scalar Input type | :struct:`ScalarAggregateOptions` | |
@@ -238,6 +240,8 @@ the input to a single output value.
+--------------------+---------+------------------+------------------------+----------------------------------+-------+
| quantile | Unary | Numeric | Scalar Numeric | :struct:`QuantileOptions` | \(10) |
+--------------------+---------+------------------+------------------------+----------------------------------+-------+
+| skew | Unary | Numeric | Scalar Float64 | :struct:`SkewOptions` | \(11) |
++--------------------+---------+------------------+------------------------+----------------------------------+-------+
| stddev | Unary | Numeric | Scalar Float64 | :struct:`VarianceOptions` | \(11) |
+--------------------+---------+------------------+------------------------+----------------------------------+-------+
| sum | Unary | Numeric | Scalar Numeric | :struct:`ScalarAggregateOptions` | \(9) |
@@ -360,6 +364,8 @@ equivalents above and reflects how they are implemented internally.
+-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
| hash_first_last | Unary | Numeric, Binary | Struct | :struct:`ScalarAggregateOptions` | \(11) |
+-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
+| hash_kurtosis | Unary | Numeric | Float64 | :struct:`SkewOptions` | \(9) |
++-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
| hash_last | Unary | Numeric, Binary | Input type | :struct:`ScalarAggregateOptions` | \(11) |
+-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
| hash_list | Unary | Any | List of input type | | \(3) |
@@ -378,6 +384,8 @@ equivalents above and reflects how they are implemented internally.
+-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
| hash_product | Unary | Numeric | Numeric | :struct:`ScalarAggregateOptions` | \(8) |
+-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
+| hash_skew | Unary | Numeric | Float64 | :struct:`SkewOptions` | \(9) |
++-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
| hash_stddev | Unary | Numeric | Float64 | :struct:`VarianceOptions` | \(9) |
+-------------------------+---------+------------------------------------+------------------------+----------------------------------+-----------+
| hash_sum | Unary | Numeric | Numeric | :struct:`ScalarAggregateOptions` | \(8) |
diff --git a/python/pyarrow/_compute.pyx b/python/pyarrow/_compute.pyx
index ee61077aa8..63370c938b 100644
--- a/python/pyarrow/_compute.pyx
+++ b/python/pyarrow/_compute.pyx
@@ -1889,6 +1889,25 @@ class VarianceOptions(_VarianceOptions):
self._set_options(ddof, skip_nulls, min_count)
+cdef class _SkewOptions(FunctionOptions):
+ def _set_options(self, skip_nulls, min_count):
+ self.wrapped.reset(new CSkewOptions(skip_nulls, min_count))
+
+
+class SkewOptions(_SkewOptions):
+ __doc__ = f"""
+ Options for the `skew` and `kurtosis` functions.
+
+ Parameters
+ ----------
+ {_skip_nulls_doc()}
+ {_min_count_doc(default=0)}
+ """
+
+ def __init__(self, *, skip_nulls=True, min_count=0):
+ self._set_options(skip_nulls, min_count)
+
+
cdef class _SplitOptions(FunctionOptions):
def _set_options(self, max_splits, reverse):
self.wrapped.reset(new CSplitOptions(max_splits, reverse))
diff --git a/python/pyarrow/compute.py b/python/pyarrow/compute.py
index e2c17ee61d..8040cf9ff0 100644
--- a/python/pyarrow/compute.py
+++ b/python/pyarrow/compute.py
@@ -67,6 +67,7 @@ from pyarrow._compute import ( # noqa
ScalarAggregateOptions,
SelectKOptions,
SetLookupOptions,
+ SkewOptions,
SliceOptions,
SortOptions,
SplitOptions,
diff --git a/python/pyarrow/includes/libarrow.pxd b/python/pyarrow/includes/libarrow.pxd
index 8e666b114b..c3ddaba88f 100644
--- a/python/pyarrow/includes/libarrow.pxd
+++ b/python/pyarrow/includes/libarrow.pxd
@@ -2621,6 +2621,12 @@ cdef extern from "arrow/compute/api.h" namespace "arrow::compute" nogil:
c_bool skip_nulls
uint32_t min_count
+ cdef cppclass CSkewOptions \
+ "arrow::compute::SkewOptions"(CFunctionOptions):
+ CSkewOptions(c_bool skip_nulls, uint32_t min_count)
+ c_bool skip_nulls
+ uint32_t min_count
+
cdef cppclass CScalarAggregateOptions \
"arrow::compute::ScalarAggregateOptions"(CFunctionOptions):
CScalarAggregateOptions(c_bool skip_nulls, uint32_t min_count)
diff --git a/python/pyarrow/tests/test_compute.py b/python/pyarrow/tests/test_compute.py
index 99b3047d66..8a756a262b 100644
--- a/python/pyarrow/tests/test_compute.py
+++ b/python/pyarrow/tests/test_compute.py
@@ -184,6 +184,7 @@ def test_option_class_equality(request):
pc.ScalarAggregateOptions(),
pc.SelectKOptions(0, sort_keys=[("b", "ascending")]),
pc.SetLookupOptions(pa.array([1])),
+ pc.SkewOptions(min_count=2),
pc.SliceOptions(0, 1, 1),
pc.SortOptions([("dummy", "descending")], null_placement="at_start"),
pc.SplitOptions(),
@@ -441,6 +442,20 @@ def test_variance():
assert pc.variance(data, ddof=1).as_py() == 6.0
+def test_skew():
+ data = [1, 1, None, 2]
+ assert pc.skew(data).as_py() == pytest.approx(0.707106781186548, rel=1e-10)
+ assert pc.skew(data, skip_nulls=False).as_py() is None
+ assert pc.skew(data, min_count=4).as_py() is None
+
+
+def test_kurtosis():
+ data = [1, 1, None, 2]
+ assert pc.kurtosis(data).as_py() == pytest.approx(-1.5, rel=1e-10)
+ assert pc.kurtosis(data, skip_nulls=False).as_py() is None
+ assert pc.kurtosis(data, min_count=4).as_py() is None
+
+
def test_count_substring():
for (ty, offset) in [(pa.string(), pa.int32()),
(pa.large_string(), pa.int64())]: