1 package StatisticsUtil;
   2 #
   3 # File: StatisticsUtil.pm
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2024 Manish Sud. All rights reserved.
   7 #
   8 # This file is part of MayaChemTools.
   9 #
  10 # MayaChemTools is free software; you can redistribute it and/or modify it under
  11 # the terms of the GNU Lesser General Public License as published by the Free
  12 # Software Foundation; either version 3 of the License, or (at your option) any
  13 # later version.
  14 #
  15 # MayaChemTools is distributed in the hope that it will be useful, but without
  16 # any warranty; without even the implied warranty of merchantability of fitness
  17 # for a particular purpose.  See the GNU Lesser General Public License for more
  18 # details.
  19 #
  20 # You should have received a copy of the GNU Lesser General Public License
  21 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
  22 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
  23 # Boston, MA, 02111-1307, USA.
  24 #
  25 
  26 use strict;
  27 use Exporter;
  28 
  29 use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
  30 
  31 @ISA = qw(Exporter);
  32 @EXPORT = qw(Average AverageDeviation Covariance Correlation Euclidean Factorial FactorialDivison GeometricMean Frequency HarmonicMean KLargest KSmallest Kurtosis Maximum Minimum Mean Median Mode PearsonCorrelation Permutations Product Range RSquare Skewness Sum SumOfSquares StandardDeviation StandardDeviationN  StandardError Standardize StandardScores StandardScoresN TrimMean Variance VarianceN);
  33 @EXPORT_OK = qw();
  34 %EXPORT_TAGS = (all  => [@EXPORT, @EXPORT_OK]);
  35 
  36 # Compute the mean of an array of numbers
  37 sub Average {
  38   my($XArrayRef) = @_;
  39   return Mean($XArrayRef);
  40 }
  41 
  42 # Compute the average of the absolute deviation of an array of numbers: SUM( ABS(x[i] - Xmean) ) / n
  43 sub AverageDeviation {
  44   my($XArrayRef) = @_;
  45 
  46   if (!@$XArrayRef) {
  47     return undef;
  48   }
  49   my($AverageDeviation, $Mean, $Value, $SumOfDeviations);
  50 
  51   $AverageDeviation = 0;
  52   $Mean = Mean($XArrayRef);
  53   foreach $Value (@$XArrayRef) {
  54     $SumOfDeviations += abs($Value - $Mean);
  55   }
  56   $AverageDeviation = $SumOfDeviations / @$XArrayRef;
  57 
  58   return $AverageDeviation;
  59 }
  60 
  61 # Compute correlation coefficient between two arrays of numbers
  62 sub Correlation {
  63   my($XArrayRef, $YArrayRef) = @_;
  64   return PearsonCorrelation($XArrayRef, $YArrayRef);
  65 }
  66 
  67 # Compute the covariance between two arrays of numbers: SUM( (x[i] - Xmean) (y[i] - Ymean) ) / n
  68 sub Covariance {
  69   my($XArrayRef, $YArrayRef) = @_;
  70 
  71   if (!(@$XArrayRef && @$YArrayRef && (@$XArrayRef == @$YArrayRef))) {
  72     return undef;
  73   }
  74   my($Covariance, $XMean, $YMean, $Index, $ProductOfDeviations);
  75 
  76   $Covariance = 0;
  77   $XMean = Mean($XArrayRef);
  78   $YMean = Mean($YArrayRef);
  79   $ProductOfDeviations = 0;
  80   for $Index (0 .. $#{@$XArrayRef}) {
  81     $ProductOfDeviations += (($XArrayRef->[$Index] - $XMean) * ($YArrayRef->[$Index] - $YMean));
  82   }
  83   $Covariance = $ProductOfDeviations / @$XArrayRef;
  84   return $Covariance;
  85 }
  86 
  87 # Compute the euclidean distance of an array of numbers: SQRT( SUM( x[i] ** 2) )
  88 sub Euclidean {
  89   my($XArrayRef) = @_;
  90 
  91   if (!@$XArrayRef) {
  92     return undef;
  93   }
  94   my($SumOfSquares);
  95 
  96   $SumOfSquares = SumOfSquares($XArrayRef);
  97 
  98   return sqrt $SumOfSquares;
  99 }
 100 
 101 # Compute factorial of a number...
 102 sub Factorial {
 103   my($Num) = @_;
 104 
 105   return _Factorial($Num, 1);
 106 }
 107 
 108 # Perform factorial division of two numbers...
 109 sub FactorialDivison {
 110   my($Numerator, $Denominator) = @_;
 111 
 112   # Only works for integer numbers...
 113   if ($Numerator <= 0 || ($Numerator != int($Numerator)) ||
 114      $Denominator <= 0 || ($Denominator != int($Denominator)) ) {
 115     return undef;
 116   }
 117   my($LargerNum, $SmallerNum, $Result);
 118   $LargerNum = ($Numerator > $Denominator) ? $Numerator : $Denominator;
 119   $SmallerNum = ($Numerator < $Denominator) ? $Numerator : $Denominator;
 120 
 121   $Result = _Factorial($LargerNum, $SmallerNum);
 122   if ($Numerator < $Denominator) {
 123     $Result = 1/$Result;
 124   }
 125   return $Result;
 126 }
 127 
 128 # Calculate factorial of a number upto a specific limit...
 129 sub _Factorial {
 130   my($Num, $Limit) = @_;
 131 
 132   # Only works for integer numbers...
 133   if ($Num <= 0 || ($Num != int($Num)) || $Limit < 1) {
 134     return undef;
 135   }
 136 
 137   my($Result) = 1;
 138 
 139   while ($Num > $Limit) {
 140     $Result *= $Num;
 141     $Num--;
 142   }
 143   return $Result;
 144 }
 145 
 146 # Generate all possible permuations or a specific permutations of items in an array
 147 # and return a reference to an array containing array references to generated permuations...
 148 #
 149 # This alogrithm is based on the example provided by Mark Jason-Dominus, and is available
 150 # at CPAN as mjd_permute standalone script.
 151 #
 152 sub Permutations {
 153   my(@DataToPermute) = @_;
 154   my($PermutationNum, $NumOfPermutations, @Permutations);
 155 
 156   if (!@DataToPermute) {
 157     return undef;
 158   }
 159 
 160   @Permutations = ();
 161   $NumOfPermutations = Factorial(scalar @DataToPermute);
 162 
 163   for ($PermutationNum = 0; $PermutationNum < $NumOfPermutations; $PermutationNum++) {
 164     my @Permutation = @DataToPermute[_PermutationNumToPermutation($PermutationNum, $#DataToPermute)];
 165     push @Permutations, \@Permutation;
 166   }
 167 
 168   return \@Permutations;
 169 }
 170 
 171 # Generte Nth permutation for a collection of specific size...
 172 #
 173 sub _PermutationNumToPermutation {
 174   my($Num, $Size) = @_;
 175 
 176   return _PatternToPermutation(_PermutationNumToPattern($Num, $Size));
 177 }
 178 
 179 # Generate Nth pattern for a collection of specific size...
 180 #
 181 sub _PermutationNumToPattern {
 182   my($Num, $Size) = @_;
 183   my($Index, @Pattern);
 184 
 185   $Index = 1;
 186 
 187   while ($Index <= $Size + 1) {
 188     push @Pattern, $Num % $Index;
 189     $Num = int($Num/$Index);
 190     $Index++;
 191   }
 192 
 193   return @Pattern;
 194 }
 195 
 196 # Generate permutation of integers from pattern...
 197 #
 198 sub _PatternToPermutation {
 199   my(@Pattern) = @_;
 200   my(@Source, @Permutation);
 201 
 202   @Source = (0 .. $#Pattern);
 203 
 204   while (@Pattern) {
 205     push @Permutation, splice(@Source, (pop @Pattern), 1);
 206   }
 207 
 208   return @Permutation;
 209 }
 210 
 211 # Compute the frequency of occurance of values in an array of numbers. Three different
 212 # invocation methods are supported:
 213 #
 214 # Frequency(\@ArrayRef) : Using the smallest and largest values, group the numbers into
 215 # 10 bins.
 216 #
 217 # Frequency(\@ArrayRef, $NumOfBins) : Using the smallest and largest values, group the
 218 # numbers into specified bins.
 219 #
 220 # Frequency(\@ArrayRef, \@BinRange): Use bin range to goup the values into different bins.
 221 #
 222 # A hash array is returned with keys and values representing range and frequency values respectively.
 223 # The frequency value for a specific key corresponds to all the values which are greater than
 224 # the previous key and less than or equal to the current key. A key value representing maximum value is
 225 # added for generating frequency distribution for specific number of bins, and whenever the maximum
 226 # array value is greater than the maximum specified in bin range, it is also added to bin range.
 227 #
 228 sub Frequency {
 229   my($XArrayRef) = @_;
 230 
 231   if (!@$XArrayRef) {
 232     return undef;
 233   }
 234 
 235   my($BinRange, $NumOfBins, $BinRangeSpecified);
 236 
 237   $BinRangeSpecified = 0;
 238   $NumOfBins = 10;
 239   if (@_ == 2) {
 240     if (ref($_[1]) eq 'ARRAY') {
 241       $BinRange = $_[1];
 242       if (!(@$BinRange && (@$BinRange > 1))) {
 243         return undef;
 244       }
 245       # Make sure the bin range contains values in increasing order...
 246       my($Index1, $Index2);
 247       for $Index1 (0 .. $#{@$BinRange}) {
 248         for $Index2 (($Index1 + 1) .. $#{@$BinRange}) {
 249           if ($BinRange->[$Index1] >= $BinRange->[$Index2]) {
 250             return undef;
 251           }
 252         }
 253       }
 254       $BinRangeSpecified = 1;
 255     }
 256     else {
 257       $NumOfBins = $_[1];
 258       if ($NumOfBins <= 1) {
 259         return undef;
 260       }
 261     }
 262   }
 263 
 264   # Setup range keys...
 265   my(@RangeKeys);
 266   @RangeKeys = ();
 267 
 268   my($MinValue, $MaxValue) = Range($XArrayRef);
 269   if ($BinRangeSpecified) {
 270     push @RangeKeys, @$BinRange;
 271     if ($MaxValue > $RangeKeys[$#RangeKeys]) {
 272       push @RangeKeys, $MaxValue;
 273     }
 274   }
 275   else {
 276     my($MinValue, $MaxValue) = Range($XArrayRef);
 277     my($Interval) = ($MaxValue - $MinValue)/$NumOfBins;
 278     my($KeyValue) = $MinValue + $Interval;
 279     while ($KeyValue  < $MaxValue) {
 280       push @RangeKeys, $KeyValue;
 281       $KeyValue += $Interval;
 282     }
 283     push @RangeKeys, $MaxValue;
 284   }
 285 
 286   #Setup frequency hash array...
 287   my(%FrequencyMap);
 288   %FrequencyMap = ();
 289 
 290   %FrequencyMap = map { $_ => 0 } @RangeKeys;
 291 
 292   # Count values...
 293   my($Key, $Value);
 294 
 295   VALUE: for $Value (@$XArrayRef) {
 296       for $Key (@RangeKeys) {
 297         if ($Value <= $Key) {
 298           $FrequencyMap{$Key} += 1;
 299           next VALUE;
 300         }
 301       }
 302   }
 303   return (%FrequencyMap);
 304 }
 305 
 306 # Compute the geometric mean of an array of numbers: NthROOT( PRODUCT(x[i]) )
 307 sub GeometricMean {
 308   my($XArrayRef) = @_;
 309 
 310   if (!@$XArrayRef) {
 311     return undef;
 312   }
 313   my($Mean, $Product, $Value);
 314   $Product = 1;
 315   foreach $Value (@$XArrayRef) {
 316     if ($Value <= 0 ) {
 317       return undef;
 318     }
 319     $Product *= $Value;
 320   }
 321   $Mean = $Product ** (1 / @$XArrayRef);
 322   return $Mean;
 323 }
 324 
 325 # Compute the harmonic mean of an array of numbers: 1 / ( SUM(1/x[i]) / n )
 326 sub HarmonicMean {
 327   my($XArrayRef) = @_;
 328 
 329   if (!@$XArrayRef) {
 330     return undef;
 331   }
 332   my($Mean, $Sum, $Value);
 333   $Sum = 0;
 334   foreach $Value (@$XArrayRef) {
 335     if ($Value <= 0 ) {
 336       return undef;
 337     }
 338     $Sum += 1/$Value;
 339   }
 340   $Mean = 1/($Sum/@$XArrayRef);
 341   return $Mean;
 342 }
 343 
 344 # Return the k-largest value from an array of numbers
 345 sub KLargest {
 346   my($XArrayRef, $K) = @_;
 347 
 348   if (!(@$XArrayRef && ($K > 0) && ($K <= @$XArrayRef))) {
 349     return undef;
 350   }
 351   my($KLargest, @SortedXArray);
 352   @SortedXArray = sort { $b <=> $a } @$XArrayRef;
 353   $KLargest = $SortedXArray[$K - 1];
 354   return $KLargest;
 355 }
 356 
 357 # Return the k-smallest value from an array of numbers
 358 sub KSmallest {
 359   my($XArrayRef, $K) = @_;
 360 
 361   if (!(@$XArrayRef && ($K > 0) && ($K <= @$XArrayRef))) {
 362     return undef;
 363   }
 364   my($KSmallest, @SortedXArray);
 365   @SortedXArray = sort { $a <=> $b } @$XArrayRef;
 366   $KSmallest = $SortedXArray[$K - 1];
 367   return $KSmallest;
 368 }
 369 
 370 # Compute the kurtosis of an array of numbers:
 371 # [ {n(n + 1)/(n - 1)(n - 2)(n - 3)}  SUM{ ((x[i] - Xmean)/STDDEV)^4 } ] - {3((n - 1)^2)}/{(n - 2)(n-3)}
 372 #
 373 sub Kurtosis {
 374   my($XArrayRef) = @_;
 375 
 376   if (!@$XArrayRef || ((@$XArrayRef - 3) <= 0)) {
 377     return undef;
 378   }
 379   my($Kurtosis, $Mean, $StandardDeviation, $Value);
 380   $Mean = Mean($XArrayRef);
 381   if (!defined $Mean) {
 382     return undef;
 383   }
 384   $StandardDeviation = StandardDeviation($XArrayRef);
 385   if (!(defined $StandardDeviation &&  $StandardDeviation != 0)) {
 386     return undef;
 387   }
 388 
 389   my($SumOfScores, $SampleSize);
 390   $SumOfScores = 0;
 391   for $Value (@$XArrayRef) {
 392     $SumOfScores += (($Value - $Mean)/$StandardDeviation) ** 4;
 393   }
 394   $SampleSize = @$XArrayRef;
 395   $Kurtosis = ((($SampleSize * ($SampleSize + 1))/(($SampleSize - 1) * ($SampleSize - 2) * ($SampleSize - 3))) * $SumOfScores) - ((3 * (($SampleSize - 1) ** 2))/(($SampleSize - 2) * ($SampleSize - 3)));
 396   return $Kurtosis;
 397 }
 398 
 399 # Return the smallest value from an array of numbers
 400 sub Minimum {
 401   my($XArrayRef) = @_;
 402   return KSmallest($XArrayRef, 1);
 403 }
 404 
 405 # Return the largest value from an array of numbers
 406 sub Maximum {
 407   my($XArrayRef) = @_;
 408   return KLargest($XArrayRef, 1);
 409 }
 410 
 411 # Compute the mean of an array of numbers: SUM( x[i] ) / n
 412 sub Mean {
 413   my($XArrayRef) = @_;
 414 
 415   if (!@$XArrayRef) {
 416     return undef;
 417   }
 418   my($Mean, $Sum, $Value);
 419   $Sum = 0;
 420   foreach $Value (@$XArrayRef) {
 421     $Sum += $Value;
 422   }
 423   $Mean = $Sum / @$XArrayRef;
 424   return $Mean;
 425 }
 426 
 427 # Compute the median value of an array of numbers. For an even number array, it's
 428 # the average of two middle values.
 429 #
 430 # For even values of n: Xsorted[(n - 1)/2 + 1]
 431 # For odd values of n: (Xsorted[n/2] + Xsorted[n/2 + 1])/2
 432 #
 433 sub Median {
 434   my($XArrayRef) = @_;
 435 
 436   if (!@$XArrayRef) {
 437     return undef;
 438   }
 439   my($Median, @SortedXArray);
 440   $Median = 0;
 441   @SortedXArray = sort { $a <=> $b } @$XArrayRef;
 442   if (@$XArrayRef % 2) {
 443     my($MidIndex);
 444     $MidIndex = int(@SortedXArray - 1)/2;
 445     $Median = $SortedXArray[$MidIndex];
 446   }
 447   else {
 448     # Even number array...
 449     my($MidPosition);
 450     $MidPosition = int(@SortedXArray / 2);
 451     $Median = ($SortedXArray[$MidPosition - 1] + $SortedXArray[$MidPosition]) / 2;
 452   }
 453   return $Median;
 454 }
 455 
 456 # Return the most frequently occuring value in an array of numbers
 457 sub Mode {
 458   my($XArrayRef) = @_;
 459 
 460   if (!@$XArrayRef) {
 461     return undef;
 462   }
 463   my($Value, %ValueToCountMap, @CountList, @SortedCountList);
 464   %ValueToCountMap = ();
 465   @CountList = ();
 466   @SortedCountList = ();
 467   for $Value (@$XArrayRef) {
 468     if (exists $ValueToCountMap{$Value}) {
 469       $ValueToCountMap{$Value} += 1;
 470     }
 471     else {
 472       $ValueToCountMap{$Value} = 1;
 473     }
 474   }
 475   for $Value (keys %ValueToCountMap) {
 476     push @CountList, $ValueToCountMap{$Value};
 477   }
 478   @SortedCountList = sort { $b <=> $a } @CountList;
 479 
 480   # Make sure the frequency of mode value is greater than one and check for
 481   # multiple modes as well...
 482   #
 483   my($ModeCount, $ModeValue);
 484   $ModeCount = $SortedCountList[0];
 485   if ($ModeCount <= 1) {
 486     return undef;
 487   }
 488   # Get the first mode value...
 489   VALUE: for $Value (keys %ValueToCountMap) {
 490     if ($ValueToCountMap{$Value} == $ModeCount) {
 491       $ModeValue = $Value;
 492       # Set it to zero to skip it next time...
 493       $ValueToCountMap{$Value} = 0;
 494       last VALUE;
 495     }
 496   }
 497 
 498   if (wantarray) {
 499     # Retrieve all the modes...
 500     my(@Modes, $Count);
 501     @Modes = ();
 502     push @Modes, $ModeValue;
 503     for $Count (@SortedCountList) {
 504       if ($Count == $ModeCount) {
 505       VALUE: for $Value (keys %ValueToCountMap) {
 506           if ($ValueToCountMap{$Value} == $ModeCount) {
 507             push @Modes, $Value;
 508             # Set it to zero to skip it next time...
 509             $ValueToCountMap{$Value} = 0;
 510             last VALUE;
 511           }
 512         }
 513       }
 514     }
 515     return sort {$b <=> $a} @Modes;
 516   }
 517   else {
 518     return $ModeValue;
 519   }
 520 }
 521 
 522 
 523 # Compute the Pearson correlation coefficient between two arrays of numbers:
 524 #
 525 #  SUM( (x[i] - Xmean)(y[i] - Ymean) ) / SQRT( SUM( (x[i] - Xmean)^2 )(SUM( (y[i] - Ymean)^2 ))   )
 526 #
 527 # It returns values in the range from -1.0 to 1.0
 528 sub PearsonCorrelation {
 529   my($XArrayRef, $YArrayRef) = @_;
 530 
 531   if (!(@$XArrayRef && @$YArrayRef && (@$XArrayRef == @$YArrayRef))) {
 532     return undef;
 533   }
 534   my($Correlation, $XMean, $YMean, $Index, $XValueDeviation, $YValueDeviation, $SquareOfXDeviations, $SquareOfYDeviations, $ProductOfDeviations);
 535 
 536   $Correlation = 0;
 537   $XMean = Mean($XArrayRef);
 538   $YMean = Mean($YArrayRef);
 539   $ProductOfDeviations = 0; $SquareOfXDeviations = 0; $SquareOfYDeviations = 0;
 540   for $Index (0 .. $#{@$XArrayRef}) {
 541     $XValueDeviation = $XArrayRef->[$Index] - $XMean;
 542     $YValueDeviation = $YArrayRef->[$Index] - $YMean;
 543     $ProductOfDeviations += ($XValueDeviation * $YValueDeviation);
 544     $SquareOfXDeviations += $XValueDeviation ** 2;
 545     $SquareOfYDeviations += $YValueDeviation ** 2;
 546   }
 547   $Correlation = $ProductOfDeviations / sqrt($SquareOfXDeviations *  $SquareOfYDeviations);
 548   return $Correlation;
 549 }
 550 
 551 # Return the smallest and largest values from an array of numbers
 552 sub Range {
 553   my($XArrayRef) = @_;
 554 
 555   if (!@$XArrayRef) {
 556     return (undef, undef);
 557   }
 558   my($Smallest, $Largest, @SortedXArray);
 559   @SortedXArray = sort { $a <=> $b } @$XArrayRef;
 560   $Smallest = $SortedXArray[0];
 561   $Largest = $SortedXArray[$#SortedXArray];
 562   return ($Smallest, $Largest);
 563 }
 564 
 565 # Compute square of the Pearson correlation coefficient between two arrays of numbers.
 566 #
 567 sub RSquare {
 568   my($XArrayRef, $YArrayRef) = @_;
 569   my($RSquare, $Correlation);
 570 
 571   $RSquare = undef;
 572   $Correlation = PearsonCorrelation($XArrayRef, $YArrayRef);
 573   if (defined $Correlation) {
 574     $RSquare = $Correlation ** 2;
 575   }
 576   return $RSquare;
 577 }
 578 
 579 # Compute the skewness of an array of numbers:
 580 #  {n/(n - 1)(n - 2)} SUM{ ((x[i] - Xmean)/STDDEV)^3 }
 581 #
 582 sub Skewness {
 583   my($XArrayRef) = @_;
 584 
 585   if (!@$XArrayRef || ((@$XArrayRef - 2) <= 0)) {
 586     return undef;
 587   }
 588   my($Skewness, $Mean, $StandardDeviation, $Value);
 589   $Mean = Mean($XArrayRef);
 590   if (!defined $Mean) {
 591     return undef;
 592   }
 593   $StandardDeviation = StandardDeviation($XArrayRef);
 594   if (!(defined $StandardDeviation &&  $StandardDeviation != 0)) {
 595     return undef;
 596   }
 597 
 598   my($SumOfScores, $SampleSize);
 599   $SumOfScores = 0;
 600   for $Value (@$XArrayRef) {
 601     $SumOfScores += (($Value - $Mean)/$StandardDeviation) ** 3;
 602   }
 603   $SampleSize = @$XArrayRef;
 604   $Skewness = ($SampleSize/(($SampleSize - 1) * ($SampleSize - 2) )) * $SumOfScores;
 605   return $Skewness;
 606 }
 607 
 608 # Compute the standard deviation of an array of numbers
 609 sub StandardDeviation {
 610   my($XArrayRef) = @_;
 611   return _CalculateStandardDeviation($XArrayRef, 2);
 612 }
 613 
 614 # Compute the standard deviation of an array of numbers representing entire population
 615 sub StandardDeviationN {
 616   my($XArrayRef) = @_;
 617   return _CalculateStandardDeviation($XArrayRef, 1);
 618 }
 619 
 620 # Compute the standard deviation of an array of numbers.
 621 # Mode 1: SQRT ( SUM( (x[i] - mean)^2 ) / n )
 622 # Mode 2: SQRT ( SUM( (x[i] - mean)^2 ) / (n - 1) )
 623 #
 624 sub _CalculateStandardDeviation {
 625   my($XArrayRef, $Mode) = @_;
 626 
 627   if (!@$XArrayRef) {
 628     return undef;
 629   }
 630   my($StandardDeviation, $Value, $SquareOfDeviations, $Mean, $N);
 631 
 632   $StandardDeviation = 0;
 633   $Mean = Mean($XArrayRef);
 634   $SquareOfDeviations = 0;
 635   foreach $Value (@$XArrayRef) {
 636     $SquareOfDeviations += ($Value - $Mean) ** 2;
 637   }
 638   $N = ($Mode == 1) ? @$XArrayRef : (@$XArrayRef - 1);
 639   $StandardDeviation =  sqrt($SquareOfDeviations / $N);
 640 
 641   return $StandardDeviation;
 642 }
 643 
 644 # Compute the standard error using standard deviation and sample size
 645 sub StandardError {
 646   my($StandardDeviation, $Count) = @_;
 647 
 648   if ($Count <= 0) {
 649     return undef;
 650   }
 651   my($StandardError);
 652   $StandardError = $StandardDeviation / sqrt($Count);
 653 
 654   return $StandardError;
 655 }
 656 
 657 # Standardize the value using mean and standard deviation
 658 sub Standardize {
 659   my($Value, $Mean, $StandardDeviation) = @_;
 660 
 661   if ($StandardDeviation <= 0) {
 662     return undef;
 663   }
 664   my($StandardizedValue);
 665   $StandardizedValue = ($Value - $Mean)/$StandardDeviation;
 666 
 667   return $StandardizedValue;
 668 }
 669 
 670 # Compute the standard deviation above the mean for an array of numbers.
 671 sub StandardScores {
 672   my($XArrayRef) = @_;
 673   return _CalculateStandardScores($XArrayRef, 2);
 674 }
 675 
 676 # Compute the standard deviation above the mean for an array of numbers representing entire population
 677 sub StandardScoresN {
 678   my($XArrayRef) = @_;
 679   return _CalculateStandardScores($XArrayRef, 1);
 680 }
 681 
 682 # Compute the standard deviation above the mean for an array of numbers.
 683 # Mode 1:  (x[i] - mean) / n
 684 # Mode 2:  (x[i] - mean) / (n - 1)
 685 #
 686 sub _CalculateStandardScores {
 687   my($XArrayRef, $Mode) = @_;
 688 
 689   if (!@$XArrayRef) {
 690     return undef;
 691   }
 692   my(@StandardScores, $Mean, $StandardDeviation, $Value);
 693 
 694   $Mean = Mean($XArrayRef);
 695   $StandardDeviation = _CalculateStandardDeviation($XArrayRef, $Mode);
 696   if (!(defined($StandardDeviation) && $StandardDeviation > 0)) {
 697     return undef;
 698   }
 699   @StandardScores = ();
 700   for $Value (@$XArrayRef) {
 701     push @StandardScores, ($Value - $Mean)/$StandardDeviation;
 702   }
 703 
 704   return @StandardScores;
 705 }
 706 
 707 # Compute the product of an array of numbers
 708 sub Product {
 709   my($XArrayRef) = @_;
 710 
 711   if (!@$XArrayRef) {
 712     return undef;
 713   }
 714   my($Product, $Value);
 715   $Product = 1;
 716   foreach $Value (@$XArrayRef) {
 717     $Product *= $Value;
 718   }
 719   return $Product;
 720 }
 721 
 722 # Compute the sum of an array of numbers
 723 sub Sum {
 724   my($XArrayRef) = @_;
 725 
 726   if (!@$XArrayRef) {
 727     return undef;
 728   }
 729   my($Sum, $Value);
 730   $Sum = 0;
 731   foreach $Value (@$XArrayRef) {
 732     $Sum += $Value;
 733   }
 734   return $Sum;
 735 }
 736 
 737 # Compute the sum of squares of an array of numbers
 738 sub SumOfSquares {
 739   my($XArrayRef) = @_;
 740 
 741   if (!@$XArrayRef) {
 742     return undef;
 743   }
 744   my($SumOfSquares, $Value);
 745   $SumOfSquares = 0;
 746   foreach $Value (@$XArrayRef) {
 747     $SumOfSquares += $Value ** 2;
 748   }
 749   return $SumOfSquares;
 750 }
 751 
 752 # Compute the mean of an array of numbers by excluding a fraction of
 753 # numbers from the top and bottom of the data set.
 754 sub TrimMean {
 755   my($XArrayRef, $FractionToExclude) = @_;
 756 
 757   if (!(@$XArrayRef && $FractionToExclude > 0 && $FractionToExclude <= 1)) {
 758     return undef;
 759   }
 760   my($NumberToExclude);
 761   $NumberToExclude = int(@$XArrayRef * $FractionToExclude);
 762   $NumberToExclude = ($NumberToExclude % 2) ? ($NumberToExclude - 1) : $NumberToExclude;
 763   if ($NumberToExclude == @$XArrayRef) {
 764     return undef;
 765   }
 766   my($Mean, $Sum, $Index, $FirstIndex, $LastIndex);
 767   $FirstIndex = $NumberToExclude/2;
 768   $LastIndex = @$XArrayRef - ($NumberToExclude/2) - 1;
 769   $Sum = 0;
 770   my(@SortedXArray);
 771   @SortedXArray = sort { $a <=> $b } @$XArrayRef;
 772   for $Index ($FirstIndex .. $LastIndex) {
 773     $Sum += $SortedXArray[$Index];
 774   }
 775   $Mean = $Sum/(@SortedXArray - $NumberToExclude);
 776   return $Mean;
 777 }
 778 
 779 # Compute the variance of an array of numbers
 780 sub Variance {
 781   my($XArrayRef) = @_;
 782   return _CalculateVariance($XArrayRef, 2);
 783 }
 784 
 785 # Compute the variance of an array of numbers representing entire population
 786 sub VarianceN {
 787   my($XArrayRef) = @_;
 788   return _CalculateVariance($XArrayRef, 1);
 789 }
 790 
 791 # Compute the variance of an array of numbers:
 792 # Mode 1: SUM( (x[i] - Xmean)^2  / n )
 793 # Mode 2: SUM( (x[i] - Xmean)^2  / (n - 1) )
 794 #
 795 sub _CalculateVariance {
 796   my($XArrayRef, $Mode) = @_;
 797 
 798   if (!@$XArrayRef) {
 799     return undef;
 800   }
 801   my($Variance, $Value, $SquareOfDeviations, $Mean, $N);
 802 
 803   $Variance = 0;
 804   $Mean = Mean($XArrayRef);
 805   $SquareOfDeviations = 0;
 806   foreach $Value (@$XArrayRef) {
 807     $SquareOfDeviations += ($Value - $Mean) ** 2;
 808   }
 809   $N = ($Mode == 1) ? @$XArrayRef : (@$XArrayRef - 1);
 810   $Variance =  $SquareOfDeviations / $N;
 811 
 812   return $Variance;
 813 }
 814