stat.cpp

# include <Siv3D.hpp> // OpenSiv3D v0.4.2

//

namespace s3d
{
    namespace Stat
    {
        /// <summary>
        /// 配列の要素数を求めます。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列の要素数
        /// </returns>
        //template <class Iterator>
        //[[nodiscard]] constexpr size_t Count(Iterator first, Iterator last)
        //{
        //    using category_t = typename std::iterator_traits<Iterator>::iterator_category;

        //    if constexpr (std::is_same_v<category_t, std::random_access_iterator_tag>)
        //    {
        //        return last - first;
        //    }

        //    size_t count = 0;
        //    while (first != last)
        //    {
        //        ++count;
        //        ++first;
        //    }
        //    return count;
        //}

        /// <summary>
        /// 数値が入った配列のイテレータから平均値を計算します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の数値の平均値
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Mean(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            using category_t = typename std::iterator_traits<Iterator>::iterator_category;
            
            if constexpr (std::is_same_v<category_t, std::random_access_iterator_tag>)
            {
                value_t sum{};
                size_t count = last - first;

                while (first != last)
                {
                    sum += *first;

                    ++first;
                }
                return sum / count;
            }

            value_t sum{};
            size_t count{};

            while (first != last)
            {
                ++count;
                sum += *first;

                ++first;
            }
            return sum / count;
        }

        /// <summary>
        /// 数値が入った配列のイテレータから分散を計算します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の数値の分散
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Variance(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            const value_t mean = Mean(first, last);
            value_t sum{};
            size_t count = 0;
            while (first != last)
            {
                value_t add = *first - mean;
                add *= add;
                sum += add;

                ++count;
                ++first;
            }
            return sum / count;
        }

        /// <summary>
        /// 数値が入った配列の最大値を出力します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の最大値
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Max(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            return *std::max_element(first, last);
        }

        /// <summary>
        /// 数値が入った配列の最小値を出力します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の最小値
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Min(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            return *std::min_element(first, last);
        }

        /// <summary>
        /// 数値が入った配列の最頻値を出力します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の最頻値
        /// </returns>
        template <class Iterator>
        [[nodiscard]] auto Mode(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            HashTable<value_t, size_t> hash;

            while (first != last)
            {
                auto itr = hash.find(*first);
                if (itr != hash.end())
                {
                    ++itr.value();
                }
                else
                {
                    hash.emplace(*first, 1);
                }

                ++first;
            }

            return std::max_element(hash.begin(), hash.end(),
                [](const auto& a, const auto& b)
                {
                    return a.second < b.second;
                }
            )->first;
        }

        /// <summary>
        /// 数値が入った配列のイテレータから標準偏差を計算します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の数値の標準偏差
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto StandardDeviation(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;
            return static_cast<value_t>(Sqrt(Variance(first, last)));
        }

        /// <summary>
        /// 数値が入った配列のイテレータから尖度を計算します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の数値の尖度
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Kurtosis(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            const value_t mean = Mean(first, last);
            value_t sum{};
            size_t count = 0;
            auto v = Variance(first, last);
            v *= v;

            while (first != last)
            {
                value_t add = *first - mean;
                add *= add * add * add;
                sum += add;

                ++count;
                ++first;
            }
            return ((sum / v) / count) - 3;
        }

        /// <summary>
        /// 数値が入った配列のイテレータから歪度を計算します。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の数値の歪度
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Skewness(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            const value_t mean = Mean(first, last);
            value_t sum{};
            size_t count = 0;
            auto sd = StandardDeviation(first, last);
            sd *= sd * sd;

            while (first != last)
            {
                value_t add = *first - mean;
                add *= add * add;
                sum += add;

                ++count;
                ++first;
            }
            return (sum / sd) / count;
        }

        /// <summary>
        /// 数値が入った配列のイテレータから中央値を求めます。
        /// </summary>
        /// <param name="first">
        /// 開始イテレータ
        /// </param>
        /// <param name="last">
        /// 終端イテレータ
        /// </param>
        /// <returns>
        /// 配列内の数値の中央値
        /// </returns>
        template <class Iterator>
        [[nodiscard]] constexpr auto Median(Iterator first, Iterator last)
        {
            using value_t = typename std::iterator_traits<Iterator>::value_type;

            if (first == last)
            {
                return static_cast<value_t>(0);
            }

            Array<value_t> data(first, last);
            const size_t median_index = data.size() / 2;

            std::nth_element(data.begin(), data.begin() + median_index, data.end());

            if (IsOdd(data.size()))
            {
                return data[median_index];
            }
            else
            {
                return (data[median_index] + *std::max_element(data.begin(), data.begin() + median_index)) / 2;
            }
        }
    }
}

void Test(double a, double b, double epsilon = 1e-5)
{
    if (std::abs(a - b) < epsilon)
    {
        Console << U"OK";
    }
    else
    {
        Console << U"Failed! expectd:{:.5f}, actual:{:.5f}"_fmt(b, a);
    }
}

void Test(double a)
{
    Console << a;
}

