jclstatistics.pas

East make Tray Icon in delphi

{**************************************************************************************************}
{                                                                                                  }
{ Project JEDI Code Library (JCL)                                                                  }
{                                                                                                  }
{ The contents of this file are subject to the Mozilla Public License Version 1.1 (the "License"); }
{ you may not use this file except in compliance with the License. You may obtain a copy of the    }
{ License at http://www.mozilla.org/MPL/                                                           }
{                                                                                                  }
{ Software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTY OF   }
{ ANY KIND, either express or implied. See the License for the specific language governing rights  }
{ and limitations under the License.                                                               }
{                                                                                                  }
{ The Original Code is JclStatistics.pas.                                                          }
{                                                                                                  }
{ The Initial Developer of the Original Code is ESB Consultancy.                                   }
{ Portions created by ESB Consultancy are Copyright ESB Consultancy. All rights reserved.          }
{                                                                                                  }
{ Contributors (in alphabetical order):                                                            }
{   ESB Consultancy                                                                                }
{   Fred Hovey                                                                                     }
{   Marcel van Brakel                                                                              }
{   Matthias Thoma                                                                                 }
{   Robert Marquardt                                                                               }
{   Robert Rossmair                                                                                }
{   Petr Vones                                                                                     }
{                                                                                                  }
{**************************************************************************************************}
{                                                                                                  }
{ Various common statistics routines to calculate, for example, the arithmetic mean, geometric     }
{ meanor median of a set of numbers.                                                               }
{                                                                                                  }
{**************************************************************************************************}

// Last modified: $Date: 2005/03/08 08:33:17 $
// For history see end of file

{ TODO : Test cases! }

unit JclStatistics;

{$I jcl.inc}

interface

uses
  JclBase, JclMath;

type
  EJclStatisticsError = class(EJclMathError);

{ Mean functions }

function ArithmeticMean(const X: TDynFloatArray): Float;
function GeometricMean(const X: TDynFloatArray): Float;
function HarmonicMean(const X: TDynFloatArray): Float;
function HeronianMean(const A, B: Float): Float;

{ Miscellanous }

function BinomialCoeff(N, R: Cardinal): Float;
function IsPositiveFloatArray(const X: TDynFloatArray): Boolean;
function MaxFloatArray(const B: TDynFloatArray): Float;
function MaxFloatArrayIndex(const B: TDynFloatArray): Integer;
function Median(const X: TDynFloatArray): Float;
function MedianUnsorted(const X: TDynFloatArray): Float;
function MinFloatArray(const B: TDynFloatArray): Float;
function MinFloatArrayIndex(const B: TDynFloatArray): Integer;
function Permutation(N, R: Cardinal): Float;
function Combinations(N, R: Cardinal): Float;
function SumOfSquares(const X: TDynFloatArray): Float;
function PopulationVariance(const X: TDynFloatArray): Float;
procedure PopulationVarianceAndMean(const X: TDynFloatArray; var Variance, Mean: Float);
function SampleVariance(const X: TDynFloatArray): Float;
procedure SampleVarianceAndMean(const X: TDynFloatArray; var Variance, Mean: Float);
function StdError(const X: TDynFloatArray): Float; overload;
function StdError(const Variance: Float; const SampleSize: Integer): Float; overload;
function SumFloatArray(const B: TDynFloatArray): Float;
function SumSquareDiffFloatArray(const B: TDynFloatArray; Diff: Float): Float;
function SumSquareFloatArray(const B: TDynFloatArray): Float;
function SumPairProductFloatArray(const X, Y: TDynFloatArray): Float;

implementation

uses
  JclLogic, JclResources, JclSysUtils;

//=== Local helpers ==========================================================

function GetDynLength(const X: TDynFloatArray): Integer;
begin
  Result := Length(X);
end;

function GetDynLengthNotNull(const X: TDynFloatArray): Integer;
begin
  Result := Length(X);
  if Result = 0 then
    raise EJclMathError.CreateRes(@RsEmptyArray);
end;

procedure InvalidSampleSize(SampleSize: Integer);
begin
  raise EJclStatisticsError.CreateResFmt(@RsInvalidSampleSize, [SampleSize]);
end;

function GetSampleSize(const Sample: TDynFloatArray; MinValidSize: Integer = 1): Integer;
begin
  Result := Length(Sample);
  if Result < MinValidSize then
    InvalidSampleSize(Result);
end;

//=== Mean Functions =========================================================

function ArithmeticMean(const X: TDynFloatArray): Float;
begin
  Result := SumFloatArray(X) / Length(X);
end;

function GeometricMean(const X: TDynFloatArray): Float;
var
  I, N: Integer;
begin
  N := GetSampleSize(X);
  Result := 1.0;
  for I := 0 to N - 1 do
  begin
    if X[I] <= PrecisionTolerance then
      raise EJclMathError.CreateRes(@RsNonPositiveArray);
    Result := Result * X[I];
  end;
  Result := Power(Result, 1 / N);
end;

function HarmonicMean(const X: TDynFloatArray): Float;
var
  I, N: Integer;
begin
  Result := 0.0;
  N := GetSampleSize(X);
  for I := 0 to N - 1 do
  begin
    if X[I] <= PrecisionTolerance then
      raise EJclMathError.CreateRes(@RsNonPositiveArray);
    Result := Result + 1 / X[I];
  end;
  Result := N / Result;
end;

function HeronianMean(const A, B: Float): Float;
begin
  Assert(A >= 0);
  Assert(B >= 0);
  Result := (A + Sqrt(A * B) + B) / 3;
end;

//=== Miscellanous ===========================================================

function BinomialCoeff(N, R: Cardinal): Float;
var
  I: Integer;
  K: LongWord;
begin
  if (N = 0) or (R > N) or (N > MaxFactorial) then
  begin
    Result := 0.0;
    Exit;
  end;
  Result := 1.0;
  if not ((R = 0) or (R = N)) then
  begin
    if R > N div 2 then
    R := N - R;
    K := 2;
    try
      for I := N - R + 1 to N do
      begin
        Result := Result * I;
        if K <= R then
        begin
          Result := Result / K;
          Inc(K);
        end;
      end;
      Result := Int(Result + 0.5);
    except
      Result := -1.0;
    end;
  end;
end;


function IsPositiveFloatArray(const X: TDynFloatArray): Boolean;
var
  I, N: Integer;
begin
  Result := False;
  N := GetDynLengthNotNull(X);
  for I := 0 to N - 1 do
    if X[I] <= PrecisionTolerance then
      Exit;
  Result := True;
end;

function MaxFloatArray(const B: TDynFloatArray): Float;
var
  I, N: Integer;
begin
  N := GetDynLengthNotNull(B);
  Result := B[0];
  for I := 1 to N - 1 do
    if B[I] > Result then
      Result := B[I];
end;

function MaxFloatArrayIndex(const B: TDynFloatArray): Integer;
var
  I, N: Integer;
  Max: Float;
begin
  Result := 0;
  N := GetDynLengthNotNull(B);
  Max := B[0];
  for I := 1 to N - 1 do
    if B[I] > Max then
    begin
      Max := B[I];
      Result := I;
    end;
end;

// The FloatArray X must be presorted so Median can calculate the correct value.
//            Y_{(n+1)/2}                     if N is odd
// Median = { 1/2 * (Y_{n/2} + Y_{1+(n/2) }   if N is even

function Median(const X: TDynFloatArray): Float;
var
  N: Integer;
begin
  N := GetSampleSize(X);
  if N = 1 then
    Result := X[0]
  else
  if Odd(N) then
    Result := X[N div 2]
  else
    Result := (X[N div 2 - 1] + X[N div 2]) / 2;
end;

function MedianUnsorted(const X: TDynFloatArray): Float;
var
  SortedList: TDynFloatArray;

begin
  // We need to sort the values first
  SortedList := Copy(X);
  // type cast to Pointer for the sake of FPC
  SortDynArray(Pointer(SortedList), SizeOf(Float),DynArrayCompareFloat);

  // and call the median function afterwards
  Result := Median(SortedList);
end;