Line data    Source code

       1             : /*************************************************************************
       2             :  *
       3             :  *  Project
       4             :  *                         _____ _____  __  __ _____
       5             :  *                        / ____|  __ \|  \/  |  __ \
       6             :  *  ___  _ __   ___ _ __ | |  __| |__) | \  / | |__) |
       7             :  * / _ \| '_ \ / _ \ '_ \| | |_ |  ___/| |\/| |  ___/
       8             :  *| (_) | |_) |  __/ | | | |__| | |    | |  | | |
       9             :  * \___/| .__/ \___|_| |_|\_____|_|    |_|  |_|_|
      10             :  *      | |
      11             :  *      |_|
      12             :  *
      13             :  * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
      14             :  *
      15             :  * This software is licensed as described in the file LICENSE, which
      16             :  * you should have received as part of this distribution. The terms
      17             :  * among other details are referenced in the official documentation
      18             :  * seen here : https://akielaries.github.io/openGPMP/ along with
      19             :  * important files seen in this project.
      20             :  *
      21             :  * You may opt to use, copy, modify, merge, publish, distribute
      22             :  * and/or sell copies of the Software, and permit persons to whom
      23             :  * the Software is furnished to do so, under the terms of the
      24             :  * LICENSE file. As this is an Open Source effort, all implementations
      25             :  * must be of the same methodology.
      26             :  *
      27             :  *
      28             :  *
      29             :  * This software is distributed on an AS IS basis, WITHOUT
      30             :  * WARRANTY OF ANY KIND, either express or implied.
      31             :  *
      32             :  ************************************************************************/
      33             : #include <algorithm>
      34             : #include <cmath>
      35             : #include <openGPMP/stats/describe.hpp>
      36             : #include <vector>
      37             : 
      38           0 : double gpmp::stats::Describe::u_stat(const std::vector<double> &sample1,
      39             :                                      const std::vector<double> &sample2) {
      40           0 :     double U = 0;
      41           0 :     for (double x1 : sample1) {
      42           0 :         for (double x2 : sample2) {
      43           0 :             if (x1 < x2) {
      44           0 :                 U++;
      45             :             }
      46             :         }
      47             :     }
      48           0 :     return U;
      49             : }
      50             : 
      51             : // Arithmetic Mean
      52          17 : double gpmp::stats::Describe::mean_arith(const std::vector<double> &data) {
      53          17 :     double sum = 0.0;
      54          83 :     for (const auto &value : data) {
      55          66 :         sum += value;
      56             :     }
      57          17 :     return sum / static_cast<double>(data.size());
      58             : }
      59             : 
      60             : // Geometric Mean
      61           2 : double gpmp::stats::Describe::mean_geo(const std::vector<double> &data) {
      62           2 :     double product = 1.0;
      63           6 :     for (const auto &value : data) {
      64           4 :         product *= value;
      65             :     }
      66           2 :     return std::pow(product, 1.0 / static_cast<double>(data.size()));
      67             : }
      68             : 
      69             : // Cubic Generalized Mean
      70           4 : double gpmp::stats::Describe::mean_cubic(const std::vector<double> &data,
      71             :                                          double p) {
      72           4 :     double sum = 0.0;
      73          15 :     for (const auto &value : data) {
      74          11 :         sum += std::pow(value, p);
      75             :     }
      76           4 :     return std::pow(sum / static_cast<double>(data.size()), 1.0 / p);
      77             : }
      78             : 
      79             : // Power Geometric Mean
      80           1 : double gpmp::stats::Describe::mean_geo_pow(const std::vector<double> &data,
      81             :                                            double p) {
      82           1 :     double product = 1.0;
      83           2 :     for (const auto &value : data) {
      84           1 :         product *= std::pow(value, p);
      85             :     }
      86           1 :     return std::pow(product, 1.0 / static_cast<double>(data.size()));
      87             : }
      88             : 
      89             : // Harmonic Mean
      90           2 : double gpmp::stats::Describe::mean_harmonic(const std::vector<double> &data) {
      91           2 :     double sum = 0.0;
      92           6 :     for (const auto &value : data) {
      93           4 :         sum += 1.0 / value;
      94             :     }
      95           2 :     return static_cast<double>(data.size()) / sum;
      96             : }
      97             : 
      98             : // Heronian Mean
      99           1 : double gpmp::stats::Describe::mean_heronian(const std::vector<double> &data) {
     100           1 :     double product = 1.0;
     101           2 :     for (const auto &value : data) {
     102           1 :         product *= std::sqrt(value);
     103             :     }
     104           1 :     return std::pow(product, 2.0 / static_cast<double>(data.size()));
     105             : }
     106             : 
     107             : // Heinz Mean
     108           0 : double gpmp::stats::Describe::mean_heinz(const std::vector<double> &data) {
     109           0 :     double sum = 0.0;
     110           0 :     for (const auto &value : data) {
     111           0 :         sum += value * std::log(value);
     112             :     }
     113           0 :     return std::exp(sum / static_cast<double>(data.size()));
     114             : }
     115             : 
     116             : // Lehmer Mean
     117           1 : double gpmp::stats::Describe::mean_lehmer(const std::vector<double> &data,
     118             :                                           double p) {
     119           1 :     double sum = 0.0;
     120           2 :     for (const auto &value : data) {
     121           1 :         sum += std::pow(value, p);
     122             :     }
     123           1 :     return sum / static_cast<double>(data.size());
     124             : }
     125             : 
     126             : // Median
     127           2 : double gpmp::stats::Describe::Median(std::vector<double> data) {
     128           2 :     std::sort(data.begin(), data.end());
     129           2 :     size_t size = data.size();
     130           2 :     if (size % 2 == 0) {
     131           1 :         return (data[size / 2 - 1] + data[size / 2]) / 2.0;
     132             :     } else {
     133           1 :         return data[size / 2];
     134             :     }
     135             : }
     136             : 
     137             : // Average Absolute Deviation
     138           2 : double gpmp::stats::Describe::avg_abs_dev(const std::vector<double> &data) {
     139           2 :     double mean = mean_arith(data);
     140           2 :     double sum = 0.0;
     141           7 :     for (const auto &value : data) {
     142           5 :         sum += std::abs(value - mean);
     143             :     }
     144           2 :     return sum / static_cast<double>(data.size());
     145             : }
     146             : 
     147             : // Coefficient of Variation
     148           2 : double gpmp::stats::Describe::var_coeff(const std::vector<double> &data) {
     149           2 :     double mean = mean_arith(data);
     150           2 :     double stddev = stdev(data, mean);
     151           2 :     return (stddev / mean) * 100.0; // Multiply by 100 for percentage
     152             : }
     153             : 
     154             : // Interquartile Range
     155           1 : double gpmp::stats::Describe::iq_range(const std::vector<double> &data) {
     156           1 :     std::vector<double> sortedData = data;
     157           1 :     std::sort(sortedData.begin(), sortedData.end());
     158             : 
     159           1 :     size_t size = sortedData.size();
     160           1 :     size_t lowerIndex = size / 4;
     161           1 :     size_t upperIndex = 3 * size / 4;
     162             : 
     163           2 :     return sortedData[upperIndex] - sortedData[lowerIndex];
     164           1 : }
     165             : 
     166             : // percentile
     167           1 : double gpmp::stats::Describe::percentile(const std::vector<double> &data,
     168             :                                          double percentile) {
     169           1 :     std::vector<double> sortedData = data;
     170           1 :     std::sort(sortedData.begin(), sortedData.end());
     171             : 
     172           1 :     size_t size = sortedData.size();
     173           1 :     size_t index = static_cast<size_t>(percentile * (size - 1));
     174           2 :     return sortedData[index];
     175           1 : }
     176             : 
     177             : // Range
     178           1 : double gpmp::stats::Describe::range(const std::vector<double> &data) {
     179           1 :     auto result = std::minmax_element(data.begin(), data.end());
     180           1 :     return *result.second - *result.first;
     181             : }
     182             : 
     183             : // Standard Deviation
     184           7 : double gpmp::stats::Describe::stdev(const std::vector<double> &data,
     185             :                                     double mean) {
     186           7 :     double sum = 0.0;
     187          32 :     for (const auto &value : data) {
     188          25 :         sum += std::pow(value - mean, 2.0);
     189             :     }
     190           7 :     return std::sqrt(sum / static_cast<double>(data.size()));
     191             : }
     192             : 
     193             : // variance
     194           3 : double gpmp::stats::Describe::variance(const std::vector<double> &data,
     195             :                                        double mean) {
     196           3 :     double sum = 0.0;
     197          13 :     for (const auto &value : data) {
     198          10 :         sum += std::pow(value - mean, 2.0);
     199             :     }
     200           3 :     return sum / static_cast<double>(data.size());
     201             : }
     202             : 
     203             : // central limit theorem
     204           1 : double gpmp::stats::Describe::clt(const std::vector<double> &data,
     205             :                                   int numSamples) {
     206           1 :     double mean = mean_arith(data);
     207           1 :     double stddev = stdev(data, mean);
     208           1 :     return stddev / std::sqrt(static_cast<double>(numSamples));
     209             : }
     210             : 
     211             : // Kurtosis
     212           1 : double gpmp::stats::Describe::kurtosis(const std::vector<double> &data,
     213             :                                        double mean) {
     214           1 :     double sum = 0.0;
     215           6 :     for (const auto &value : data) {
     216           5 :         sum += std::pow(value - mean, 4.0);
     217             :     }
     218           1 :     double var = variance(data, mean);
     219           1 :     return sum / (data.size() * std::pow(var, 2.0)) - 3.0;
     220             : }
     221             : 
     222             : // l-moments (first two)
     223           1 : double gpmp::stats::Describe::lmoment1(const std::vector<double> &data,
     224             :                                        double mean) {
     225           1 :     double sum = 0.0;
     226           6 :     for (const auto &value : data) {
     227           5 :         sum += std::pow(value - mean, 3.0);
     228             :     }
     229           1 :     return sum / data.size();
     230             : }
     231             : 
     232           1 : double gpmp::stats::Describe::lmoment2(const std::vector<double> &data,
     233             :                                        double mean) {
     234           1 :     double sum = 0.0;
     235           6 :     for (const auto &value : data) {
     236           5 :         sum += std::pow(value - mean, 4.0);
     237             :     }
     238           1 :     return sum / data.size();
     239             : }
     240             : 
     241             : // skewness
     242           1 : double gpmp::stats::Describe::skewness(const std::vector<double> &data,
     243             :                                        double mean,
     244             :                                        double stddev) {
     245           1 :     double sum = 0.0;
     246           6 :     for (const auto &value : data) {
     247           5 :         sum += std::pow((value - mean) / stddev, 3.0);
     248             :     }
     249           1 :     return sum / static_cast<double>(data.size());
     250             : }
     251             : 
     252             : std::vector<size_t>
     253           1 : gpmp::stats::Describe::rank_data(const std::vector<double> &data) {
     254           1 :     std::vector<size_t> ranks(data.size());
     255             : 
     256           6 :     for (size_t i = 0; i < data.size(); ++i) {
     257           5 :         size_t rank = 1;
     258          30 :         for (size_t j = 0; j < data.size(); ++j) {
     259          25 :             if (j != i && data[j] < data[i]) {
     260          10 :                 rank++;
     261             :             }
     262             :         }
     263           5 :         ranks[i] = rank;
     264             :     }
     265             : 
     266           1 :     return ranks;
     267             : }
     268             : 
     269           0 : double gpmp::stats::Describe::partial_corr(const std::vector<double> &x,
     270             :                                            const std::vector<double> &y,
     271             :                                            const std::vector<double> &z) {
     272           0 :     double r_xy = ppmc(x, y);
     273           0 :     double r_xz = ppmc(x, z);
     274           0 :     double r_yz = ppmc(y, z);
     275             : 
     276           0 :     return (r_xy - (r_xz * r_yz)) /
     277           0 :            std::sqrt((1.0 - std::pow(r_xz, 2.0)) * (1.0 - std::pow(r_yz, 2.0)));
     278             : }
     279             : 
     280             : // Pearson Product-Moment Correlation
     281           1 : double gpmp::stats::Describe::ppmc(const std::vector<double> &x,
     282             :                                    const std::vector<double> &y) {
     283           1 :     double mean_x = mean_arith(x);
     284           1 :     double mean_y = mean_arith(y);
     285             : 
     286           1 :     double numerator = 0.0;
     287           1 :     double denominator_x = 0.0;
     288           1 :     double denominator_y = 0.0;
     289             : 
     290           6 :     for (size_t i = 0; i < x.size(); ++i) {
     291           5 :         numerator += (x[i] - mean_x) * (y[i] - mean_y);
     292           5 :         denominator_x += std::pow(x[i] - mean_x, 2.0);
     293           5 :         denominator_y += std::pow(y[i] - mean_y, 2.0);
     294             :     }
     295             : 
     296           1 :     return numerator / std::sqrt(denominator_x * denominator_y);
     297             : }
     298             : 
     299             : // Kendall's Tau Rank Correlation
     300           0 : double gpmp::stats::Describe::kendalls_tau(const std::vector<double> &x,
     301             :                                            const std::vector<double> &y) {
     302           0 :     size_t concordant = 0;
     303           0 :     size_t discordant = 0;
     304             : 
     305           0 :     for (size_t i = 0; i < x.size() - 1; ++i) {
     306           0 :         for (size_t j = i + 1; j < x.size(); ++j) {
     307           0 :             if ((x[i] < x[j] && y[i] < y[j]) || (x[i] > x[j] && y[i] > y[j])) {
     308           0 :                 concordant++;
     309           0 :             } else if ((x[i] < x[j] && y[i] > y[j]) ||
     310           0 :                        (x[i] > x[j] && y[i] < y[j])) {
     311           0 :                 discordant++;
     312             :             }
     313             :         }
     314             :     }
     315             : 
     316           0 :     return static_cast<double>(concordant - discordant) /
     317           0 :            std::sqrt(static_cast<double>((concordant + discordant) *
     318           0 :                                          (x.size() * (x.size() - 1)) / 2));
     319             : }
     320             : 
     321             : // Spearman's Rank Correlation
     322           0 : double gpmp::stats::Describe::spearmans_rho(const std::vector<double> &x,
     323             :                                             const std::vector<double> &y) {
     324           0 :     std::vector<size_t> ranks_x = rank_data(x);
     325           0 :     std::vector<size_t> ranks_y = rank_data(y);
     326             : 
     327           0 :     double d_squared = 0.0;
     328           0 :     for (size_t i = 0; i < x.size(); ++i) {
     329           0 :         d_squared += std::pow(ranks_x[i] - ranks_y[i], 2.0);
     330             :     }
     331             : 
     332             :     return 1.0 -
     333           0 :            (6.0 * d_squared) / (x.size() * (std::pow(x.size(), 2.0) - 1.0));
     334           0 : }