Merge pull request #5200 from hczhcz/patch-5

Add statistics functions: skewSamp, skewPop, kurtSamp, and kurtPop

dbms/src/AggregateFunctions/AggregateFunctionStatisticsSimple.cpp

@@ -56,6 +56,10 @@ void registerAggregateFunctionsStatisticsSimple(AggregateFunctionFactory & facto
     factory.registerFunction("varPop", createAggregateFunctionStatisticsUnary<AggregateFunctionVarPopSimple>);
     factory.registerFunction("stddevSamp", createAggregateFunctionStatisticsUnary<AggregateFunctionStddevSampSimple>);
     factory.registerFunction("stddevPop", createAggregateFunctionStatisticsUnary<AggregateFunctionStddevPopSimple>);
+    factory.registerFunction("skewSamp", createAggregateFunctionStatisticsUnary<AggregateFunctionSkewSampSimple>);
+    factory.registerFunction("skewPop", createAggregateFunctionStatisticsUnary<AggregateFunctionSkewPopSimple>);
+    factory.registerFunction("kurtSamp", createAggregateFunctionStatisticsUnary<AggregateFunctionKurtSampSimple>);
+    factory.registerFunction("kurtPop", createAggregateFunctionStatisticsUnary<AggregateFunctionKurtPopSimple>);
 
     factory.registerFunction("covarSamp", createAggregateFunctionStatisticsBinary<AggregateFunctionCovarSampSimple>);
     factory.registerFunction("covarPop", createAggregateFunctionStatisticsBinary<AggregateFunctionCovarPopSimple>);

dbms/src/AggregateFunctions/AggregateFunctionStatisticsSimple.h

@@ -32,30 +32,42 @@ namespace DB
 
 namespace ErrorCodes
 {
-    extern const int LOGICAL_ERROR;
     extern const int DECIMAL_OVERFLOW;
 }
 
 
-template <typename T>
+/**
+    Calculating univariate central moments
+    Levels:
+        level 2 (pop & samp): var, stddev
+        level 3: skewness
+        level 4: kurtosis
+    References:
+        https://en.wikipedia.org/wiki/Moment_(mathematics)
+        https://en.wikipedia.org/wiki/Skewness
+        https://en.wikipedia.org/wiki/Kurtosis
+*/
+template <typename T, size_t _level>
 struct VarMoments
 {
-    T m0{};
-    T m1{};
-    T m2{};
+    T m[_level + 1]{};
 
     void add(T x)
     {
-        ++m0;
-        m1 += x;
-        m2 += x * x;
+        ++m[0];
+        m[1] += x;
+        m[2] += x * x;
+        if constexpr (_level >= 3) m[3] += x * x * x;
+        if constexpr (_level >= 4) m[4] += x * x * x * x;
     }
 
     void merge(const VarMoments & rhs)
     {
-        m0 += rhs.m0;
-        m1 += rhs.m1;
-        m2 += rhs.m2;
+        m[0] += rhs.m[0];
+        m[1] += rhs.m[1];
+        m[2] += rhs.m[2];
+        if constexpr (_level >= 3) m[3] += rhs.m[3];
+        if constexpr (_level >= 4) m[4] += rhs.m[4];
     }
 
     void write(WriteBuffer & buf) const
@@ -70,45 +82,90 @@ struct VarMoments
 
     T NO_SANITIZE_UNDEFINED getPopulation() const
     {
-        return (m2 - m1 * m1 / m0) / m0;
+        return (m[2] - m[1] * m[1] / m[0]) / m[0];
     }
 
     T NO_SANITIZE_UNDEFINED getSample() const
     {
-        if (m0 == 0)
+        if (m[0] == 0)
             return std::numeric_limits<T>::quiet_NaN();
-        return (m2 - m1 * m1 / m0) / (m0 - 1);
+        return (m[2] - m[1] * m[1] / m[0]) / (m[0] - 1);
     }
 
-    T get() const { throw Exception("Unexpected call", ErrorCodes::LOGICAL_ERROR); }
+    T NO_SANITIZE_UNDEFINED getMoment3() const
+    {
+        // to avoid accuracy problem
+        if (m[0] == 1)
+            return 0;
+        return (m[3]
+            - (3 * m[2]
+                - 2 * m[1] * m[1] / m[0]
+            ) * m[1] / m[0]
+        ) / m[0];
+    }
+
+    T NO_SANITIZE_UNDEFINED getMoment4() const
+    {
+        // to avoid accuracy problem
+        if (m[0] == 1)
+            return 0;
+        return (m[4]
+            - (4 * m[3]
+                - (6 * m[2]
+                    - 3 * m[1] * m[1] / m[0]
+                ) * m[1] / m[0]
+            ) * m[1] / m[0]
+        ) / m[0];
+    }
 };
 
-template <typename T>
+template <typename T, size_t _level>
 struct VarMomentsDecimal
 {
     using NativeType = typename T::NativeType;
 
     UInt64 m0{};
-    NativeType m1{};
-    NativeType m2{};
+    NativeType m[_level]{};
+
+    NativeType & getM(size_t i)
+    {
+        return m[i - 1];
+    }
+
+    const NativeType & getM(size_t i) const
+    {
+        return m[i - 1];
+    }
 
     void add(NativeType x)
     {
         ++m0;
-        m1 += x;
+        getM(1) += x;
 
-        NativeType tmp; /// scale' = 2 * scale
-        if (common::mulOverflow(x, x, tmp) || common::addOverflow(m2, tmp, m2))
+        NativeType tmp;
+        if (common::mulOverflow(x, x, tmp) || common::addOverflow(getM(2), tmp, getM(2)))
             throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
+        if constexpr (_level >= 3)
+            if (common::mulOverflow(tmp, x, tmp) || common::addOverflow(getM(3), tmp, getM(3)))
+                throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
+        if constexpr (_level >= 4)
+            if (common::mulOverflow(tmp, x, tmp) || common::addOverflow(getM(4), tmp, getM(4)))
+                throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
     }
 
     void merge(const VarMomentsDecimal & rhs)
     {
         m0 += rhs.m0;
-        m1 += rhs.m1;
+        getM(1) += rhs.getM(1);
 
-        if (common::addOverflow(m2, rhs.m2, m2))
+        if (common::addOverflow(getM(2), rhs.getM(2), getM(2)))
             throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
+        if constexpr (_level >= 3)
+            if (common::addOverflow(getM(3), rhs.getM(3), getM(3)))
+                throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
+        if constexpr (_level >= 4)
+            if (common::addOverflow(getM(4), rhs.getM(4), getM(4)))
+                throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
     }
 
     void write(WriteBuffer & buf) const { writePODBinary(*this, buf); }
@@ -120,8 +177,8 @@ struct VarMomentsDecimal
             return std::numeric_limits<Float64>::infinity();
 
         NativeType tmp;
-        if (common::mulOverflow(m1, m1, tmp) ||
-            common::subOverflow(m2, NativeType(tmp / m0), tmp))
+        if (common::mulOverflow(getM(1), getM(1), tmp) ||
+            common::subOverflow(getM(2), NativeType(tmp / m0), tmp))
             throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
         return convertFromDecimal<DataTypeDecimal<T>, DataTypeNumber<Float64>>(tmp / m0, scale);
     }
@@ -134,15 +191,50 @@ struct VarMomentsDecimal
             return std::numeric_limits<Float64>::infinity();
 
         NativeType tmp;
-        if (common::mulOverflow(m1, m1, tmp) ||
-            common::subOverflow(m2, NativeType(tmp / m0), tmp))
+        if (common::mulOverflow(getM(1), getM(1), tmp) ||
+            common::subOverflow(getM(2), NativeType(tmp / m0), tmp))
             throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
         return convertFromDecimal<DataTypeDecimal<T>, DataTypeNumber<Float64>>(tmp / (m0 - 1), scale);
     }
 