void Main()
{
    const Array<double> v = { 1.0, 2.0, 3.0, 4.0 };
    const std::array<float, 3> a = { 1.1f, 2.2f, 3.3f };
    const float b[2] = { 1.01f, 3.03f };
    const Array<double> z;

    Console << U"[Min]";
    {

        Test(Stat::Min(std::begin(v), std::end(v)), 1);
        Test(Stat::Min(std::begin(a), std::end(a)), 1.1);
        Test(Stat::Min(std::begin(b), std::end(b)), 1.01);
        Test(Stat::Min(std::begin(z), std::end(z)), 0);
    }

    Console << U"[Max]";
    {
        Test(Stat::Max(std::begin(v), std::end(v)), 4.0);
        Test(Stat::Max(std::begin(a), std::end(a)), 3.3);
        Test(Stat::Max(std::begin(b), std::end(b)), 3.03);
        Test(Stat::Max(std::begin(z), std::end(z)), 0);
    }

    Console << U"[Mean]";
    {
        Test(Stat::Mean(std::begin(v), std::end(v)), 2.5);
        Test(Stat::Mean(std::begin(a), std::end(a)), 2.2);
        Test(Stat::Mean(std::begin(b), std::end(b)), 2.02);
        Test(Stat::Mean(std::begin(z), std::end(z)), 0);

        std::list<double> c = { 3, 4, 6, 9 };
        Test(Stat::Mean(std::begin(c), std::end(c)), 5.5);
    }

    Console << U"[Mode]";
    {
        Test(Stat::Mode(std::begin(v), std::end(v)), 1.0);
        Test(Stat::Mode(std::begin(a), std::end(a)), 1.1);
        Test(Stat::Mode(std::begin(b), std::end(b)), 1.01);
        Test(Stat::Mode(std::begin(z), std::end(z)), 0);
    }

    Console << U"[Variance]";
    {
        // https://www.wolframalpha.com/input/?i=population+variance+1+2+3+4&lang=en
        // https://www.wolframalpha.com/input/?i=population+variance+1.1+2.2+3.3&lang=en
        // https://www.wolframalpha.com/input/?i=population+variance+1.01+3.03&lang=en
        Test(Stat::Variance(std::begin(v), std::end(v)), 1.25);
        Test(Stat::Variance(std::begin(a), std::end(a)), 0.806667);
        Test(Stat::Variance(std::begin(b), std::end(b)), 1.0201);
        Test(Stat::Variance(std::begin(z), std::end(z)), 0);
    }

    Console << U"[Median]";
    {
        const Array<double> c = { -3.0, -4.0, -5.0, -6.0 };
        const Array<float> f = { 5.f, 4.f, 3.f, 2.f };
        const std::array<float, 3> d = { 4.2f, 4.2f, 4.2f };
        const float e[1] = { 3.14f };

        // https://www.wolframalpha.com/input/?i=median+1+2+3+4&lang=en
        // https://www.wolframalpha.com/input/?i=median+1.1+2.2+3.3&lang=en
        // https://www.wolframalpha.com/input/?i=median+1.01+3.03&lang=en
        // https://www.wolframalpha.com/input/?i=median+%7B-3%2C+-4%2C+-5%2C+-6%7D&lang=en
        // https://www.wolframalpha.com/input/?i=median+%7B4.2%2C+4.2%2C+4.2%7D&lang=en
        // https://www.wolframalpha.com/input/?i=median+%7B3.14%7D&lang=en
        // https://www.wolframalpha.com/input/?i=median+5+4+3+2&lang=en
        Test(Stat::Median(std::begin(v), std::end(v)), 2.5);
        Test(Stat::Median(std::begin(a), std::end(a)), 2.2);
        Test(Stat::Median(std::begin(b), std::end(b)), 2.02);
        Test(Stat::Median(std::begin(c), std::end(c)), -4.5);
        Test(Stat::Median(std::begin(d), std::end(d)), 4.2);
        Test(Stat::Median(std::begin(e), std::end(e)), 3.14);
        Test(Stat::Median(std::begin(f), std::end(f)), 3.5);
        Test(Stat::Median(std::begin(z), std::end(z)), 0);
    }
    Console << U"[Kurtosis]";
    {
        const Array<double> c = { 2.0, 7.0, 4.0, 9.0, 3.0 };
        const std::array<float, 3> d = { 1.1f, 1.1f, 3.3f };
        const float e[2] = { 1.01f, 3.03f };

        //https://www.wolframalpha.com/input/?i=Kurtosis+%7B2.0%2C+7.0%2C+4.0%2C+9.0%2C+3.0%7D&lang=en
        //https://www.wolframalpha.com/input/?i=Kurtosis+%7B1.1%2C+1.1%2C+3.3%7D&lang=en
        //https://www.wolframalpha.com/input/?i=Kurtosis+%7B1.01%2C+3.03%7D&lang=en
        Test(Stat::Kurtosis(std::begin(c), std::end(c)), 1.60035 - 3);
        Test(Stat::Kurtosis(std::begin(d), std::end(d)), 1.5 - 3);
        Test(Stat::Kurtosis(std::begin(e), std::end(e)), 1 - 3);
        Test(Stat::Kurtosis(std::begin(z), std::end(z)), 0);
    }
    Console << U"[Skewness]";
    {
        const Array<double> c = { 2.0, 7.0, 4.0, 9.0, 3.0 };
        const std::array<float, 3> d = { 1.1f, 1.1f, 3.3f };
        const float e[2] = { 1.01f, 3.03f };

        //https://www.wolframalpha.com/input/?i=Skewness+%7B2.0%2C+7.0%2C+4.0%2C+9.0%2C+3.0%7D&lang=en
        //https://www.wolframalpha.com/input/?i=Skewness+%7B1.1%2C+1.1%2C+3.3%7D&lang=en
        //https://www.wolframalpha.com/input/?i=Skewness+%7B1.01%2C+3.03%7D&lang=en
        Test(Stat::Skewness(std::begin(c), std::end(c)), 0.40604);
        Test(Stat::Skewness(std::begin(d), std::end(d)), 0.707107);
        Test(Stat::Skewness(std::begin(e), std::end(e)), -3.23271 * 1e-16);   //-3.23271×10^-16
        Test(Stat::Skewness(std::begin(z), std::end(z)), 0);
    }

    while (System::Update())
    {

    }
}