-    Float64 get() const { throw Exception("Unexpected call", ErrorCodes::LOGICAL_ERROR); }
+    Float64 getMoment3(UInt32 scale) const
+    {
+        if (m0 == 0)
+            return std::numeric_limits<Float64>::infinity();
+
+        NativeType tmp;
+        if (common::mulOverflow(2 * getM(1), getM(1), tmp) ||
+            common::subOverflow(3 * getM(2), NativeType(tmp / m0), tmp) ||
+            common::mulOverflow(tmp, getM(1), tmp) ||
+            common::subOverflow(getM(3), NativeType(tmp / m0), tmp))
+            throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
+        return convertFromDecimal<DataTypeDecimal<T>, DataTypeNumber<Float64>>(tmp / m0, scale);
+    }
+
+    Float64 getMoment4(UInt32 scale) const
+    {
+        if (m0 == 0)
+            return std::numeric_limits<Float64>::infinity();
+
+        NativeType tmp;
+        if (common::mulOverflow(3 * getM(1), getM(1), tmp) ||
+            common::subOverflow(6 * getM(2), NativeType(tmp / m0), tmp) ||
+            common::mulOverflow(tmp, getM(1), tmp) ||
+            common::subOverflow(4 * getM(3), NativeType(tmp / m0), tmp) ||
+            common::mulOverflow(tmp, getM(1), tmp) ||
+            common::subOverflow(getM(4), NativeType(tmp / m0), tmp))
+            throw Exception("Decimal math overflow", ErrorCodes::DECIMAL_OVERFLOW);
+        return convertFromDecimal<DataTypeDecimal<T>, DataTypeNumber<Float64>>(tmp / m0, scale);
+    }
 };
 
+/**
+    Calculating multivariate central moments
+    Levels:
+        level 2 (pop & samp): covar
+    References:
+        https://en.wikipedia.org/wiki/Moment_(mathematics)
+*/
 template <typename T>
 struct CovarMoments
 {
@@ -188,8 +280,6 @@ struct CovarMoments
             return std::numeric_limits<T>::quiet_NaN();
         return (xy - x1 * y1 / m0) / (m0 - 1);
     }
-
-    T get() const { throw Exception("Unexpected call", ErrorCodes::LOGICAL_ERROR); }
 };
 
 template <typename T>
@@ -236,9 +326,6 @@ struct CorrMoments
     {
         return (m0 * xy - x1 * y1) / sqrt((m0 * x2 - x1 * x1) * (m0 * y2 - y1 * y1));
     }
-
-    T getPopulation() const { throw Exception("Unexpected call", ErrorCodes::LOGICAL_ERROR); }
-    T getSample() const { throw Exception("Unexpected call", ErrorCodes::LOGICAL_ERROR); }
 };
 
 
@@ -246,18 +333,20 @@ enum class StatisticsFunctionKind
 {
     varPop, varSamp,
     stddevPop, stddevSamp,
+    skewPop, skewSamp,
+    kurtPop, kurtSamp,
     covarPop, covarSamp,
     corr
 };
 
 
-template <typename T, StatisticsFunctionKind _kind>
+template <typename T, StatisticsFunctionKind _kind, size_t _level>
 struct StatFuncOneArg
 {
     using Type1 = T;
     using Type2 = T;
     using ResultType = std::conditional_t<std::is_same_v<T, Float32>, Float32, Float64>;
-    using Data = std::conditional_t<IsDecimalNumber<T>, VarMomentsDecimal<Decimal128>, VarMoments<ResultType>>;
+    using Data = std::conditional_t<IsDecimalNumber<T>, VarMomentsDecimal<Decimal128, _level>, VarMoments<ResultType, _level>>;
 
     static constexpr StatisticsFunctionKind kind = _kind;
     static constexpr UInt32 num_args = 1;
@@ -300,17 +389,28 @@ public:
 
     String getName() const override
     {
-        switch (StatFunc::kind)
-        {
-            case StatisticsFunctionKind::varPop: return "varPop";
-            case StatisticsFunctionKind::varSamp: return "varSamp";
-            case StatisticsFunctionKind::stddevPop: return "stddevPop";
-            case StatisticsFunctionKind::stddevSamp: return "stddevSamp";
-            case StatisticsFunctionKind::covarPop: return "covarPop";
-            case StatisticsFunctionKind::covarSamp: return "covarSamp";
-            case StatisticsFunctionKind::corr: return "corr";
-        }
-        __builtin_unreachable();
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::varPop)
+            return "varPop";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::varSamp)
+            return "varSamp";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::stddevPop)
+            return "stddevPop";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::stddevSamp)
+            return "stddevSamp";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::skewPop)
+            return "skewPop";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::skewSamp)
+            return "skewSamp";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::kurtPop)
+            return "kurtPop";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::kurtSamp)
+            return "kurtSamp";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::covarPop)
+            return "covarPop";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::covarSamp)
+            return "covarSamp";
+        if constexpr (StatFunc::kind == StatisticsFunctionKind::corr)
+            return "corr";
     }
 
     DataTypePtr getReturnType() const override
@@ -351,28 +451,103 @@ public:
 
         if constexpr (IsDecimalNumber<T1>)
         {
-            switch (StatFunc::kind)
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::varPop)
+                dst.push_back(data.getPopulation(src_scale * 2));
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::varSamp)
+                dst.push_back(data.getSample(src_scale * 2));
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::stddevPop)
+                dst.push_back(sqrt(data.getPopulation(src_scale * 2)));
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::stddevSamp)
+                dst.push_back(sqrt(data.getSample(src_scale * 2)));
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::skewPop)
             {
-                case StatisticsFunctionKind::varPop: dst.push_back(data.getPopulation(src_scale * 2)); break;
-                case StatisticsFunctionKind::varSamp: dst.push_back(data.getSample(src_scale * 2)); break;
-                case StatisticsFunctionKind::stddevPop: dst.push_back(sqrt(data.getPopulation(src_scale * 2))); break;
-                case StatisticsFunctionKind::stddevSamp: dst.push_back(sqrt(data.getSample(src_scale * 2))); break;
-                default:
-                    __builtin_unreachable();
+                Float64 var_value = data.getPopulation(src_scale * 2);
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment3(src_scale * 3) / pow(var_value, 1.5));
+                else
+                    dst.push_back(std::numeric_limits<Float64>::quiet_NaN());
+            }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::skewSamp)
+            {
+                Float64 var_value = data.getSample(src_scale * 2);
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment3(src_scale * 3) / pow(var_value, 1.5));
+                else
+                    dst.push_back(std::numeric_limits<Float64>::quiet_NaN());
+            }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::kurtPop)
+            {
+                Float64 var_value = data.getPopulation(src_scale * 2);
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment4(src_scale * 4) / pow(var_value, 2));
+                else
+                    dst.push_back(std::numeric_limits<Float64>::quiet_NaN());
+            }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::kurtSamp)
+            {
+                Float64 var_value = data.getSample(src_scale * 2);
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment4(src_scale * 4) / pow(var_value, 2));
+                else
+                    dst.push_back(std::numeric_limits<Float64>::quiet_NaN());
             }
         }
         else
         {
-            switch (StatFunc::kind)
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::varPop)
+                dst.push_back(data.getPopulation());
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::varSamp)
+                dst.push_back(data.getSample());
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::stddevPop)
+                dst.push_back(sqrt(data.getPopulation()));
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::stddevSamp)
+                dst.push_back(sqrt(data.getSample()));
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::skewPop)
             {
-                case StatisticsFunctionKind::varPop: dst.push_back(data.getPopulation()); break;
-                case StatisticsFunctionKind::varSamp: dst.push_back(data.getSample()); break;
-                case StatisticsFunctionKind::stddevPop: dst.push_back(sqrt(data.getPopulation())); break;
-                case StatisticsFunctionKind::stddevSamp: dst.push_back(sqrt(data.getSample())); break;
-                case StatisticsFunctionKind::covarPop: dst.push_back(data.getPopulation()); break;
-                case StatisticsFunctionKind::covarSamp: dst.push_back(data.getSample()); break;
-                case StatisticsFunctionKind::corr: dst.push_back(data.get()); break;
+                ResultType var_value = data.getPopulation();
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment3() / pow(var_value, 1.5));
+                else
+                    dst.push_back(std::numeric_limits<ResultType>::quiet_NaN());
             }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::skewSamp)
+            {
+                ResultType var_value = data.getSample();
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment3() / pow(var_value, 1.5));
+                else
+                    dst.push_back(std::numeric_limits<ResultType>::quiet_NaN());
+            }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::kurtPop)
+            {
+                ResultType var_value = data.getPopulation();
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment4() / pow(var_value, 2));
+                else
+                    dst.push_back(std::numeric_limits<ResultType>::quiet_NaN());
+            }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::kurtSamp)
+            {
+                ResultType var_value = data.getSample();
+
+                if (var_value > 0)
+                    dst.push_back(data.getMoment4() / pow(var_value, 2));
+                else
+                    dst.push_back(std::numeric_limits<ResultType>::quiet_NaN());
+            }
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::covarPop)
+                dst.push_back(data.getPopulation());
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::covarSamp)
+                dst.push_back(data.getSample());
+            if constexpr (StatFunc::kind == StatisticsFunctionKind::corr)
+                dst.push_back(data.get());
         }
     }
 
@@ -383,10 +558,14 @@ private:
 };
 
 
-template <typename T> using AggregateFunctionVarPopSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::varPop>>;
-template <typename T> using AggregateFunctionVarSampSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::varSamp>>;
-template <typename T> using AggregateFunctionStddevPopSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::stddevPop>>;
-template <typename T> using AggregateFunctionStddevSampSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::stddevSamp>>;
+template <typename T> using AggregateFunctionVarPopSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::varPop, 2>>;
+template <typename T> using AggregateFunctionVarSampSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::varSamp, 2>>;
+template <typename T> using AggregateFunctionStddevPopSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::stddevPop, 2>>;
+template <typename T> using AggregateFunctionStddevSampSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::stddevSamp, 2>>;
+template <typename T> using AggregateFunctionSkewPopSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::skewPop, 3>>;
+template <typename T> using AggregateFunctionSkewSampSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::skewSamp, 3>>;
+template <typename T> using AggregateFunctionKurtPopSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::kurtPop, 4>>;
+template <typename T> using AggregateFunctionKurtSampSimple = AggregateFunctionVarianceSimple<StatFuncOneArg<T, StatisticsFunctionKind::kurtSamp, 4>>;
 template <typename T1, typename T2> using AggregateFunctionCovarPopSimple = AggregateFunctionVarianceSimple<StatFuncTwoArg<T1, T2, StatisticsFunctionKind::covarPop>>;
 template <typename T1, typename T2> using AggregateFunctionCovarSampSimple = AggregateFunctionVarianceSimple<StatFuncTwoArg<T1, T2, StatisticsFunctionKind::covarSamp>>;
 template <typename T1, typename T2> using AggregateFunctionCorrSimple = AggregateFunctionVarianceSimple<StatFuncTwoArg<T1, T2, StatisticsFunctionKind::corr>>;

dbms/tests/queries/0_stateless/00181_aggregate_functions_statistics.reference

@@ -5,6 +5,12 @@ nan
 nan
 0
 nan
+nan
+0
+nan
+nan
+0
+nan
 0
 0
 nan
@@ -14,6 +20,12 @@ nan
 nan
 0
 nan
+nan
+0
+nan
+nan
+0
+nan
 0
 0
 nan

dbms/tests/queries/0_stateless/00181_aggregate_functions_statistics.sql

@@ -30,6 +30,41 @@ SELECT
 FROM series
 );
 
+/* skewSamp */
+
+SELECT skewSamp(x_value) FROM (SELECT x_value FROM series LIMIT 0);
+SELECT skewSamp(x_value) FROM (SELECT x_value FROM series LIMIT 1);
+
+SELECT round(abs(res1 - res2), 6) FROM
+(
+SELECT
+    skewSamp(x_value) AS res1,
+    (
+        sum(x_value * x_value * x_value) / count()
+        - 3 * sum(x_value * x_value) / count() * sum(x_value) / count()
+        + 2 * sum(x_value) / count() * sum(x_value) / count() * sum(x_value) / count()
+    ) / pow((sum(x_value * x_value) - ((sum(x_value) * sum(x_value)) / count())) / (count() - 1), 1.5) AS res2
+FROM series
+);
+
+/* kurtSamp */
+
+SELECT kurtSamp(x_value) FROM (SELECT x_value FROM series LIMIT 0);
+SELECT kurtSamp(x_value) FROM (SELECT x_value FROM series LIMIT 1);
+
+SELECT round(abs(res1 - res2), 6) FROM
+(
+SELECT
+    kurtSamp(x_value) AS res1,
+    (
+        sum(x_value * x_value * x_value * x_value) / count()
+        - 4 * sum(x_value * x_value * x_value) / count() * sum(x_value) / count()
+        + 6 * sum(x_value * x_value) / count() * sum(x_value) / count() * sum(x_value) / count()
+        - 3 * sum(x_value) / count() * sum(x_value) / count() * sum(x_value) / count() * sum(x_value) / count()
+    ) / pow((sum(x_value * x_value) - ((sum(x_value) * sum(x_value)) / count())) / (count() - 1), 2) AS res2
+FROM series
+);
+
 /* varPop */
 
 SELECT varPop(x_value) FROM (SELECT x_value FROM series LIMIT 0);
@@ -56,6 +91,41 @@ SELECT
 FROM series
 );
 
+/* skewPop */
+
+SELECT skewPop(x_value) FROM (SELECT x_value FROM series LIMIT 0);
+SELECT skewPop(x_value) FROM (SELECT x_value FROM series LIMIT 1);
+
+SELECT round(abs(res1 - res2), 6) FROM
+(
+SELECT
+    skewPop(x_value) AS res1,
+    (
+        sum(x_value * x_value * x_value) / count()
+        - 3 * sum(x_value * x_value) / count() * sum(x_value) / count()
+        + 2 * sum(x_value) / count() * sum(x_value) / count() * sum(x_value) / count()
+    ) / pow((sum(x_value * x_value) - ((sum(x_value) * sum(x_value)) / count())) / count(), 1.5) AS res2
+FROM series
+);
+
+/* kurtPop */
+
+SELECT kurtPop(x_value) FROM (SELECT x_value FROM series LIMIT 0);
+SELECT kurtPop(x_value) FROM (SELECT x_value FROM series LIMIT 1);
+
+SELECT round(abs(res1 - res2), 6) FROM
+(
+SELECT
+    kurtPop(x_value) AS res1,
+    (
+        sum(x_value * x_value * x_value * x_value) / count()
+        - 4 * sum(x_value * x_value * x_value) / count() * sum(x_value) / count()
+        + 6 * sum(x_value * x_value) / count() * sum(x_value) / count() * sum(x_value) / count()
+        - 3 * sum(x_value) / count() * sum(x_value) / count() * sum(x_value) / count() * sum(x_value) / count()
+    ) / pow((sum(x_value * x_value) - ((sum(x_value) * sum(x_value)) / count())) / count(), 2) AS res2
+FROM series
+);
+
 /* covarSamp */
 
 SELECT covarSamp(x_value, y_value) FROM (SELECT x_value, y_value FROM series LIMIT 0);