diff --git a/components/systools/examples/random/exrandom.lpi b/components/systools/examples/random/exrandom.lpi
new file mode 100644
index 000000000..b890f81bf
--- /dev/null
+++ b/components/systools/examples/random/exrandom.lpi
@@ -0,0 +1,84 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/components/systools/examples/random/exrandom.lpr b/components/systools/examples/random/exrandom.lpr
new file mode 100644
index 000000000..bd7f33ca2
--- /dev/null
+++ b/components/systools/examples/random/exrandom.lpr
@@ -0,0 +1,43 @@
+(* ***** BEGIN LICENSE BLOCK *****
+ * Version: MPL 1.1
+ *
+ * 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 TurboPower SysTools
+ *
+ * The Initial Developer of the Original Code is
+ * TurboPower Software
+ *
+ * Portions created by the Initial Developer are Copyright (C) 1996-2002
+ * the Initial Developer. All Rights Reserved.
+ *
+ * Contributor(s):
+ *
+ * ***** END LICENSE BLOCK ***** *)
+
+{$IFDEF FPC}
+ {$mode DELPHI}
+{$ENDIF}
+
+program ExRandom;
+
+uses
+ Interfaces,
+ Forms,
+ ExRndU in 'ExRndU.pas' {Form1};
+
+{$R *.res}
+
+begin
+ Application.Initialize;
+ Application.CreateForm(TForm1, Form1);
+ Application.Run;
+end.
diff --git a/components/systools/examples/random/exrndu.lfm b/components/systools/examples/random/exrndu.lfm
new file mode 100644
index 000000000..c23b931e2
--- /dev/null
+++ b/components/systools/examples/random/exrndu.lfm
@@ -0,0 +1,146 @@
+object Form1: TForm1
+ Left = 192
+ Height = 498
+ Top = 114
+ Width = 462
+ Caption = 'Random Distributions '
+ ClientHeight = 498
+ ClientWidth = 462
+ Color = clBtnFace
+ Font.Color = clWindowText
+ OnCreate = FormCreate
+ OnDestroy = FormDestroy
+ LCLVersion = '1.9.0.0'
+ object imgGraph: TImage
+ Left = 32
+ Height = 250
+ Top = 200
+ Width = 400
+ end
+ object lblPrompt: TLabel
+ Left = 32
+ Height = 15
+ Top = 24
+ Width = 159
+ Caption = 'Select the distribution to view:'
+ ParentColor = False
+ end
+ object lblGraphTitle: TLabel
+ Left = 32
+ Height = 1
+ Top = 184
+ Width = 1
+ ParentColor = False
+ end
+ object lblParms: TLabel
+ Left = 32
+ Height = 15
+ Top = 56
+ Width = 62
+ Caption = 'Parameters:'
+ ParentColor = False
+ end
+ object lblParm1: TLabel
+ Left = 48
+ Height = 1
+ Top = 80
+ Width = 1
+ ParentColor = False
+ end
+ object lblParm2: TLabel
+ Left = 48
+ Height = 1
+ Top = 104
+ Width = 1
+ ParentColor = False
+ end
+ object lblLeft: TLabel
+ Left = 32
+ Height = 1
+ Top = 456
+ Width = 1
+ ParentColor = False
+ end
+ object lblRight: TLabel
+ Left = 429
+ Height = 1
+ Top = 456
+ Width = 1
+ Alignment = taRightJustify
+ ParentColor = False
+ end
+ object lblMaxY: TLabel
+ Left = 8
+ Height = 1
+ Top = 184
+ Width = 1
+ ParentColor = False
+ end
+ object cboDist: TComboBox
+ Left = 200
+ Height = 23
+ Top = 21
+ Width = 145
+ ItemHeight = 15
+ OnChange = cboDistChange
+ Sorted = True
+ Style = csDropDownList
+ TabOrder = 0
+ end
+ object btnGenerate: TButton
+ Left = 32
+ Height = 25
+ Top = 144
+ Width = 145
+ Caption = 'Generate graph'
+ Default = True
+ OnClick = btnGenerateClick
+ TabOrder = 3
+ end
+ object prgGenProgress: TProgressBar
+ Left = 32
+ Height = 16
+ Top = 472
+ Width = 400
+ TabOrder = 6
+ end
+ object edtParm1: TEdit
+ Left = 200
+ Height = 23
+ Top = 74
+ Width = 121
+ MaxLength = 10
+ TabOrder = 1
+ end
+ object edtParm2: TEdit
+ Left = 200
+ Height = 23
+ Top = 98
+ Width = 121
+ MaxLength = 10
+ TabOrder = 2
+ end
+ object updRight: TUpDown
+ Left = 440
+ Height = 24
+ Top = 448
+ Width = 16
+ Min = 1
+ OnClick = updRightClick
+ Position = 1
+ TabOrder = 5
+ Wrap = False
+ end
+ object updLeft: TUpDown
+ Left = 8
+ Height = 24
+ Top = 448
+ Width = 16
+ Max = 0
+ Min = -100
+ OnClick = updLeftClick
+ Position = 0
+ TabOrder = 4
+ Wrap = False
+ end
+end
diff --git a/components/systools/examples/random/exrndu.pas b/components/systools/examples/random/exrndu.pas
new file mode 100644
index 000000000..ba87aa7bb
--- /dev/null
+++ b/components/systools/examples/random/exrndu.pas
@@ -0,0 +1,523 @@
+(* ***** BEGIN LICENSE BLOCK *****
+ * Version: MPL 1.1
+ *
+ * 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 TurboPower SysTools
+ *
+ * The Initial Developer of the Original Code is
+ * TurboPower Software
+ *
+ * Portions created by the Initial Developer are Copyright (C) 1996-2002
+ * the Initial Developer. All Rights Reserved.
+ *
+ * Contributor(s):
+ *
+ * ***** END LICENSE BLOCK ***** *)
+
+{$IFDEF FPC}
+ {$mode DELPHI}
+{$ENDIF}
+
+unit ExRndU;
+
+interface
+
+uses
+ {$IFNDEF FPC}
+ Windows, Messages,
+ {$ENDIF}
+ SysUtils, Classes, Graphics, Controls, Forms,
+ Dialogs, ComCtrls, StdCtrls, ExtCtrls,
+
+ StRandom;
+
+type
+ TGetRandom = function : double of object;
+
+type
+ TForm1 = class(TForm)
+ imgGraph: TImage;
+ cboDist: TComboBox;
+ lblPrompt: TLabel;
+ btnGenerate: TButton;
+ prgGenProgress: TProgressBar;
+ lblGraphTitle: TLabel;
+ lblParms: TLabel;
+ lblParm1: TLabel;
+ lblParm2: TLabel;
+ edtParm1: TEdit;
+ edtParm2: TEdit;
+ lblLeft: TLabel;
+ lblRight: TLabel;
+ updRight: TUpDown;
+ updLeft: TUpDown;
+ lblMaxY: TLabel;
+ procedure btnGenerateClick(Sender: TObject);
+ procedure cboDistChange(Sender: TObject);
+ procedure FormCreate(Sender: TObject);
+ procedure updRightClick(Sender: TObject; Button: TUDBtnType);
+ procedure updLeftClick(Sender: TObject; Button: TUDBtnType);
+ procedure FormDestroy(Sender: TObject);
+ private
+ { Private declarations }
+ public
+ { Public declarations }
+ GraphLeft : double;
+ GraphRight : double;
+ Value1 : double;
+ Value2 : double;
+ PRNG : TStRandomBase;
+ GetRandom : TGetRandom;
+
+ procedure GenerateGraph(aDistInx : integer);
+
+ procedure PrepForBeta;
+ procedure PrepForCauchy;
+ procedure PrepForChiSquared;
+ procedure PrepForErlang;
+ procedure PrepForExponential;
+ procedure PrepForF;
+ procedure PrepForGamma;
+ procedure PrepForLogNormal;
+ procedure PrepForNormal;
+ procedure PrepForT;
+ procedure PrepForUniform;
+ procedure PrepForWeibull;
+
+ function GetBeta : double;
+ function GetCauchy : double;
+ function GetChiSquared : double;
+ function GetErlang : double;
+ function GetExponential : double;
+ function GetF : double;
+ function GetGamma : double;
+ function GetLogNormal : double;
+ function GetNormal : double;
+ function GetT : double;
+ function GetUniform : double;
+ function GetWeibull : double;
+
+ end;
+
+var
+ Form1: TForm1;
+
+implementation
+
+{$R *.lfm}
+
+const
+ DistNames : array [0..11] of string = (
+ 'Beta', 'Cauchy', 'ChiSquared', 'Erlang', 'Exponential',
+ 'F', 'Gamma', 'LogNormal', 'Normal', 'Student''s t',
+ 'Uniform', 'Weibull');
+
+const
+ RandomCount = 1000000;
+
+procedure TForm1.GenerateGraph(aDistInx : integer);
+var
+ Buckets : array[0..400] of integer;
+ i : integer;
+ R : double;
+ Inx : integer;
+ MaxHt : integer;
+ MaxLineFactor : double;
+ GraphWidth : double;
+ OldPercent : integer;
+ NewPercent : integer;
+begin
+ {zero out the buckets}
+ FillChar(Buckets, sizeof(Buckets), 0);
+
+ {calculate random numbers according to distribution, convert to a
+ bucket index, and increment that bucket count}
+ OldPercent := -1;
+ GraphWidth := imgGraph.Width;
+ for i := 1 to RandomCount do begin
+ NewPercent := (i * 100) div RandomCount;
+ if (NewPercent <> OldPercent) then begin
+ prgGenProgress.Position := NewPercent;
+ OldPercent := NewPercent;
+ end;
+ R := GetRandom;
+ if (GraphLeft <= R) and (R <= GraphRight) then begin
+ Inx := trunc((R - GraphLeft) * GraphWidth / (GraphRight - GraphLeft));
+ if (0 <= Inx) and (Inx <= 400) then
+ inc(Buckets[Inx]);
+ end;
+ end;
+
+ {calculate the largest bucket}
+ MaxHt := 1;
+ for i := 0 to 400 do
+ if (MaxHt < Buckets[i]) then
+ MaxHt := Buckets[i];
+
+ {draw the graph}
+ imgGraph.Canvas.Lock;
+ try
+ imgGraph.Canvas.FillRect(Rect(0, 0, imgGraph.Width, imgGraph.Height));
+ MaxLineFactor := imgGraph.Height / MaxHt;
+ imgGraph.Canvas.Pen.Color := clRed;
+ for i := 0 to 400 do begin
+ imgGraph.Canvas.PenPos := Point(i, imgGraph.Height);
+ imgGraph.Canvas.LineTo(i, imgGraph.Height - trunc(Buckets[i] * MaxLineFactor));
+ end;
+ finally
+ imgGraph.Canvas.Unlock;
+ end;
+
+ lblMaxY.Caption := Format('Max: %8.6f', [MaxHt / RandomCount]);
+end;
+
+procedure TForm1.btnGenerateClick(Sender: TObject);
+begin
+ if (edtParm1.Text = '') then
+ Value1 := 0.0
+ else
+ Value1 := StrToFloat(edtParm1.Text);
+ if (edtParm2.Text = '') then
+ Value2 := 0.0
+ else
+ Value2 := StrToFloat(edtParm2.Text);
+ GenerateGraph(cboDist.ItemIndex);
+end;
+
+procedure TForm1.cboDistChange(Sender: TObject);
+begin
+ case cboDist.ItemIndex of
+ 0 : PrepForBeta;
+ 1 : PrepForCauchy;
+ 2 : PrepForChiSquared;
+ 3 : PrepForErlang;
+ 4 : PrepForExponential;
+ 5 : PrepForF;
+ 6 : PrepForGamma;
+ 7 : PrepForLogNormal;
+ 8 : PrepForNormal;
+ 9 : PrepForT;
+ 10: PrepForUniform;
+ 11: PrepForWeibull
+ end;
+ updRightClick(Self, btNext);
+ updLeftClick(Self, btNext);
+ edtParm1.Text := FloatToStr(Value1);
+ edtParm2.Text := FloatToStr(Value2);
+end;
+
+procedure TForm1.PrepForBeta;
+begin
+ lblParm1.Caption := 'Shape 1:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Shape 2:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := 0;
+ updRight.Position := 1;
+ Value1 := 2.0;
+ Value2 := 4.0;
+ GetRandom := GetBeta;
+end;
+
+procedure TForm1.PrepForCauchy;
+begin
+ lblParm1.Caption := '(none)';
+ lblParm1.Visible := true;
+ lblParm2.Visible := false;
+ edtParm1.Visible := false;
+ edtParm1.Enabled := false;
+ edtParm2.Visible := false;
+ edtParm2.Enabled := false;
+ updLeft.Position := -5;
+ updRight.Position := 5;
+ Value1 := 0.0;
+ Value2 := 0.0;
+ GetRandom := GetCauchy;
+end;
+
+procedure TForm1.PrepForChiSquared;
+begin
+ lblParm1.Caption := 'Degrees of freedom:';
+ lblParm1.Visible := true;
+ lblParm2.Visible := false;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := false;
+ edtParm2.Enabled := false;
+ updLeft.Position := 0;
+ updRight.Position := 20;
+ Value1 := 5.0;
+ Value2 := 0.0;
+ GetRandom := GetChiSquared;
+end;
+
+procedure TForm1.PrepForErlang;
+begin
+ lblParm1.Caption := 'Mean:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Order:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := 0;
+ updRight.Position := 5;
+ Value1 := 1.0;
+ Value2 := 4.0;
+ GetRandom := GetErlang;
+end;
+
+procedure TForm1.PrepForExponential;
+begin
+ lblParm1.Caption := 'Mean:';
+ lblParm1.Visible := true;
+ lblParm2.Visible := false;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := false;
+ edtParm2.Enabled := false;
+ updLeft.Position := 0;
+ updRight.Position := 10;
+ Value1 := 1.0;
+ Value2 := 0.0;
+ GetRandom := GetExponential;
+end;
+
+procedure TForm1.PrepForF;
+begin
+ lblParm1.Caption := 'Degrees of freedom 1:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Degrees of freedom 2:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := 0;
+ updRight.Position := 20;
+ Value1 := 10.0;
+ Value2 := 5.0;
+ GetRandom := GetF;
+end;
+
+procedure TForm1.PrepForGamma;
+begin
+ lblParm1.Caption := 'Shape:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Scale:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := 0;
+ updRight.Position := 10;
+ Value1 := 2.0;
+ Value2 := 1.0;
+ GetRandom := GetGamma;
+end;
+
+procedure TForm1.PrepForLogNormal;
+begin
+ lblParm1.Caption := 'Mean:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Standard deviation:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := 0;
+ updRight.Position := 10;
+ Value1 := 0.0;
+ Value2 := 1.0;
+ GetRandom := GetLogNormal;
+end;
+
+procedure TForm1.PrepForNormal;
+begin
+ lblParm1.Caption := 'Mean:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Standard deviation:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := -5;
+ updRight.Position := 5;
+ Value1 := 0.0;
+ Value2 := 1.0;
+ GetRandom := GetNormal;
+end;
+
+procedure TForm1.PrepForT;
+begin
+ lblParm1.Caption := 'Degrees of freedom:';
+ lblParm1.Visible := true;
+ lblParm2.Visible := false;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := false;
+ edtParm2.Enabled := false;
+ updLeft.Position := -10;
+ updRight.Position := 10;
+ Value1 := 10.0;
+ Value2 := 0.0;
+ GetRandom := GetT;
+end;
+
+procedure TForm1.PrepForUniform;
+begin
+ lblParm1.Caption := '(none)';
+ lblParm1.Visible := true;
+ lblParm2.Visible := false;
+ edtParm1.Visible := false;
+ edtParm1.Enabled := false;
+ edtParm2.Visible := false;
+ edtParm2.Enabled := false;
+ updLeft.Position := 0;
+ updRight.Position := 1;
+ Value1 := 0.0;
+ Value2 := 0.0;
+ GetRandom := GetUniform;
+end;
+
+procedure TForm1.PrepForWeibull;
+begin
+ lblParm1.Caption := 'Shape:';
+ lblParm1.Visible := true;
+ lblParm2.Caption := 'Scale:';
+ lblParm2.Visible := true;
+ edtParm1.Visible := true;
+ edtParm1.Enabled := true;
+ edtParm2.Visible := true;
+ edtParm2.Enabled := true;
+ updLeft.Position := 0;
+ updRight.Position := 10;
+ Value1 := 2.0;
+ Value2 := 3.0;
+ GetRandom := GetWeibull;
+end;
+
+function TForm1.GetBeta : double;
+begin
+ Result := PRNG.AsBeta(Value1, Value2)
+end;
+
+function TForm1.GetCauchy : double;
+begin
+ Result := PRNG.AsCauchy
+end;
+
+function TForm1.GetChiSquared : double;
+begin
+ if (Value1 > 65535.0) then
+ raise Exception.Create(
+ 'TForm1.GetChiSquared: the degrees of freedom value 1 is too large for this example program');
+ Result := PRNG.AsChiSquared(trunc(Value1))
+end;
+
+function TForm1.GetErlang : double;
+begin
+ Result := PRNG.AsErlang(Value1, trunc(Value2))
+end;
+
+function TForm1.GetExponential : double;
+begin
+ Result := PRNG.AsExponential(Value1)
+end;
+
+function TForm1.GetF : double;
+begin
+ if (Value1 > 65535.0) then
+ raise Exception.Create(
+ 'TForm1.GetF: the degrees of freedom value 1 is too large for this example program');
+ if (Value2 > 65535.0) then
+ raise Exception.Create(
+ 'TForm1.GetF: the degrees of freedom value 2 is too large for this example program');
+ Result := PRNG.AsF(trunc(Value1), trunc(Value2))
+end;
+
+function TForm1.GetGamma : double;
+begin
+ Result := PRNG.AsGamma(Value1, Value2)
+end;
+
+function TForm1.GetLogNormal : double;
+begin
+ Result := PRNG.AsLogNormal(Value1, Value2)
+end;
+
+function TForm1.GetNormal : double;
+begin
+ Result := PRNG.AsNormal(Value1, Value2)
+end;
+
+function TForm1.GetT : double;
+begin
+ if (Value1 > 65535.0) then
+ raise Exception.Create(
+ 'TForm1.GetT: the degrees of freedom value is too large for this example program');
+ Result := PRNG.AsT(trunc(Value1))
+end;
+
+function TForm1.GetUniform : double;
+begin
+ Result := PRNG.AsFloat
+end;
+
+function TForm1.GetWeibull : double;
+begin
+ Result := PRNG.AsWeibull(Value1, Value2)
+end;
+
+procedure TForm1.FormCreate(Sender: TObject);
+var
+ i : integer;
+ UniformInx : integer;
+begin
+ cboDist.Items.Clear;
+ UniformInx := -1;
+ for i := 0 to high(DistNames) do begin
+ cboDist.Items.Add(DistNames[i]);
+ if (Copy(DistNames[i], 1, 7) = 'Uniform') then
+ UniformInx := i;
+ end;
+ cboDist.ItemIndex := UniformInx;
+ cboDistChange(Self);
+ PRNG := TStRandomSystem.Create(0);
+end;
+
+procedure TForm1.updRightClick(Sender: TObject; Button: TUDBtnType);
+begin
+ lblRight.Caption := IntToStr(updRight.Position);
+ GraphRight := updRight.Position;
+end;
+
+procedure TForm1.updLeftClick(Sender: TObject; Button: TUDBtnType);
+begin
+ lblLeft.Caption := IntToStr(updLeft.Position);
+ GraphLeft := updLeft.Position;
+end;
+
+procedure TForm1.FormDestroy(Sender: TObject);
+begin
+ PRNG.Free;
+end;
+
+end.
diff --git a/components/systools/laz_systools.lpk b/components/systools/laz_systools.lpk
index 0513f5cbe..435012ff9 100644
--- a/components/systools/laz_systools.lpk
+++ b/components/systools/laz_systools.lpk
@@ -16,7 +16,7 @@
-
+
@@ -101,6 +101,14 @@
+
+
+
+
+
+
+
+
diff --git a/components/systools/laz_systools.pas b/components/systools/laz_systools.pas
index 0050e464b..e0f03b538 100644
--- a/components/systools/laz_systools.pas
+++ b/components/systools/laz_systools.pas
@@ -10,7 +10,7 @@ interface
uses
StBarC, StBase, StConst, StBarPN, StStrL, St2DBarC, StDate, StUtils, StCRC,
StHASH, StToHTML, StStrms, StDict, StIniStm, StDecMth, StExpr, StMath,
- StFIN, StDateSt, StMoney;
+ StFIN, StDateSt, StMoney, StRandom, StStat;
implementation
diff --git a/components/systools/source/design/StReg.pas b/components/systools/source/design/StReg.pas
index 68e1b5534..572503752 100644
--- a/components/systools/source/design/StReg.pas
+++ b/components/systools/source/design/StReg.pas
@@ -126,10 +126,12 @@ uses
StRegIni,
StSaturn,
StSort,
+ *)
StStat,
StStrL,
StStrms,
StStrS,
+ (*
StStrW,
StStrZ,
StText,
@@ -141,15 +143,17 @@ uses
StVArr,
StVenus,
{ new units in ver 4: }
+ *)
StIniStm,
+ (*
StMerge,
StSystem,
StTxtDat,
StDecMth,
*)
- StMoney
+ StMoney,
+ StRandom
(*
- StRandom,
StNTLog,
{ !!! StExpEng unit designed to handle problem with initialization }
{ section in C++Builder; should NOT be included in Registration unit }
diff --git a/components/systools/source/run/strandom.pas b/components/systools/source/run/strandom.pas
new file mode 100644
index 000000000..efbd9a88c
--- /dev/null
+++ b/components/systools/source/run/strandom.pas
@@ -0,0 +1,735 @@
+// Upgraded to Delphi 2009: Sebastian Zierer
+
+(* ***** BEGIN LICENSE BLOCK *****
+ * Version: MPL 1.1
+ *
+ * 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 TurboPower SysTools
+ *
+ * The Initial Developer of the Original Code is
+ * TurboPower Software
+ *
+ * Portions created by the Initial Developer are Copyright (C) 1996-2002
+ * the Initial Developer. All Rights Reserved.
+ *
+ * Contributor(s):
+ *
+ * ***** END LICENSE BLOCK ***** *)
+
+{*********************************************************}
+{* SysTools: StRandom.pas 4.04 *}
+{*********************************************************}
+{* SysTools: Classes for random number distributions *}
+{*********************************************************}
+
+{$IFDEF FPC}
+ {$mode DELPHI}
+{$ENDIF}
+
+//{$I StDefine.inc}
+
+unit StRandom;
+
+interface
+
+uses
+ {$IFNDEF FPC}
+ Windows,
+ {$ENDIF}
+ SysUtils, Classes,
+ StBase;
+
+type
+ TStRandomBase = class
+ private
+ protected
+ function rbMarsagliaGamma(aShape : double) : double;
+ function rbMontyPythonNormal : double;
+ public
+ {uniform distributions}
+ function AsFloat : double; virtual; abstract;
+ function AsInt(aUpperLimit : integer) : integer;
+ function AsIntInRange(aLowerLimit : integer;
+ aUpperLimit : integer) : integer;
+
+ {continuous non-uniform distributions}
+ function AsBeta(aShape1, aShape2 : double) : double;
+ function AsCauchy : double;
+ function AsChiSquared(aFreedom : integer) : double;
+ function AsErlang(aMean : double;
+ aOrder : integer) : double;
+ function AsExponential(aMean : double) : double;
+ function AsF(aFreedom1 : integer;
+ aFreedom2 : integer) : double;
+ function AsGamma(aShape : double; aScale : double) : double;
+ function AsLogNormal(aMean : double;
+ aStdDev : double) : double;
+ function AsNormal(aMean : double;
+ aStdDev : double) : double;
+ function AsT(aFreedom : integer) : double;
+ function AsWeibull(aShape : double;
+ aScale : double) : double;
+ end;
+
+ TStRandomSystem = class(TStRandomBase)
+ private
+ FSeed : integer;
+ protected
+ procedure rsSetSeed(aValue : integer);
+ public
+ constructor Create(aSeed : integer);
+ function AsFloat : double; override;
+ property Seed : integer read FSeed write rsSetSeed;
+ end;
+
+ TStRandomCombined = class(TStRandomBase)
+ private
+ FSeed1 : integer;
+ FSeed2 : integer;
+ protected
+ procedure rcSetSeed1(aValue : integer);
+ procedure rcSetSeed2(aValue : integer);
+ public
+ constructor Create(aSeed1, aSeed2 : integer);
+ function AsFloat : double; override;
+ property Seed1 : integer read FSeed1 write rcSetSeed1;
+ property Seed2 : integer read FSeed2 write rcSetSeed2;
+ end;
+
+ TStRandomMother = class(TStRandomBase)
+ private
+ FNminus4 : integer;
+ FNminus3 : integer;
+ FNminus2 : integer;
+ FNminus1 : integer;
+ FC : integer;
+ protected
+ procedure rsSetSeed(aValue : integer);
+ public
+ constructor Create(aSeed : integer);
+ function AsFloat : double; override;
+ property Seed : integer write rsSetSeed;
+ end;
+
+implementation
+
+uses
+ StConst;
+
+var
+ Root2Pi : double;
+ InvRoot2Pi : double;
+ RootLn4 : double;
+ Ln2 : double;
+ MPN_s : double;
+ Ln2MPN_s : double;
+ MPN_sPlus1 : double;
+
+ Mum1 : integer;
+ Mum2 : integer;
+ Mum3 : integer;
+ Mum4 : integer;
+
+{===Helper routines==================================================}
+function GetRandomSeed : integer;
+var
+ Hash : integer;
+ SystemTime: TSystemTime;
+ G : integer;
+begin
+ {start with the tick count}
+ Hash := integer(GetTickCount);
+
+ {get the current time}
+ GetLocalTime(SystemTime);
+
+ {hash in the milliseconds}
+ Hash := (Hash shl 4) + SystemTime.wMilliseconds;
+ G := Hash and longint($F0000000);
+ if (G <> 0) then
+ Hash := (Hash xor (G shr 24)) xor G;
+
+ {hash in the second}
+ Hash := (Hash shl 4) + SystemTime.wSecond;
+ G := Hash and longint($F0000000);
+ if (G <> 0) then
+ Hash := (Hash xor (G shr 24)) xor G;
+
+ {hash in the minute}
+ Hash := (Hash shl 4) + SystemTime.wMinute;
+ G := Hash and longint($F0000000);
+ if (G <> 0) then
+ Hash := (Hash xor (G shr 24)) xor G;
+
+ {hash in the hour}
+ Hash := (Hash shl 3) + SystemTime.wHour;
+ G := Hash and longint($F0000000);
+ if (G <> 0) then
+ Hash := (Hash xor (G shr 24)) xor G;
+
+ {return the hash}
+ Result := Hash;
+end;
+{====================================================================}
+
+
+{===TStRandomBase====================================================}
+function TStRandomBase.AsBeta(aShape1, aShape2 : double) : double;
+var
+ R1, R2 : double;
+begin
+ if not ((aShape1 > 0.0) and (aShape2 > 0.0)) then
+ raise EStPRNGError.Create(stscPRNGBetaShapeS);
+
+ if (aShape2 = 1.0) then begin
+ repeat
+ R1 := AsFloat;
+ until R1 <> 0.0;
+ Result := exp(ln(R1) / aShape1);
+ end
+ else if (aShape1 = 1.0) then begin
+ repeat
+ R1 := AsFloat;
+ until R1 <> 0.0;
+ Result := 1.0 - exp(ln(R1) / aShape1);
+ end
+ else begin
+ R1 := AsGamma(aShape1, 1.0);
+ R2 := AsGamma(aShape2, 1.0);
+ Result := R1 / (R1 + R2);
+ end;
+end;
+{--------}
+function TStRandomBase.AsCauchy : double;
+var
+ x : double;
+ y : double;
+begin
+ repeat
+ repeat
+ x := AsFloat;
+ until (x <> 0.0);
+ y := (AsFloat * 2.0) - 1.0;
+ until sqr(x) + sqr(y) < 1.0;
+ Result := y / x;
+end;
+{--------}
+function TStRandomBase.AsChiSquared(aFreedom : integer) : double;
+begin
+ if not (aFreedom > 0) then
+ raise EStPRNGError.Create(stscPRNGDegFreedomS);
+
+ Result := AsGamma(aFreedom * 0.5, 2.0)
+end;
+{--------}
+function TStRandomBase.AsErlang(aMean : double;
+ aOrder : integer) : double;
+var
+ Product : double;
+ i : integer;
+begin
+ if not (aMean > 0.0) then
+ raise EStPRNGError.Create(stscPRNGMeanS);
+ if not (aOrder > 0) then
+ raise EStPRNGError.Create(stscPRNGErlangOrderS);
+
+ if (aOrder < 10) then begin
+ Product := 1.0;
+ for i := 1 to aOrder do
+ Product := Product * AsFloat;
+ Result := -aMean * ln(Product) / aOrder;
+ end
+ else begin
+ Result := AsGamma(aOrder, aMean);
+ end;
+end;
+{--------}
+function TStRandomBase.AsExponential(aMean : double) : double;
+var
+ R : double;
+begin
+ if not (aMean > 0.0) then
+ raise EStPRNGError.Create(stscPRNGMeanS);
+
+ repeat
+ R := AsFloat;
+ until (R <> 0.0);
+ Result := -aMean * ln(R);
+end;
+{--------}
+function TStRandomBase.AsF(aFreedom1 : integer;
+ aFreedom2 : integer) : double;
+begin
+ Result := (AsChiSquared(aFreedom1) * aFreedom1) /
+ (AsChiSquared(aFreedom2) * aFreedom2);
+end;
+{--------}
+function TStRandomBase.AsGamma(aShape : double; aScale : double) : double;
+var
+ R : double;
+begin
+ if not (aShape > 0.0) then
+ raise EStPRNGError.Create(stscPRNGGammaShapeS);
+ if not (aScale > 0.0) then
+ raise EStPRNGError.Create(stscPRNGGammaScaleS);
+
+ {there are three cases:
+ ..0.0 < shape < 1.0, use Marsaglia's technique of
+ Gamma(shape) = Gamma(shape+1) * uniform^(1/shape)}
+ if (aShape < 1.0) then begin
+ repeat
+ R := AsFloat;
+ until (R <> 0.0);
+ Result := aScale * rbMarsagliaGamma(aShape + 1.0) *
+ exp(ln(R) / aShape);
+ end
+
+ {..shape = 1.0: this is the same as exponential}
+ else if (aShape = 1.0) then begin
+ repeat
+ R := AsFloat;
+ until (R <> 0.0);
+ Result := aScale * -ln(R);
+ end
+
+ {..shape > 1.0: use Marsaglia./Tsang algorithm}
+ else begin
+ Result := aScale * rbMarsagliaGamma(aShape);
+ end;
+end;
+{--------}
+function TStRandomBase.AsInt(aUpperLimit : integer) : integer;
+begin
+ if not (aUpperLimit > 0) then
+ raise EStPRNGError.Create(stscPRNGLimitS);
+
+ Result := Trunc(AsFloat * aUpperLimit);
+end;
+{--------}
+function TStRandomBase.AsIntInRange(aLowerLimit : integer;
+ aUpperLimit : integer) : integer;
+begin
+ if not (aLowerLimit < aUpperLimit) then
+ raise EStPRNGError.Create(stscPRNGUpperLimitS);
+
+ Result := Trunc(AsFloat * (aUpperLimit - aLowerLimit)) + ALowerLimit;
+end;
+{--------}
+function TStRandomBase.AsLogNormal(aMean : double;
+ aStdDev : double) : double;
+begin
+ Result := exp(AsNormal(aMean, aStdDev));
+end;
+{--------}
+function TStRandomBase.AsNormal(aMean : double;
+ aStdDev : double) : double;
+begin
+ if not (aStdDev > 0.0) then
+ raise EStPRNGError.Create(stscPRNGStdDevS);
+
+ Result := (rbMontyPythonNormal * aStdDev) + aMean;
+
+(*** alternative: The Box-Muller transformation
+var
+ R1, R2 : double;
+ RadiusSqrd : double;
+begin
+ {get two random numbers that define a point in the unit circle}
+ repeat
+ R1 := (2.0 * aRandGen.AsFloat) - 1.0;
+ R2 := (2.0 * aRandGen.AsFloat) - 1.0;
+ RadiusSqrd := sqr(R1) + sqr(R2);
+ until (RadiusSqrd < 1.0) and (RadiusSqrd > 0.0);
+
+ {apply Box-Muller transformation}
+ Result := (R1 * sqrt(-2.0 * ln(RadiusSqrd) / RadiusSqrd) * aStdDev)
+ + aMean;
+ ***)
+end;
+{--------}
+function TStRandomBase.AsT(aFreedom : integer) : double;
+begin
+ if not (aFreedom > 0) then
+ raise EStPRNGError.Create(stscPRNGDegFreedomS);
+
+ Result := rbMontyPythonNormal / sqrt(AsChiSquared(aFreedom) / aFreedom);
+end;
+{--------}
+function TStRandomBase.AsWeibull(aShape : double;
+ aScale : double) : double;
+var
+ R : double;
+begin
+ if not (aShape > 0) then
+ raise EStPRNGError.Create(stscPRNGWeibullShapeS);
+ if not (aScale > 0) then
+ raise EStPRNGError.Create(stscPRNGWeibullScaleS);
+
+ repeat
+ R := AsFloat;
+ until (R <> 0.0);
+ Result := exp(ln(-ln(R)) / aShape) * aScale;
+end;
+{--------}
+
+function TStRandomBase.rbMarsagliaGamma(aShape : double) : double;
+var
+ d : double;
+ c : double;
+ x : double;
+ v : double;
+ u : double;
+ Done : boolean;
+begin
+ {Notes: implements the Marsaglia/Tsang method of generating random
+ numbers belonging to the gamma distribution:
+
+ Marsaglia & Tsang, "A Simple Method for Generating Gamma
+ Variables", ACM Transactions on Mathematical Software,
+ Vol. 26, No. 3, September 2000, Pages 363-372
+
+ It is pointless to try and work out what's going on in this
+ routine without reading this paper :-)
+ }
+
+ d := aShape - (1.0 / 3.0);
+ c := 1.0 / sqrt(9.0 * d);
+ Done := false;
+ {$IFDEF SuppressWarnings}
+ v := 0.0;
+ {$ENDIF}
+
+ while not Done do begin
+ repeat
+ x := rbMontyPythonNormal;
+ v := 1.0 + (c * x);
+ until (v > 0.0);
+
+ v := v * v * v;
+ u := AsFloat;
+
+ Done := u < (1.0 - 0.0331 * sqr(sqr(x)));
+
+ if not Done then
+ Done := ln(u) < (0.5 * sqr(x)) + d * (1.0 - v + ln(v))
+ end;
+
+ Result := d * v;
+end;
+{--------}
+function TStRandomBase.rbMontyPythonNormal : double;
+var
+ x : double;
+ y : double;
+ v : double;
+ NonZeroRandom : double;
+begin
+ {Notes: implements the Monty Python method of generating random
+ numbers belonging to the Normal (Gaussian) distribution:
+
+ Marsaglia & Tsang, "The Monty Python Method for Generating
+ Random Variables", ACM Transactions on Mathematical
+ Software, Vol. 24, No. 3, September 1998, Pages 341-350
+
+ It is pointless to try and work out what's going on in this
+ routine without reading this paper :-)
+
+ Some constants:
+ a = sqrt(ln(4))
+ b = sqrt(2 * pi)
+ s = a / (b - a)
+ }
+
+ {step 1: generate a random number x between +/- sqrt(2*Pi) and
+ return it if its absolute value is less than sqrt(ln(4));
+ note that this exit will happen about 47% of the time}
+ x := ((AsFloat * 2.0) - 1.0) * Root2Pi;
+ if (abs(x) < RootLn4) then begin
+ Result := x;
+ Exit;
+ end;
+
+ {step 2a: generate another random number y strictly between 0 and 1}
+ repeat
+ y := AsFloat;
+ until (y <> 0.0);
+
+ {step 2b: the first quadratic pretest avoids ln() calculation
+ calculate v = 2.8658 - |x| * (2.0213 - 0.3605 * |x|)
+ return x if y < v}
+ v := 2.8658 - Abs(x) * (2.0213 - 0.3605 * Abs(x));
+ if (y < v) then begin
+ Result := x;
+ Exit;
+ end;
+
+ {step 2c: the second quadratic pretest again avoids ln() calculation
+ return s * (b - x) if y > v + 0.0506}
+ if (y > v + 0.0506) then begin
+ if (x > 0) then
+ Result := MPN_s * (Root2Pi - x)
+ else
+ Result := -MPN_s * (Root2Pi + x);
+ Exit;
+ end;
+
+ {step 2d: return x if y < f(x) or
+ ln(y) < ln(2) - (0.5 * x * x) }
+ if (ln(y) < (Ln2 - (0.5 * x * x))) then begin
+ Result := x;
+ Exit;
+ end;
+
+ {step 3: translate x to s * (b - x) and return it if y > g(x) or
+ ln(1 + s - y) < ln(2 * s) - (0.5 * x * x) }
+ if (x > 0) then
+ x := MPN_s * (Root2Pi - x)
+ else
+ x := -MPN_s * (Root2Pi + x);
+ if (ln(MPN_sPlus1 - y) < (Ln2MPN_s - (0.5 * x * x))) then begin
+ Result := x;
+ Exit;
+ end;
+
+ {step 4: the iterative process}
+ repeat
+ repeat
+ NonZeroRandom := AsFloat;
+ until (NonZeroRandom <> 0.0);
+ x := -ln(NonZeroRandom) * InvRoot2Pi;
+ repeat
+ NonZeroRandom := AsFloat;
+ until (NonZeroRandom <> 0.0);
+ y := -ln(NonZeroRandom);
+ until (y + y) > (x * x);
+ if (NonZeroRandom < 0.5) then
+ Result := -(Root2Pi + x)
+ else
+ Result := Root2Pi + x;
+end;
+{====================================================================}
+
+
+{===TStRandomSystem==================================================}
+constructor TStRandomSystem.Create(aSeed : integer);
+begin
+ inherited Create;
+ Seed := aSeed;
+end;
+{--------}
+function TStRandomSystem.AsFloat : double;
+var
+ SaveSeed : integer;
+begin
+ SaveSeed := RandSeed;
+ RandSeed := FSeed;
+ Result := System.Random;
+ FSeed := RandSeed;
+ RandSeed := SaveSeed;
+end;
+{--------}
+procedure TStRandomSystem.rsSetSeed(aValue : integer);
+begin
+ if (aValue = 0) then
+ FSeed := GetRandomSeed
+ else
+ FSeed := aValue;
+end;
+{====================================================================}
+
+
+{===TStRandomCombined================================================}
+const
+ m1 = 2147483563;
+ m2 = 2147483399;
+{--------}
+constructor TStRandomCombined.Create(aSeed1, aSeed2 : integer);
+begin
+ inherited Create;
+ Seed1 := aSeed1;
+ if (aSeed1 = 0) and (aSeed2 = 0) then
+ Sleep(10); // a small delay to enable seed to change
+ Seed2 := aSeed2;
+end;
+{--------}
+function TStRandomCombined.AsFloat : double;
+const
+ a1 = 40014;
+ q1 = 53668; {equals m1 div a1}
+ r1 = 12211; {equals m1 mod a1}
+
+ a2 = 40692;
+ q2 = 52774; {equals m2 div a2}
+ r2 = 3791; {equals m2 mod a2}
+
+ OneOverM1 : double = 1.0 / m1;
+var
+ k : longint;
+ Z : longint;
+begin
+ {advance first PRNG}
+ k := FSeed1 div q1;
+ FSeed1 := (a1 * (FSeed1 - (k * q1))) - (k * r1);
+ if (FSeed1 < 0) then
+ inc(FSeed1, m1);
+
+ {advance second PRNG}
+ k := FSeed2 div q2;
+ FSeed2 := (a2 * (FSeed2 - (k * q2))) - (k * r2);
+ if (FSeed2 < 0) then
+ inc(FSeed2, m2);
+
+ {combine the two seeds}
+ Z := FSeed1 - FSeed2;
+ if (Z <= 0) then
+ Z := Z + m1 - 1;
+ Result := Z * OneOverM1;
+end;
+{--------}
+procedure TStRandomCombined.rcSetSeed1(aValue : integer);
+begin
+ if (aValue = 0) then
+ FSeed1 := GetRandomSeed
+ else
+ FSeed1 := aValue;
+end;
+{--------}
+procedure TStRandomCombined.rcSetSeed2(aValue : integer);
+begin
+ if (aValue = 0) then
+ FSeed2 := GetRandomSeed
+ else
+ FSeed2 := aValue;
+end;
+{====================================================================}
+
+
+{===TStRandomMother==================================================}
+constructor TStRandomMother.Create(aSeed : integer);
+begin
+ inherited Create;
+ Seed := aSeed;
+end;
+{--------}
+function TStRandomMother.AsFloat : double;
+const
+ TwoM31 : double = 1.0 / $7FFFFFFF;
+begin
+ asm
+ push esi
+ push edi
+ push ebx
+
+ {get around a compiler bug where it doesn't notice that edx is
+ being changed in the asm code !!! D5 bug}
+ push edx
+
+ {set ebx to point to self}
+ mov ebx, eax
+
+ {multiply X(n-4) by 21111111}
+ mov eax, [ebx].TStRandomMother.FNMinus4
+ mul [Mum1]
+ mov edi, eax
+ mov esi, edx
+
+ {multiply X(n-3) by 1492 (save it in X(n-4) before though)}
+ mov eax, [ebx].TStRandomMother.FNMinus3
+ mov [ebx].TStRandomMother.FNMinus4, eax
+ mul [Mum2]
+ add edi, eax
+ adc esi, edx
+
+ {multiply X(n-2) by 1776 (save it in X(n-3) before though)}
+ mov eax, [ebx].TStRandomMother.FNMinus2
+ mov [ebx].TStRandomMother.FNMinus3, eax
+ mul [Mum3]
+ add edi, eax
+ adc esi, edx
+
+ {multiply X(n-1) by 5115 (save it in X(n-2) before though)}
+ mov eax, [ebx].TStRandomMother.FNMinus1
+ mov [ebx].TStRandomMother.FNMinus2, eax
+ mul [Mum4]
+ add edi, eax
+ adc esi, edx
+
+ {add in the remainder}
+ add edi, [ebx].TStRandomMother.FC
+ adc esi, 0;
+
+ {save the lower 32 bits in X(n-1), the upper into the remainder}
+ mov [ebx].TStRandomMother.FNMinus1, edi
+ mov [ebx].TStRandomMother.FC, esi
+
+ {get around a compiler bug where it doesn't notice that edx was
+ changed in the asm code !!! D5 bug}
+ pop edx
+
+ pop ebx
+ pop edi
+ pop esi
+ end;
+ Result := (FNMinus1 shr 1) * TwoM31;
+end;
+{--------}
+{$IFOPT Q+}
+{note: TStRandomMother.rsSetSeed expressly overflows integers (it's
+ equivalent to calculating mod 2^32), so we have to force
+ overflow checks off}
+{$DEFINE SaveQPlus}
+{$Q-}
+{$ENDIF}
+procedure TStRandomMother.rsSetSeed(aValue : integer);
+begin
+ if (aValue = 0) then
+ aValue := GetRandomSeed;
+ FNminus4 := aValue;
+ {note: the following code uses the generator
+ Xn := (69069 * Xn-1) mod 2^32
+ from D.E.Knuth, The Art of Computer Programming, Vol. 2
+ (second edition), Addison-Wesley, 1981, pp.102}
+ FNminus3 := 69069 * FNminus4;
+ FNminus2 := 69069 * FNminus3;
+ FNminus1 := 69069 * FNminus2;
+ FC := 69069 * FNminus1;
+end;
+{$IFDEF SaveQPlus}
+{$Q+}
+{$ENDIF}
+{====================================================================}
+
+
+{====================================================================}
+procedure CalcConstants;
+begin
+ {for the normal variates}
+ Root2Pi := sqrt(2 * Pi);
+ InvRoot2Pi := 1.0 / Root2Pi;
+ RootLn4 := sqrt(ln(4.0));
+ Ln2 := ln(2.0);
+ MPN_s := RootLn4 / (Root2Pi - RootLn4);
+ Ln2MPN_s := ln(2.0 * MPN_s);
+ MPN_sPlus1 := MPN_s + 1.0;
+
+ Mum1 := 2111111111;
+ Mum2 := 1492;
+ Mum3 := 1776;
+ Mum4 := 5115;
+end;
+{====================================================================}
+
+
+initialization
+ CalcConstants;
+
+end.
diff --git a/components/systools/source/run/ststat.pas b/components/systools/source/run/ststat.pas
new file mode 100644
index 000000000..90388a2d9
--- /dev/null
+++ b/components/systools/source/run/ststat.pas
@@ -0,0 +1,2338 @@
+// Upgraded to Delphi 2009: Sebastian Zierer
+
+(* ***** BEGIN LICENSE BLOCK *****
+ * Version: MPL 1.1
+ *
+ * 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 TurboPower SysTools
+ *
+ * The Initial Developer of the Original Code is
+ * TurboPower Software
+ *
+ * Portions created by the Initial Developer are Copyright (C) 1996-2002
+ * the Initial Developer. All Rights Reserved.
+ *
+ * Contributor(s):
+ *
+ * ***** END LICENSE BLOCK ***** *)
+
+{*********************************************************}
+{* SysTools: StStat.pas 4.04 *}
+{*********************************************************}
+{* SysTools: Statistical math functions modeled on *}
+{* those in Excel *}
+{*********************************************************}
+
+{$IFDEF FPC}
+ {$mode DELPHI}
+{$ENDIF}
+
+//{$I StDefine.inc}
+
+unit StStat;
+
+{ The statistical distribution functions return results in singles }
+{ since the fractional accuracy of these is typically about 3e-7. }
+
+interface
+
+uses
+ Windows,
+ {$IFDEF UseMathUnit}
+ Math,
+ {$ELSE}
+ StMath,
+ {$ENDIF}
+ SysUtils, StConst, StBase;
+
+{AVEDEV}
+function AveDev(const Data: array of Double) : Double;
+function AveDev16(const Data; NData : Integer) : Double;
+ {-Returns the average of the absolute deviations of data points from their
+ mean. AveDev is a measure of the variability in a data set.
+ }
+
+{CONFIDENCE}
+function Confidence(Alpha, StandardDev : Double; Size : LongInt) : Double;
+ {-Returns the confidence interval for a population mean.
+ The confidence interval is a range on either side of a sample mean.
+ Alpha is the significance level used to compute the confidence level.
+ The confidence level equals 100*(1 - Alpha)%, or in other words,
+ an Alpha of 0.05 indicates a 95 percent confidence level.
+ StandardDev is the population standard deviation for the data range
+ and is assumed to be known.
+ Size is the sample Size.
+ }
+
+{CORREL}
+function Correlation(const Data1, Data2 : array of Double) : Double;
+function Correlation16(const Data1, Data2; NData : Integer) : Double;
+ {-Returns the correlation coefficient of the Data1 and Data2 arrays.
+ Use the correlation coefficient to determine the relationship between
+ two data sets.
+ This function also returns the same value as the PEARSON function
+ in Excel.
+ }
+
+{COVAR}
+function Covariance(const Data1, Data2 : array of Double) : Double;
+function Covariance16(const Data1, Data2; NData : Integer) : Double;
+ {-Returns covariance, the average of products of deviations, for the Data1
+ and Data2 arrays. Use covariance to determine the relationship between
+ two data sets.
+ }
+
+{DEVSQ}
+function DevSq(const Data : array of Double) : Double;
+function DevSq16(const Data; NData : Integer) : Double;
+ {-Returns the sum of squares of deviations of data points from their
+ sample mean.
+ }
+
+{FREQUENCY}
+procedure Frequency(const Data : array of Double; const Bins : array of Double;
+ var Counts : array of LongInt);
+procedure Frequency16(const Data; NData : Integer; const Bins; NBins : Integer;
+ var Counts);
+ {-Calculates how often values occur within an array of data,
+ and then returns an array of counts.
+ Data is an array of values for which you want to count frequencies.
+ Bins is an array of intervals into which you want to group the values in
+ Data.
+ Counts is an array into which the frequency counts are returned. Counts
+ must have one more element than does Bins. The first element of Counts
+ has all items less than the first number in Bins. The next element of
+ Counts is the items that fall between Bins[0] and Bins[1]. The last
+ element of Counts has all items greater than the last number in Bins.
+ }
+
+{GEOMEAN}
+function GeometricMean(const Data : array of Double) : Double;
+function GeometricMean16(const Data; NData : Integer) : Double;
+ {-Returns the geometric mean of an array of positive data. The geometric
+ mean is the n'th root of the product of n positive numbers.}
+
+{HARMEAN}
+function HarmonicMean(const Data : array of Double) : Double;
+function HarmonicMean16(const Data; NData : Integer) : Double;
+ {-Returns the harmonic mean of an array of data. The harmonic mean is the
+ reciprocal of the arithmetic mean of reciprocals.}
+
+{LARGE}
+function Largest(const Data : array of Double; K : Integer) : Double;
+function Largest16(const Data; NData : Integer; K : Integer) : Double;
+ {-Returns the K'th largest value in an array of data. You can use this
+ function to select a value based on its relative standing.}
+
+{MEDIAN}
+function Median(const Data : array of Double) : Double;
+function Median16(const Data; NData : Integer) : Double;
+ {-Returns the median of the given numbers. The median is the number in the
+ middle of a set of numbers; that is, half the numbers have values that
+ are greater than the median, and half have values that are less.
+ If there is an even number of numbers in the set, MEDIAN calculates
+ the average of the two numbers in the middle.
+ }
+
+{MODE}
+function Mode(const Data: array of Double) : Double;
+function Mode16(const Data; NData : Integer) : Double;
+ {-Returns the most frequently occurring, or repetitive, value in an array
+ of data. In case of duplicate frequencies it returns the smallest
+ such value.}
+
+{PERCENTILE}
+function Percentile(const Data : array of Double; K : Double) : Double;
+function Percentile16(const Data; NData : Integer; K : Double) : Double;
+ {-Returns the value of the K'th percentile of an array of data.
+ K is the percentile value, a number between 0 and 1. If K is not a
+ multiple of 1/(n-1) where n is the size of Data, Percentile
+ interpolates between the closest bins.
+ }
+
+{PERCENTRANK}
+function PercentRank(const Data : array of Double; X : Double) : Double;
+function PercentRank16(const Data; NData : Integer; X : Double) : Double;
+ {-Returns the percentile position of a value within an array of data.
+ X is a data value. If X is not found within the array, PercentRank
+ interpolates between the closest data points.
+ }
+
+{PERMUT}
+function Permutations(Number, NumberChosen : Integer) : Extended;
+ {-Returns the number of permutations for a given Number of objects that
+ can be selected from Number objects. A permutation is any set or subset
+ of objects or events where internal order is significant.
+ Permutations are different from combinations, for which the internal order
+ is not significant.
+ Number is an Integer that describes the number of objects.
+ NumberChosen is an Integer that describes the number of objects in
+ each permutation.
+ }
+
+{COMBIN}
+function Combinations(Number, NumberChosen : Integer) : Extended;
+ {-Returns the number of combinations for a given Number of objects that
+ can be selected from Number objects. A combination is any set or subset
+ of objects or events where internal order is not significant.
+ Number is an Integer that describes the number of objects.
+ NumberChosen is an Integer that describes the number of objects in
+ each permutation.
+ }
+
+{FACT}
+function Factorial(N : Integer) : Extended;
+ {-Returns N! as a floating point number.
+ Extended is used for range, not accuracy.
+ }
+
+{RANK}
+function Rank(Number : Double; const Data : array of Double;
+ Ascending : Boolean) : Integer;
+function Rank16(Number : Double; const Data; NData : Integer;
+ Ascending : Boolean) : Integer;
+ {-Returns the rank of a number in a list of numbers. If you were to sort
+ a list that contained no duplicates, the rank of the Number would be its
+ position within the sorted list.
+ Number is the number whose rank you want.
+ Data is the list of numbers.
+ If Ascending is True the rank is measured from the beginning of the
+ array; otherwise from the end of the array.
+ If the Number is not found in the array, 0 is returned. Numbers that
+ appear multiple times all have the same rank, but they affect the
+ rank of numbers appearing after them. For example, in a list of
+ Integers, if the number 10 appears twice and has a rank of 5,
+ then 11 would have a rank of 7 (no number would have a rank of 6).
+ Be sure to sort Data before calling this routine if you want an
+ unambiguous ranking.
+ }
+
+{SMALL}
+function Smallest(const Data : array of Double; K : Integer) : Double;
+function Smallest16(const Data; NData : Integer; K : Integer) : Double;
+ {-Returns the K'th smallest value in an array of data. You can use this
+ function to select a value based on its relative standing.}
+
+{TRIMMEAN}
+function TrimMean(const Data : array of Double; Percent : Double) : Double;
+function TrimMean16(const Data; NData : Integer; Percent : Double) : Double;
+ {-Returns the mean of Data after trimming Percent points from the data set.
+ If Percent is 0.2 and there are 20 points in Data, the 2 largest and 2
+ smallest points would be dropped before computing the mean. Percent must
+ be a number between 0 and 1.
+ }
+
+{--------------------------------------------------------------------------}
+
+type
+ {full statistics for a linear regression}
+ TStLinEst = record
+ B0, B1 : Double; {model coefficients}
+ seB0, seB1 : Double; {standard error of model coefficients}
+ R2 : Double; {coefficient of determination}
+ sigma : Double; {standard error of regression}
+ SSr, SSe : Double; {elements for ANOVA table}
+ F0 : Double; {F-statistic to test B1=0}
+ df : Integer; {denominator degrees of freedom for F-statistic}
+ end;
+
+{LINEST}
+procedure LinEst(const KnownY : array of Double;
+ const KnownX : array of Double; var LF : TStLinEst; ErrorStats : Boolean);
+procedure LinEst16(const KnownY; const KnownX; NData : Integer;
+ var LF : TStLinEst; ErrorStats : Boolean);
+ {-Performs linear fit to data and returns coefficients and error
+ statistics.
+ KnownY is the dependent array of known data points.
+ KnownX is the independent array of known data points.
+ NData must be greater than 2.
+ If ErrorStats is FALSE, only B0 and B1 are computed; the other fields
+ of TStLinEst are set to 0.0.
+ See declaration of TStLinEst for returned data.
+
+ }
+
+{LOGEST}
+procedure LogEst(const KnownY : array of Double;
+ const KnownX : array of Double; var LF : TStLinEst; ErrorStats : Boolean);
+procedure LogEst16(const KnownY; const KnownX; NData : Integer;
+ var LF : TStLinEst; ErrorStats : Boolean);
+ {-Performs log-linear fit to data and returns coefficients and error
+ statistics.
+ KnownY is the dependent array of known data points.
+ KnownX is the independent array of known data points.
+ NData must be greater than 2.
+ KnownY is transformed using ln(KnownY) before fitting:
+ y = B0*B1^x implies that ln(y) = ln(B0)+X*ln(B1)
+ In the returned LF, B0 and B1 are returned as shown. Other values in LF
+ are in terms of the log transformation.
+ If ErrorStats is FALSE, only B0 and B1 are computed; the other fields
+ of TStLinEst are set to 0.0.
+ See declaration of TStLinEst for returned data.
+ }
+
+{FORECAST}
+function Forecast(X : Double; const KnownY: array of Double;
+ const KnownX : array of Double) : Double;
+function Forecast16(X : Double; const KnownY; const KnownX;
+ NData : Integer) : Double;
+ {-Calculates a future value by using existing values. The predicted value
+ is a y-value for a given X-value. The known values are existing X-values
+ and y-values, and the new value is predicted by using linear regression.
+ X is the data point for which you want to predict a value.
+ KnownY is the dependent array of known data points.
+ KnownX is the independent array of known data points.
+ }
+
+{similar to GROWTH but more consistent with FORECAST}
+function ForecastExponential(X : Double; const KnownY : array of Double;
+ const KnownX : array of Double) : Double;
+function ForecastExponential16(X : Double; const KnownY; const KnownX;
+ NData : Integer) : Double;
+ {-Calculates a future value by using existing values. The predicted value
+ is a y-value for a given X-value. The known values are existing X-values
+ and y-values, and the new value is predicted by using linear regression
+ to an exponential growth model, y = B0*B1^X.
+ X is the data point for which you want to predict a value.
+ KnownY is the dependent array of known data points.
+ KnownX is the independent array of known data points.
+ }
+
+{INTERCEPT}
+function Intercept(const KnownY : array of Double;
+ const KnownX : array of Double) : Double;
+function Intercept16(const KnownY; const KnownX; NData : Integer) : Double;
+ {-Calculates the point at which a line will intersect the y-axis by using
+ existing X-values and y-values. The intercept point is based on a best-fit
+ regression line plotted through the known X-values and known y-values.
+ Use the intercept when you want to determine the value of the dependent
+ variable when the independent variable is 0 (zero).
+ }
+
+{RSQ}
+function RSquared(const KnownY : array of Double;
+ const KnownX : array of Double) : Double;
+function RSquared16(const KnownY; const KnownX; NData : Integer) : Double;
+ {-Returns the square of the Pearson product moment correlation coefficient
+ through data points in KnownY's and KnownX's. The r-squared value can
+ be interpreted as the proportion of the variance in y attributable to
+ the variance in X.
+ }
+
+{SLOPE}
+function Slope(const KnownY : array of Double;
+ const KnownX : array of Double) : Double;
+function Slope16(const KnownY; const KnownX; NData : Integer) : Double;
+ {-Returns the slope of the linear regression line through data points in
+ KnownY's and KnownX's. The slope is the vertical distance divided by the
+ horizontal distance between any two points on the line, which is the rate
+ of change along the regression line.
+ }
+
+{STEYX}
+function StandardErrorY(const KnownY : array of Double;
+ const KnownX : array of Double) : Double;
+function StandardErrorY16(const KnownY; const KnownX;
+ NData : Integer) : Double;
+ {-Returns the standard error of the predicted y-value for each X in a linear
+ regression. The standard error is a measure of the amount of error in the
+ prediction of y for an individual X.
+ }
+
+{--------------------------------------------------------------------------}
+
+{BETADIST}
+function BetaDist(X, Alpha, Beta, A, B : Single) : Single;
+ {-Returns the cumulative beta probability density function.
+ The cumulative beta probability density function is commonly used to
+ study variation in the percentage of something across samples.
+ X is the value at which to evaluate the function, A <= X <= B.
+ Alpha is a parameter to the distribution.
+ Beta is a parameter to the distribution.
+ A is the lower bound to the interval of X.
+ B is the upper bound to the interval of X.
+ The standard beta distribution has A=0 and B=1.
+ Fractional error less than 3.0e-7.
+ }
+
+{BETAINV}
+function BetaInv(Probability, Alpha, Beta, A, B : Single) : Single;
+ {-Returns the inverse of the cumulative beta probability density function.
+ That is, if Probability = BetaDist(X,...), then
+ BetaInv(Probability,...) = X.
+ Probability is a probability (0 <= p <= 1) associated with the
+ beta distribution.
+ Alpha is a parameter to the distribution.
+ Beta is a parameter to the distribution.
+ A is the lower bound to the interval of the distribution.
+ B is the upper bound to the interval of the distribution.
+ Fractional error less than 3.0e-7.
+ }
+
+{BINOMDIST}
+function BinomDist(NumberS, Trials : Integer; ProbabilityS : Single;
+ Cumulative : Boolean) : Single;
+ {-Returns the individual term binomial distribution probability.
+ Use BinomDist in problems with a fixed number of tests or Trials,
+ when the outcomes of any trial are only success or failure,
+ when Trials are independent, and when the probability of success is
+ constant throughout the experiment.
+ NumberS is the number of successes in Trials.
+ Trials is the number of independent trials.
+ ProbabilityS is the probability of success on each trial.
+ Cumulative is a logical value that determines the form of the function.
+ If Cumulative is TRUE, then BinomDist returns the cumulative
+ distribution function, which is the probability that there are at most
+ NumberS successes; if FALSE, it returns the probability mass function,
+ which is the probability that there are NumberS successes.
+ }
+
+{CRITBINOM}
+function CritBinom(Trials : Integer; ProbabilityS, Alpha : Single) : Integer;
+ {-Returns the smallest value for which the cumulative binomial distribution
+ is greater than or equal to a criterion value.
+ Trials is the number of Bernoulli trials.
+ ProbabilityS is the probability of a success on each trial.
+ Alpha is the criterion value.
+ }
+
+{CHIDIST}
+function ChiDist(X : Single; DegreesFreedom : Integer) : Single;
+ {-Returns the one-tailed probability of the chi-squared distribution.
+ The chi-squared distribution is associated with a chi-squared test.
+ Use the chi-squared test to compare observed and expected values.
+ X is the value at which you want to evaluate the distribution.
+ DegreesFreedom is the number of degrees of freedom.
+ }
+
+{CHIINV}
+function ChiInv(Probability : Single; DegreesFreedom : Integer) : Single;
+ {-Returns the inverse of the one-tailed probability of the chi-squared
+ distribution. If Probability = ChiDist(X,...), then
+ ChiInv(Probability,...) = X.
+ Probability is a probability associated with the chi-squared
+ distribution.
+ DegreesFreedom is the number of degrees of freedom.
+ }
+
+{EXPONDIST}
+function ExponDist(X, Lambda : Single; Cumulative : Boolean) : Single;
+ {-Returns the exponential distribution. Use ExponDist to model the time
+ between events.
+ X is the value of the function.
+ Lambda is the parameter value.
+ Cumulative is a logical value that indicates which form of the
+ exponential function to provide. If Cumulative is TRUE, ExponDist
+ returns the cumulative distribution function; if FALSE, it returns
+ the probability density function.
+ }
+
+{FDIST}
+function FDist(X : Single; DegreesFreedom1, DegreesFreedom2 : Integer) : Single;
+ {-Returns the F probability distribution. You can use this function to
+ determine whether two data sets have different degrees of diversity.
+ X is the value at which to evaluate the function.
+ DegreesFreedom1 is the numerator degrees of freedom.
+ DegreesFreedom2 is the denominator degrees of freedom.
+ Fractional error less than 3.0e-7.
+ }
+
+{FINV}
+function FInv(Probability : Single;
+ DegreesFreedom1, DegreesFreedom2 : Integer) : Single;
+ {-Returns the inverse of the F probability distribution. If
+ p = FDist(X,...), then FInv(p,...) = X.
+ Probability is a probability associated with the F cumulative
+ distribution.
+ DegreesFreedom1 is the numerator degrees of freedom.
+ DegreesFreedom2 is the denominator degrees of freedom.
+ Fractional error less than 3.0e-7.
+ }
+
+{LOGNORMDIST}
+function LogNormDist(X, Mean, StandardDev : Single) : Single;
+ {-Returns the cumulative lognormal distribution of X, where ln(X) is
+ normally distributed with parameters Mean and StandardDev.
+ Use this function to analyze data that has been logarithmically
+ transformed.
+ X is the value at which to evaluate the function.
+ Mean is the mean of ln(X).
+ StandardDev is the standard deviation of ln(X).
+ }
+
+{LOGINV}
+function LogInv(Probability, Mean, StandardDev : Single) : Single;
+ {-Returns the inverse of the lognormal cumulative distribution function of
+ X, where ln(X) is normally distributed with parameters Mean and
+ StandardDev. If p = LogNormDist(X,...) then LogInv(p,...) = X.
+ Probability is a probability associated with the lognormal distribution.
+ Mean is the mean of ln(X).
+ StandardDev is the standard deviation of ln(X).
+ }
+
+{NORMDIST}
+function NormDist(X, Mean, StandardDev : Single; Cumulative : Boolean) : Single;
+ {-Returns the normal cumulative distribution for the specified Mean and
+ standard deviation.
+ X is the value for which you want the distribution.
+ Mean is the arithmetic mean of the distribution.
+ StandardDev is the standard deviation of the distribution.
+ Cumulative is a logical value that determines the form of the function.
+ If Cumulative is TRUE, NormDist returns the cumulative distribution
+ function; if FALSE, it returns the probability density function.
+ }
+
+{NORMINV}
+function NormInv(Probability, Mean, StandardDev : Single) : Single;
+ {-Returns the inverse of the normal cumulative distribution for the
+ specified mean and standard deviation.
+ Probability is a probability corresponding to the normal distribution.
+ Mean is the arithmetic mean of the distribution.
+ StandardDev is the standard deviation of the distribution.
+ }
+
+{NORMSDIST}
+function NormSDist(Z : Single) : Single;
+ {-Returns the standard normal cumulative distribution function.
+ The distribution has a mean of 0 (zero) and a standard deviation of one.
+ Z is the value for which you want the distribution.
+ }
+
+{NORMSINV}
+function NormSInv(Probability : Single) : Single;
+ {-Returns the inverse of the standard normal cumulative distribution.
+ The distribution has a mean of zero and a standard deviation of one.
+ Probability is a probability corresponding to the normal distribution.
+ }
+
+{POISSON}
+function Poisson(X : Integer; Mean : Single; Cumulative : Boolean) : Single;
+ {-Returns the Poisson distribution.
+ X is the number of events.
+ Mean is the expected numeric value.
+ Cumulative is a logical value that determines the form of the
+ probability distribution returned. If Cumulative is TRUE, Poisson
+ returns the cumulative Poisson probability that the number of random
+ events occurring will be between zero and X inclusive; if FALSE,
+ it returns the Poisson probability mass function that the number of
+ events occurring will be exactly X.
+ }
+
+{TDIST}
+function TDist(X : Single; DegreesFreedom : Integer; TwoTails : Boolean) : Single;
+ {-Returns the Student's t-distribution. The t-distribution is used in the
+ hypothesis testing of small sample data sets. Use this function in place
+ of a table of critical values for the t-distribution.
+ X is the numeric value at which to evaluate the distribution.
+ DegreesFreedom is an Integer indicating the number of degrees of freedom.
+ TwoTails if a logical value that indicates the number of distribution
+ tails to return. If FALSE, TDist returns the one-tailed distribution;
+ otherwise it returns the two-tailed distribution.
+ }
+
+{TINV}
+function TInv(Probability : Single; DegreesFreedom : Integer) : Single;
+ {-Returns the inverse of the Student's t-distribution for the specified
+ degrees of freedom.
+ Probability is the probability associated with the two-tailed Student's
+ t-distribution.
+ DegreesFreedom is the number of degrees of freedom to characterize
+ the distribution.
+ }
+
+{--------------------------------------------------------------------------}
+{undocumented functions you can call if you need}
+
+function Erfc(X : Single) : Single;
+ {-Returns the complementary error function with fractional error
+ everywhere less than 1.2e-7. X is any finite value.
+ }
+
+function GammaLn(X : Single) : Single;
+ {-Returns ln(Gamma(X)) where X > 0.0.}
+
+{--------------------------------------------------------------------------}
+
+{.$DEFINE Debug}
+{$IFDEF Debug}
+{like Largest and Smallest but using slower simpler algorithm}
+function LargestSort(const Data: array of Double; K : Integer) : Double;
+function SmallestSort(const Data: array of double; K : Integer) : Double;
+{$ENDIF}
+
+implementation
+
+procedure RaiseStatError(Code : LongInt);
+ {-Generate a statistics exception}
+var
+ E : EStStatError;
+begin
+ E := EStStatError.CreateResTP(Code, 0);
+ E.ErrorCode := Code;
+ raise E;
+end;
+
+procedure DoubleArraySort(var Data; NData : Integer);
+ {-Heapsort an array of Doubles into Ascending order}
+type
+ TDoubleArray1 = array[1..StMaxBlockSize div SizeOf(Double)] of Double;
+var
+ i : Integer;
+ T : Double;
+ DA : TDoubleArray1 absolute Data;
+
+ procedure Adjust(i, N : Integer);
+ var
+ j : Integer;
+ S : Double;
+ begin
+ {j is left child of i}
+ j := 2*i;
+ {save i'th element temporarily}
+ S := DA[i];
+ while (j <= N) do begin
+ {compare left and right child}
+ if (j < N) and (DA[j] < DA[j+1]) then
+ {j indexes larger child}
+ inc(j);
+ if (S >= DA[j]) then
+ {a position for item is found}
+ break;
+ {move the larger child up a level}
+ DA[j shr 1] := DA[j];
+ {look at left child of j}
+ j := j shl 1;
+ end;
+ {store saved item}
+ DA[j shr 1] := S;
+ end;
+
+begin
+ {transform the elements into a heap}
+ for i := (NData shr 1) downto 1 do
+ Adjust(i, NData);
+
+ {repeatedly exchange the maximum at top of heap with the last element}
+ for i := NData downto 2 do begin
+ T := DA[1];
+ DA[1] := DA[i];
+ DA[i] := T;
+ {update the heap for the remaining elements}
+ Adjust(1, i-1);
+ end;
+end;
+
+function CopyAndSort(const Data; NData : Integer;
+ var SD : PDoubleArray) : Cardinal;
+ {-Allocates heap space for an array copy, moves data, sorts, and returns size}
+var
+ Size : LongInt;
+begin
+ Size := LongInt(NData)*sizeof(Double);
+ {if (Size > MaxBlockSize) then}
+ { RaiseStatError(stscStatBadCount);}
+ Result := Size;
+ getmem(SD, Size); {raises exception if insufficient memory}
+ try
+ move(Data, SD^, Size);
+ DoubleArraySort(SD^, NData);
+ except
+ freemem(SD, Size);
+ raise;
+ end;
+end;
+
+function AveDev(const Data: array of Double) : Double;
+begin
+ Result := AveDev16(Data, High(Data)+1);
+end;
+
+function Mean(const Data; NData : Integer) : Extended;
+ {-Computes the mean of an array of Doubles}
+var
+ i : Integer;
+ s : Extended;
+begin
+ s := 0.0;
+ for I := 0 to NData-1 do
+ s := s+TDoubleArray(Data)[I];
+ Result := s/NData;
+end;
+
+function AveDev16(const Data; NData : Integer) : Double;
+var
+ i : Integer;
+ m, s : Extended;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ {compute sum of absolute deviations}
+ m := Mean(Data, NData);
+ s := 0.0;
+ for i := 0 to NData-1 do
+ s := s+abs(TDoubleArray(Data)[i]-m);
+
+ Result := s/NData;
+end;
+
+function Confidence(Alpha, StandardDev : Double; Size : LongInt) : Double;
+begin
+ if (StandardDev <= 0) or (Size < 1) then
+ RaiseStatError(stscStatBadParam);
+ Result := NormSInv(1.0-Alpha/2.0)*StandardDev/sqrt(Size);
+end;
+
+function Correlation(const Data1, Data2 : array of Double) : Double;
+begin
+ if (High(Data1) <> High(Data2)) then
+ RaiseStatError(stscStatBadCount);
+ Result := Correlation16(Data1, Data2, High(Data1)+1);
+end;
+
+function Correlation16(const Data1, Data2; NData : Integer) : Double;
+var
+ sx, sy, xmean, ymean, sxx, sxy, syy, x, y : Extended;
+ i : Integer;
+begin
+ if (NData <= 1) then
+ RaiseStatError(stscStatBadCount);
+
+ {compute basic sums}
+ sx := 0.0;
+ sy := 0.0;
+ sxx := 0.0;
+ sxy := 0.0;
+ syy := 0.0;
+ for i := 0 to NData-1 do begin
+ x := TDoubleArray(Data1)[i];
+ y := TDoubleArray(Data2)[i];
+ sx := sx+x;
+ sy := sy+y;
+ sxx := sxx+x*x;
+ syy := syy+y*y;
+ sxy := sxy+x*y;
+ end;
+ xmean := sx/NData;
+ ymean := sy/NData;
+ sxx := sxx-NData*xmean*xmean;
+ syy := syy-NData*ymean*ymean;
+ sxy := sxy-NData*xmean*ymean;
+
+ Result := sxy/sqrt(sxx*syy);
+end;
+
+function Covariance(const Data1, Data2 : array of Double) : Double;
+begin
+ if (High(Data1) <> High(Data2)) then
+ RaiseStatError(stscStatBadCount);
+ Result := Covariance16(Data1, Data2, High(Data1)+1);
+end;
+
+function Covariance16(const Data1, Data2; NData : Integer) : Double;
+var
+ sx, sy, xmean, ymean, sxy, x, y : Extended;
+ i : Integer;
+begin
+ if (NData <= 1) then
+ RaiseStatError(stscStatBadCount);
+
+ {compute basic sums}
+ sx := 0.0;
+ sy := 0.0;
+ sxy := 0.0;
+ for i := 0 to NData-1 do begin
+ x := TDoubleArray(Data1)[i];
+ y := TDoubleArray(Data2)[i];
+ sx := sx+x;
+ sy := sy+y;
+ sxy := sxy+x*y;
+ end;
+ xmean := sx/NData;
+ ymean := sy/NData;
+ sxy := sxy-NData*xmean*ymean;
+
+ Result := sxy/NData;
+end;
+
+function DevSq(const Data: array of Double) : Double;
+begin
+ Result := DevSq16(Data, High(Data)+1);
+end;
+
+function DevSq16(const Data; NData : Integer) : Double;
+var
+ i : Integer;
+ sx, sxx, x : Extended;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ sx := 0.0;
+ sxx := 0.0;
+ for i := 0 to NData-1 do begin
+ x := TDoubleArray(Data)[i];
+ sx := sx+x;
+ sxx := sxx+x*x;
+ end;
+ Result := sxx-sqr(sx)/NData;
+end;
+
+procedure Frequency(const Data: array of Double; const Bins: array of Double;
+ var Counts: array of LongInt);
+begin
+ if (High(Counts) <= High(Bins)) then
+ RaiseStatError(stscStatBadCount);
+ Frequency16(Data, High(Data)+1, Bins, High(Bins)+1, Counts);
+end;
+
+procedure Frequency16(const Data; NData : Integer; const Bins; NBins : Integer;
+ var Counts);
+var
+ b, i : Integer;
+ Size : Cardinal;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) or (NBins <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {initialize all counts to zero}
+ fillchar(Counts, (NBins+1)*sizeof(Integer), 0);
+
+ {scan all data elements, putting into correct bin}
+ b := 0;
+ i := 0;
+ while (i < NData) do begin
+ if (SD^[i] <= TDoubleArray(Bins)[b]) then begin
+ {current data element falls into this bin}
+ inc(TIntArray(Counts)[b]);
+ inc(i);
+ end else begin
+ {move to next bin that collects data}
+ repeat
+ inc(b);
+ until (b = NBins) or (TDoubleArray(Bins)[b] > SD^[i]);
+ if (b = NBins) then begin
+ {add remaining elements to the catchall bin}
+ inc(TIntArray(Counts)[b], NData-i);
+ i := NData;
+ end;
+ end;
+ end;
+
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+function GeometricMean(const Data: array of Double) : Double;
+begin
+ Result := GeometricMean16(Data, High(Data)+1);
+end;
+
+function GeometricMean16(const Data; NData : Integer) : Double;
+var
+ i : Integer;
+ s, t : Extended;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ s := 1.0;
+ for i := 0 to NData-1 do begin
+ t := TDoubleArray(Data)[i];
+ if (t <= 0.0) then
+ RaiseStatError(stscStatBadData);
+ s := s*t;
+ end;
+
+ Result := Power(s, 1.0/NData);
+end;
+
+function HarmonicMean(const Data: array of Double) : Double;
+begin
+ Result := HarmonicMean16(Data, High(Data)+1);
+end;
+
+function HarmonicMean16(const Data; NData : Integer) : Double;
+var
+ i : Integer;
+ s, t : Extended;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ s := 0.0;
+ for i := 0 to NData-1 do begin
+ t := TDoubleArray(Data)[i];
+ if (t = 0.0) then
+ RaiseStatError(stscStatBadData);
+ s := s+(1.0/t);
+ end;
+ Result := NData/s;
+end;
+
+function Largest(const Data: array of Double; K : Integer) : Double;
+begin
+ Result := Largest16(Data, High(Data)+1, K);
+end;
+
+function Largest16(const Data; NData : Integer; K : Integer) : Double;
+var
+ b, t, i, j : integer;
+ Size : LongInt;
+ temp, pval : Double;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+ if (K <= 0) or (K > NData) then
+ RaiseStatError(stscStatBadParam);
+
+ Size := LongInt(NData)*sizeof(Double);
+ {if (Size > MaxBlockSize) then}
+ { RaiseStatError(stscStatBadCount);}
+ getmem(SD, Size); {raises exception if insufficient memory}
+ try
+ move(Data, SD^, Size);
+
+ {make K 0-based}
+ dec(K);
+
+ {use quicksort-like selection}
+ b := 0;
+ t := NData-1;
+ while (t > b) do begin
+ {use random pivot in case of already-sorted data}
+ pval := SD^[b+random(t-b+1)];
+ i := b;
+ j := t;
+ repeat
+ while (SD^[i] > pval) do
+ inc(i);
+ while (pval > SD^[j]) do
+ dec(j);
+ if (i <= j) then begin
+ temp := SD^[i];
+ SD^[i] := SD^[j];
+ SD^[j] := temp;
+ inc(i);
+ dec(j);
+ end;
+ until (i > j);
+ if (j < K) then
+ b := i;
+ if (K < i) then
+ t := j;
+ end;
+ Result := SD^[K];
+
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+{debug version of Largest: slower but simpler}
+{$IFDEF Debug}
+function LargestSort(const Data: array of Double; K : Integer) : Double;
+var
+ Size : Cardinal;
+ NData : Integer;
+ SD : PDoubleArray;
+begin
+ NData := High(Data)+1;
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+ if (K <= 0) or (K > NData) then
+ RaiseStatError(stscStatBadParam);
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {K=1 returns largest value, K=NData returns smallest}
+ Result := SD^[NData-K];
+ finally
+ freemem(SD, Size);
+ end;
+end;
+{$ENDIF}
+
+function Median(const Data: array of Double) : Double;
+begin
+ Result := Median16(Data, High(Data)+1);
+end;
+
+function Median16(const Data; NData : Integer) : Double;
+var
+ m : Integer;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ m := NData shr 1;
+ if odd(NData) then
+ Result := Largest16(Data, NData, m+1)
+ else
+ Result := (Largest16(Data, NData, m+1)+Largest16(Data, NData, m))/2.0;
+end;
+
+function Mode(const Data: array of Double) : Double;
+begin
+ Result := Mode16(Data, High(Data)+1);
+end;
+
+function Mode16(const Data; NData : Integer) : Double;
+var
+ maxf, i, f : Integer;
+ Size : Cardinal;
+ last, max : Double;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {find the value with highest frequency}
+ last := SD^[0];
+ max := last;
+ f := 1;
+ maxf := f;
+
+ for i := 1 to NData-1 do begin
+ if SD^[i] = last then
+ {keep count of identical values}
+ inc(f)
+ else begin
+ {start a new series}
+ if f > maxf then begin
+ max := last;
+ maxf := f;
+ end;
+ last := SD^[i];
+ f := 1;
+ end;
+ end;
+
+ {test last group}
+ if f > maxf then
+ max := last;
+
+ Result := max;
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+function Percentile(const Data: array of Double; K : Double) : Double;
+begin
+ Result := Percentile16(Data, High(Data)+1, K);
+end;
+
+function Percentile16(const Data; NData : Integer; K : Double) : Double;
+const
+ eps = 1.0e-10;
+var
+ ibin : Integer;
+ Size : Cardinal;
+ rbin, l, h : Double;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+ if (K < 0.0) or (K > 1.0) then
+ RaiseStatError(stscStatBadParam);
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {find nearest bins}
+ rbin := K*(NData-1);
+ ibin := Trunc(rbin);
+ if Frac(rbin) < eps then
+ {very close to array index below}
+ Result := SD^[ibin]
+ else if (Int(rbin)+1.0-rbin) < eps then
+ {very close to array index above}
+ Result := SD^[ibin+1]
+ else begin
+ {need to interpolate between two bins}
+ l := SD^[ibin];
+ h := SD^[ibin+1];
+ Result := l+(h-l)*(K*(NData-1)-ibin);
+ end;
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+function PercentRank(const Data: array of Double; X : Double) : Double;
+begin
+ Result := PercentRank16(Data, High(Data)+1, X);
+end;
+
+function PercentRank16(const Data; NData : Integer; X : Double) : Double;
+var
+ b, t, m : Integer;
+ Size : Cardinal;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {test end conditions}
+ if (X < SD^[0]) or (X > SD^[NData-1]) then
+ RaiseStatError(stscStatBadParam);
+
+ {find nearby bins using binary search}
+ b := 0;
+ t := NData-1;
+ while (t-b) > 1 do begin
+ m := (b+t) shr 1;
+ if (X >= SD^[m]) then
+ {search upper half}
+ b := m
+ else
+ {search lower half}
+ t := m;
+ end;
+
+ {now X is known to be between b (inclusive) and b+1}
+ {handle duplicate elements below b}
+ while (b > 0) and (SD^[b-1] = X) do
+ dec(b);
+
+ if (SD^[b] = X) then
+ {an exact match}
+ Result := b/(NData-1)
+ else
+ {interpolate}
+ Result := (b+(X-SD^[b])/(SD^[b+1]-SD^[b]))/(NData-1);
+
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+const
+ sqrt2pi = 2.5066282746310005; {sqrt(2*pi)}
+
+function GammaLn(X : Single) : Single;
+ {-Returns ln(Gamma(X)) where X > 0}
+const
+ cof : array[0..5] of Double = (
+ 76.18009172947146, -86.50532032941677, 24.01409824083091,
+ -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5);
+var
+ y, tmp, ser : Double;
+ j : Integer;
+begin
+ if (X <= 0) then
+ RaiseStatError(stscStatBadParam);
+
+ y := X;
+ tmp := X+5.5;
+ tmp := tmp-(X+0.5)*ln(tmp);
+ ser := 1.000000000190015;
+ for j := low(cof) to high(cof) do begin
+ y := y+1.0;
+ ser := ser+cof[j]/y;
+ end;
+ Result := -tmp+ln(sqrt2pi*ser/X);
+end;
+
+const
+ MFactLnA = 65;
+var
+ FactLnA : array[2..MFactLna] of Single; {lookup table of FactLn values}
+
+function FactLn(N : Integer) : Single;
+ {-Returns ln(N!) for N >= 0}
+begin
+ if (N <= 1) then
+ Result := 0.0
+ else if (N <= MFactLnA) then
+ {use lookup table}
+ Result := FactLnA[N]
+ else
+ {compute each time}
+ Result := GammaLn(N+1.0);
+end;
+
+const
+ MFactA = 33;
+var
+ FactA : array[2..MFactA] of Double; {lookup table of factorial values}
+
+function Factorial(N : Integer) : Extended;
+begin
+ if (N < 0) then
+ RaiseStatError(stscStatBadParam);
+
+ if (N <= 1) then
+ Result := 1.0
+ else if (N <= MFactA) then
+ {use lookup table}
+ Result := FactA[N]
+ else
+ {bigger than lookup table allows. may overflow!}
+ Result := exp(GammaLn(N+1.0))
+end;
+
+function Permutations(Number, NumberChosen : Integer) : Extended;
+begin
+ if (Number < 0) or (NumberChosen < 0) or (Number < NumberChosen) then
+ RaiseStatError(stscStatBadParam);
+ {the 0.5 and Int function clean up roundoff error for smaller N and K}
+ Result := Int(0.5+exp(FactLn(Number)-FactLn(Number-NumberChosen)));
+end;
+
+function Combinations(Number, NumberChosen : Integer) : Extended;
+begin
+ if (Number < 0) or (NumberChosen < 0) or (Number < NumberChosen) then
+ RaiseStatError(stscStatBadParam);
+ {the 0.5 and Int function clean up roundoff error for smaller N and K}
+ Result := Int(0.5+exp(FactLn(Number)-FactLn(NumberChosen)
+ -FactLn(Number-NumberChosen)));
+end;
+
+function Rank(Number: Double; const Data: array of Double;
+ Ascending: Boolean) : Integer;
+begin
+ Result := Rank16(Number, Data, High(Data)+1, Ascending);
+end;
+
+function Rank16(Number: Double; const Data; NData : Integer;
+ Ascending: Boolean) : Integer;
+var
+ r : Integer;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+
+ {Data not known to be sorted, so must search linearly}
+ if (Ascending) then begin
+ for r := 0 to NData-1 do
+ if TDoubleArray(Data)[r] = Number then begin
+ Result := r+1;
+ exit;
+ end;
+ end else begin
+ for r := NData-1 downto 0 do
+ if TDoubleArray(Data)[r] = Number then begin
+ Result := NData-r;
+ exit;
+ end;
+ end;
+ Result := 0;
+end;
+
+function Smallest(const Data: array of Double; K : Integer) : Double;
+begin
+ Result := Smallest16(Data, High(Data)+1, K);
+end;
+
+function Smallest16(const Data; NData : Integer; K : Integer) : Double;
+var
+ b, t, i, j : integer;
+ Size : LongInt;
+ temp, pval : Double;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+ if (K <= 0) or (K > NData) then
+ RaiseStatError(stscStatBadParam);
+
+ Size := LongInt(NData)*sizeof(Double);
+ {if (Size > MaxBlockSize) then}
+ { RaiseStatError(stscStatBadCount);}
+ getmem(SD, Size); {raises exception if insufficient memory}
+ try
+ move(Data, SD^, Size);
+
+ {make K 0-based}
+ dec(K);
+ {use quicksort-like selection}
+ b := 0;
+ t := NData-1;
+ while (t > b) do begin
+ {use random pivot in case of already-sorted data}
+ pval := SD^[b+random(t-b+1)];
+ i := b;
+ j := t;
+ repeat
+ while (SD^[i] < pval) do
+ inc(i);
+ while (pval < SD^[j]) do
+ dec(j);
+ if (i <= j) then begin
+ temp := SD^[i];
+ SD^[i] := SD^[j];
+ SD^[j] := temp;
+ inc(i);
+ dec(j);
+ end;
+ until (i > j);
+
+ if (j < K) then
+ b := i;
+ if (K < i) then
+ t := j;
+ end;
+ Result := SD^[K];
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+{debug version of Smallest: slower but simpler}
+{$IFDEF Debug}
+function SmallestSort(const Data: array of double; K : Integer) : Double;
+var
+ Size : Cardinal;
+ NData : Integer;
+ SD : PDoubleArray;
+begin
+ NData := High(Data)+1;
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+ if (K <= 0) or (K > NData) then
+ RaiseStatError(stscStatBadParam);
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {K=1 returns smallest value, K=NData returns largest}
+ Result := SD^[K-1];
+ finally
+ freemem(SD, Size);
+ end;
+end;
+{$ENDIF}
+
+function TrimMean(const Data: array of Double; Percent : Double) : Double;
+begin
+ Result := TrimMean16(Data, High(Data)+1, Percent);
+end;
+
+function TrimMean16(const Data; NData : Integer; Percent : Double) : Double;
+var
+ ntrim : Integer;
+ Size : Cardinal;
+ SD : PDoubleArray;
+begin
+ if (NData <= 0) then
+ RaiseStatError(stscStatBadCount);
+ if (Percent < 0.0) or (Percent > 1.0) then
+ RaiseStatError(stscStatBadParam);
+
+ {compute total number of trimmed points, rounding down to an even number}
+ ntrim := trunc(Percent*NData);
+ if odd(ntrim) then
+ dec(ntrim);
+
+ {take the easy way out when possible}
+ if (ntrim = 0) then begin
+ Result := Mean(Data, NData);
+ exit;
+ end;
+
+ {copy and sort array}
+ Size := CopyAndSort(Data, NData, SD);
+ try
+ {return Mean of remaining data points}
+ Result := Mean(SD^[ntrim shr 1], NData-ntrim);
+ finally
+ freemem(SD, Size);
+ end;
+end;
+
+{------------------------------------------------------------------------}
+
+procedure LinEst(const KnownY: array of Double;
+ const KnownX: array of Double; var LF : TStLinEst; ErrorStats : Boolean);
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ LinEst16(KnownY, KnownX, High(KnownY)+1, LF, ErrorStats);
+end;
+
+procedure LinEst16(const KnownY; const KnownX; NData : Integer;
+ var LF : TStLinEst; ErrorStats : Boolean);
+var
+ i : Integer;
+ sx, sy, xmean, ymean, sxx, sxy, syy, x, y : Extended;
+begin
+ if (NData <= 2) then
+ RaiseStatError(stscStatBadCount);
+
+ {compute basic sums}
+ sx := 0.0;
+ sy := 0.0;
+ sxx := 0.0;
+ sxy := 0.0;
+ syy := 0.0;
+ for i := 0 to NData-1 do begin
+ x := TDoubleArray(KnownX)[i];
+ y := TDoubleArray(KnownY)[i];
+ sx := sx+x;
+ sy := sy+y;
+ sxx := sxx+x*x;
+ syy := syy+y*y;
+ sxy := sxy+x*y;
+ end;
+ xmean := sx/NData;
+ ymean := sy/NData;
+ sxx := sxx-NData*xmean*xmean;
+ syy := syy-NData*ymean*ymean;
+ sxy := sxy-NData*xmean*ymean;
+
+ {check for zero variance}
+ if (sxx <= 0.0) or (syy <= 0.0) then
+ RaiseStatError(stscStatBadData);
+
+ {initialize returned parameters}
+ fillchar(LF, sizeof(LF), 0);
+
+ {regression coefficients}
+ LF.B1 := sxy/sxx;
+ LF.B0 := ymean-LF.B1*xmean;
+
+ {error statistics}
+ if (ErrorStats) then begin
+ LF.ssr := LF.B1*sxy;
+ LF.sse := syy-LF.ssr;
+ LF.R2 := LF.ssr/syy;
+ LF.df := NData-2;
+ LF.sigma := sqrt(LF.sse/LF.df);
+ if LF.sse = 0.0 then
+ {pick an arbitrarily large number for perfect fit}
+ LF.F0 := 1.7e+308
+ else
+ LF.F0 := (LF.ssr*LF.df)/LF.sse;
+ LF.seB1 := LF.sigma/sqrt(sxx);
+ LF.seB0 := LF.sigma*sqrt((1.0/NData)+(xmean*xmean/sxx));
+ end;
+end;
+
+procedure LogEst(const KnownY: array of Double;
+ const KnownX: array of Double; var LF : TStLinEst; ErrorStats : Boolean);
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ LogEst16(KnownY, KnownX, High(KnownY)+1, LF, ErrorStats);
+end;
+
+procedure LogEst16(const KnownY; const KnownX; NData : Integer;
+ var LF : TStLinEst; ErrorStats : Boolean);
+var
+ i : Integer;
+ Size : LongInt;
+ lny : PDoubleArray;
+begin
+ if (NData <= 2) then
+ RaiseStatError(stscStatBadCount);
+
+ {allocate array for the log-transformed data}
+ Size := LongInt(NData)*sizeof(Double);
+ {f (Size > MaxBlockSize) then}
+ { RaiseStatError(stscStatBadCount);}
+ getmem(lny, Size);
+ try
+ {initialize transformed data}
+ for i := 0 to NData-1 do
+ lny^[i] := ln(TDoubleArray(KnownY)[i]);
+
+ {fit transformed data}
+ LinEst16(lny^, KnownX, NData, LF, ErrorStats);
+
+ {return values for B0 and B1 in exponential model y=B0*B1^x}
+ LF.B0 := exp(LF.B0);
+ LF.B1 := exp(LF.B1);
+ {leave other values in LF in log form}
+ finally
+ freemem(lny, Size);
+ end;
+end;
+
+function Forecast(X : Double; const KnownY: array of Double;
+ const KnownX: array of Double) : Double;
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ Result := Forecast16(X, KnownY, KnownX, High(KnownY)+1);
+end;
+
+function Forecast16(X : Double; const KnownY; const KnownX;
+ NData : Integer) : Double;
+var
+ LF : TStLinEst;
+begin
+ LinEst16(KnownY, KnownX, NData, LF, false);
+ Result := LF.B0+LF.B1*X;
+end;
+
+function ForecastExponential(X : Double; const KnownY: array of Double;
+ const KnownX: array of Double) : Double;
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ Result := ForecastExponential16(X, KnownY, KnownX, High(KnownY)+1);
+end;
+
+function ForecastExponential16(X : Double; const KnownY; const KnownX;
+ NData : Integer) : Double;
+var
+ LF : TStLinEst;
+begin
+ LogEst16(KnownY, KnownX, NData, LF, false);
+ Result := LF.B0*Power(LF.B1, X);
+end;
+
+function Intercept(const KnownY: array of Double;
+ const KnownX: array of Double) : Double;
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ Result := Intercept16(KnownY, KnownX, High(KnownY)+1);
+end;
+
+function Intercept16(const KnownY; const KnownX; NData : Integer) : Double;
+var
+ LF : TStLinEst;
+begin
+ LinEst16(KnownY, KnownX, NData, LF, false);
+ Result := LF.B0;
+end;
+
+function RSquared(const KnownY: array of Double;
+ const KnownX: array of Double) : Double;
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ Result := RSquared16(KnownY, KnownX, High(KnownY)+1);
+end;
+
+function RSquared16(const KnownY; const KnownX; NData : Integer) : Double;
+var
+ LF : TStLinEst;
+begin
+ LinEst16(KnownY, KnownX, NData, LF, true);
+ Result := LF.R2;
+end;
+
+function Slope(const KnownY: array of Double;
+ const KnownX: array of Double) : Double;
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ Result := Slope16(KnownY, KnownX, High(KnownY)+1);
+end;
+
+function Slope16(const KnownY; const KnownX; NData : Integer) : Double;
+var
+ LF : TStLinEst;
+begin
+ LinEst16(KnownY, KnownX, NData, LF, false);
+ Result := LF.B1;
+end;
+
+function StandardErrorY(const KnownY: array of Double;
+ const KnownX: array of Double) : Double;
+begin
+ if (High(KnownY) <> High(KnownX)) then
+ RaiseStatError(stscStatBadCount);
+ Result := StandardErrorY16(KnownY, KnownX, High(KnownY)+1);
+end;
+
+function StandardErrorY16(const KnownY; const KnownX;
+ NData : Integer) : Double;
+var
+ LF : TStLinEst;
+begin
+ LinEst16(KnownY, KnownX, NData, LF, true);
+ Result := LF.sigma;
+end;
+
+{------------------------------------------------------------------------}
+
+function BetaCf(a, b, x : Single) : Single;
+ {-Evaluates continued fraction for incomplete beta function}
+const
+ MAXIT = 100;
+ EPS = 3.0E-7;
+ FPMIN = 1.0E-30;
+var
+ m, m2 : Integer;
+ aa, c, d, del, h, qab, qam, qap : Double;
+begin
+ qab := a+b;
+ qap := a+1.0;
+ qam := a-1.0;
+ c := 1.0;
+ d := 1.0-qab*x/qap;
+ if (abs(d) < FPMIN) then
+ d := FPMIN;
+ d := 1.0/d;
+ h := d;
+ for m := 1 to MAXIT do begin
+ m2 := 2*m;
+ aa := m*(b-m)*x/((qam+m2)*(a+m2));
+ d := 1.0+aa*d;
+ if (abs(d) < FPMIN) then
+ d := FPMIN;
+ c := 1.0+aa/c;
+ if (abs(c) < FPMIN) then
+ c := FPMIN;
+ d := 1.0/d;
+ h := h*d*c;
+ aa := -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
+
+ d := 1.0+aa*d;
+ if (abs(d) < FPMIN) then
+ d := FPMIN;
+ c := 1.0+aa/c;
+ if (abs(c) < FPMIN) then
+ c := FPMIN;
+ d := 1.0/d;
+ del := d*c;
+ h := h*del;
+
+ {check for convergence}
+ if (abs(del-1.0) < EPS) then
+ break;
+ if m = MAXIT then
+ RaiseStatError(stscStatNoConverge);
+ end;
+ Result := h;
+end;
+
+function BetaDist(X, Alpha, Beta, A, B : Single) : Single;
+var
+ bt : Double;
+begin
+ if (X < A) or (X > B) or (A = B) or (Alpha <= 0.0) or (Beta <= 0.0) then
+ RaiseStatError(stscStatBadParam);
+
+ {normalize X}
+ X := (X-A)/(B-A);
+
+ {compute factors in front of continued fraction expansion}
+ if (X = 0.0) or (X = 1.0) then
+ bt := 0.0
+ else
+ bt := exp(GammaLn(Alpha+Beta)-GammaLn(Alpha)-GammaLn(Beta)+
+ Alpha*ln(X)+Beta*ln(1.0-X));
+
+ {use symmetry relationship to make continued fraction converge quickly}
+ if (X < (Alpha+1.0)/(Alpha+Beta+2.0)) then
+ Result := bt*BetaCf(Alpha, Beta, X)/Alpha
+ else
+ Result := 1.0-bt*BetaCf(Beta, Alpha, 1.0-X)/Beta;
+end;
+
+function Sign(a, b : Double) : Double;
+ {-Returns abs(a) if b >= 0.0, else -abs(a)}
+begin
+ if (b >= 0.0) then
+ Result := abs(a)
+ else
+ Result := -abs(a);
+end;
+
+function BetaInv(Probability, Alpha, Beta, A, B : Single) : Single;
+const
+ MAXIT = 100;
+ UNUSED = -1.11e30;
+ XACC = 3e-7;
+var
+ j : Integer;
+ ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) or
+ (A = B) or (Alpha <= 0.0) or (Beta <= 0.0) then
+ RaiseStatError(stscStatBadParam);
+
+ if (Probability = 0.0) then
+ Result := A
+ else if (Probability = 1.0) then
+ Result := B
+ else begin
+ {Ridder's method of finding the root of BetaDist = Probability}
+ fl := -Probability; {BetaDist(A, Alpha, Beta, A, B)-Probability}
+ fh := 1.0-Probability; {BetaDist(B, Alpha, Beta, A, B)-Probability}
+ xl := A;
+ xh := B;
+ ans := UNUSED;
+
+ for j := 1 to MAXIT do begin
+ xm := 0.5*(xl+xh);
+ fm := BetaDist(xm, Alpha, Beta, A, B)-Probability;
+ s := sqrt(fm*fm-fl*fh);
+ if (s = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+ if (fl >= fh) then
+ dsign := 1.0
+ else
+ dsign := -1.0;
+ xnew := xm+(xm-xl)*(dsign*fm/s);
+ if (abs(xnew-ans) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ ans := xnew;
+
+ fnew := BetaDist(ans, Alpha, Beta, A, B)-Probability;
+ if (fnew = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+
+ {keep root bracketed on next iteration}
+ if (Sign(fm, fnew) <> fm) then begin
+ xl := xm;
+ fl := fm;
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fl, fnew) <> fl) then begin
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fh, fnew) <> fh) then begin
+ xl := ans;
+ fl := fnew;
+ end else
+ {shouldn't get here}
+ RaiseStatError(stscStatNoConverge);
+
+ if (abs(xh-xl) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ end;
+ BetaInv := A; {avoid compiler warning}
+ RaiseStatError(stscStatNoConverge);
+ end;
+end;
+
+function BinomDistPmf(N, K : Integer; p : Extended) : Double;
+ {-Returns the Probability mass function of the binomial distribution}
+begin
+ Result := Combinations(N, K)*IntPower(p, K)*IntPower(1.0-p, N-K);
+end;
+
+function BinomDist(NumberS, Trials : Integer; ProbabilityS : Single;
+ Cumulative : Boolean) : Single;
+begin
+ if (Trials < 0) or (NumberS < 0) or (NumberS > Trials) or
+ (ProbabilityS < 0.0) or (ProbabilityS > 1.0) then
+ RaiseStatError(stscStatBadParam);
+
+ if (Cumulative) then
+ Result := 1.0+BinomDistPmf(Trials, NumberS, ProbabilityS)-
+ BetaDist(ProbabilityS, NumberS, Trials-NumberS+1, 0.0, 1.0)
+ else
+ Result := BinomDistPmf(Trials, NumberS, ProbabilityS);
+end;
+
+function CritBinom(Trials : Integer; ProbabilityS, Alpha : Single) : Integer;
+var
+ s : Integer;
+ B : Double;
+begin
+ if (Trials < 0) or (ProbabilityS < 0.0) or (ProbabilityS > 1.0) or
+ (Alpha < 0.0) or (Alpha > 1.0) then
+ RaiseStatError(stscStatBadParam);
+
+ B := 0.0;
+ for s := 0 to Trials do begin
+ B := B+BinomDistPmf(Trials, s, ProbabilityS);
+ if (B >= Alpha) then begin
+ Result := s;
+ exit;
+ end;
+ end;
+ {in case of roundoff problems}
+ Result := Trials;
+end;
+
+function GammSer(a, x : Single) : Single;
+ {-Returns the series computation for GammP}
+const
+ MAXIT = 100;
+ EPS = 3.0E-7;
+var
+ N : Integer;
+ sum, del, ap : Double;
+begin
+ Result := 0.0;
+ if (x > 0.0) then begin
+ {x shouldn't be < 0.0, tested by caller}
+ ap := a;
+ sum := 1.0/a;
+ del := sum;
+ for N := 1 to MAXIT do begin
+ ap := ap+1;
+ del := del*x/ap;
+ sum := sum+del;
+ if (abs(del) < abs(sum)*EPS) then begin
+ Result := sum*exp(-X+a*ln(X)-GammaLn(a));
+ exit;
+ end;
+ end;
+ RaiseStatError(stscStatNoConverge);
+ end;
+end;
+
+function GammCf(a, x : Single) : Single;
+ {-Returns the continued fraction computation for GammP}
+const
+ MAXIT = 100;
+ EPS = 3.0e-7;
+ FPMIN = 1.0e-30;
+var
+ i : Integer;
+ an, b, c, d, del, h : Double;
+begin
+ b := x+1.0-a;
+ c := 1.0/FPMIN;
+ d := 1.0/b;
+ h := d;
+ for i := 1 to MAXIT do begin
+ an := -i*(i-a);
+ b := b+2.0;
+ d := an*d+b;
+ if (abs(d) < FPMIN) then
+ d := FPMIN;
+ c := b+an/c;
+ if (abs(c) < FPMIN) then
+ c := FPMIN;
+ d := 1.0/d;
+ del := d*c;
+ h := h*del;
+ if (abs(del-1.0) < EPS) then
+ break;
+ if i = MAXIT then
+ RaiseStatError(stscStatNoConverge);
+ end;
+ Result := h*exp(-x+a*ln(x)-GammaLn(a));
+end;
+
+function GammP(a, x : Single) : Single;
+ {-Returns the incomplete gamma function P(a, x)}
+begin
+ if (x < 0.0) or (a <= 0.0) then
+ RaiseStatError(stscStatBadParam);
+ if (x < a+1.0) then
+ {use the series representation}
+ Result := GammSer(a, x)
+ else
+ {use the continued fraction representation}
+ Result := 1.0-GammCf(a, x);
+end;
+
+function ChiDist(X : Single; DegreesFreedom : Integer) : Single;
+begin
+ if (DegreesFreedom < 1) or (X < 0.0) then
+ RaiseStatError(stscStatBadParam);
+ Result := 1.0-GammP(DegreesFreedom/2.0, X/2.0);
+end;
+
+function ChiInv(Probability : Single; DegreesFreedom : Integer) : Single;
+const
+ MAXIT = 100;
+ UNUSED = -1.11e30;
+ XACC = 3e-7;
+ FACTOR = 1.6;
+var
+ j : Integer;
+ ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) or (DegreesFreedom < 1) then
+ RaiseStatError(stscStatBadParam);
+
+ if (Probability = 0.0) then
+ Result := 0.0
+ else begin
+ {bracket the interval}
+ xl := 0.0;
+ xh := 100.0;
+ fl := ChiDist(xl, DegreesFreedom)-Probability;
+ fh := ChiDist(xh, DegreesFreedom)-Probability;
+ for j := 1 to MAXIT do begin
+ if (fl*fh < 0.0) then
+ {bracketed the root}
+ break;
+ if (abs(fl) < abs(fh)) then begin
+ xl := xl+FACTOR*(xl-xh);
+ fl := ChiDist(xl, DegreesFreedom)-Probability;
+ end else begin
+ xh := xh+FACTOR*(xh-xl);
+ fh := ChiDist(xh, DegreesFreedom)-Probability;
+ end;
+ end;
+ if (fl*fh >= 0.0) then
+ {couldn't bracket it, means Probability too close to 1.0}
+ RaiseStatError(stscStatNoConverge);
+
+ {Ridder's method of finding the root of ChiDist = Probability}
+ ans := UNUSED;
+
+ for j := 1 to MAXIT do begin
+ xm := 0.5*(xl+xh);
+ fm := ChiDist(xm, DegreesFreedom)-Probability;
+ s := sqrt(fm*fm-fl*fh);
+ if (s = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+ if (fl >= fh) then
+ dsign := 1.0
+ else
+ dsign := -1.0;
+ xnew := xm+(xm-xl)*(dsign*fm/s);
+ if (abs(xnew-ans) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ ans := xnew;
+
+ fnew := ChiDist(ans, DegreesFreedom)-Probability;
+ if (fnew = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+
+ {keep root bracketed on next iteration}
+ if (Sign(fm, fnew) <> fm) then begin
+ xl := xm;
+ fl := fm;
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fl, fnew) <> fl) then begin
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fh, fnew) <> fh) then begin
+ xl := ans;
+ fl := fnew;
+ end else
+ {shouldn't get here}
+ RaiseStatError(stscStatNoConverge);
+
+ if (abs(xh-xl) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ end;
+ Result := 0.0; {avoid compiler warning}
+ RaiseStatError(stscStatNoConverge);
+ end;
+end;
+
+function ExponDist(X, Lambda : Single; Cumulative : Boolean) : Single;
+begin
+ if (X < 0.0) or (Lambda <= 0.0) then
+ RaiseStatError(stscStatBadParam);
+
+ if (Cumulative) then
+ Result := 1.0-exp(-Lambda*X)
+ else
+ Result := Lambda*exp(-Lambda*X);
+end;
+
+function FDist(X : Single;
+ DegreesFreedom1, DegreesFreedom2 : Integer) : Single;
+begin
+ if (X < 0) or (DegreesFreedom1 < 1) or (DegreesFreedom2 < 1) then
+ RaiseStatError(stscStatBadParam);
+
+ Result := BetaDist(DegreesFreedom2/(DegreesFreedom2+DegreesFreedom1*X),
+ DegreesFreedom2/2.0, DegreesFreedom1/2.0, 0.0, 1.0);
+end;
+
+function FInv(Probability : Single;
+ DegreesFreedom1, DegreesFreedom2 : Integer) : Single;
+const
+ MAXIT = 100;
+ UNUSED = -1.11e30;
+ XACC = 3e-7;
+ FACTOR = 1.6;
+var
+ j : Integer;
+ ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) or
+ (DegreesFreedom1 < 1) or (DegreesFreedom2 < 1) then
+ RaiseStatError(stscStatBadParam);
+
+ if (Probability = 1.0) then
+ Result := 0.0
+ else begin
+ {bracket the interval}
+ xl := 0.0;
+ xh := 100.0;
+ fl := FDist(xl, DegreesFreedom1, DegreesFreedom2)-Probability;
+ fh := FDist(xh, DegreesFreedom1, DegreesFreedom2)-Probability;
+ for j := 1 to MAXIT do begin
+ if (fl*fh < 0.0) then
+ {bracketed the root}
+ break;
+ if (abs(fl) < abs(fh)) then begin
+ xl := xl+FACTOR*(xl-xh);
+ fl := FDist(xl, DegreesFreedom1, DegreesFreedom2)-Probability;
+ end else begin
+ xh := xh+FACTOR*(xh-xl);
+ fh := FDist(xh, DegreesFreedom1, DegreesFreedom2)-Probability;
+ end;
+ end;
+ if (fl*fh >= 0.0) then
+ {couldn't bracket it, means Probability too close to 0.0}
+ RaiseStatError(stscStatNoConverge);
+
+ {Ridder's method of finding the root of FDist = Probability}
+ ans := UNUSED;
+
+ for j := 1 to MAXIT do begin
+ xm := 0.5*(xl+xh);
+ fm := FDist(xm, DegreesFreedom1, DegreesFreedom2)-Probability;
+ s := sqrt(fm*fm-fl*fh);
+ if (s = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+ if (fl >= fh) then
+ dsign := 1.0
+ else
+ dsign := -1.0;
+ xnew := xm+(xm-xl)*(dsign*fm/s);
+ if (abs(xnew-ans) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ ans := xnew;
+
+ fnew := FDist(ans, DegreesFreedom1, DegreesFreedom2)-Probability;
+ if (fnew = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+
+ {keep root bracketed on next iteration}
+ if (Sign(fm, fnew) <> fm) then begin
+ xl := xm;
+ fl := fm;
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fl, fnew) <> fl) then begin
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fh, fnew) <> fh) then begin
+ xl := ans;
+ fl := fnew;
+ end else
+ {shouldn't get here}
+ RaiseStatError(stscStatNoConverge);
+
+ if (abs(xh-xl) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ end;
+ Result := 0.0; {avoid compiler warning}
+ RaiseStatError(stscStatNoConverge);
+ end;
+end;
+
+function LogNormDist(X, Mean, StandardDev : Single) : Single;
+begin
+ if (X <= 0.0) or (StandardDev <= 0) then
+ RaiseStatError(stscStatBadParam);
+ Result := NormSDist((ln(X)-Mean)/StandardDev);
+end;
+
+function LogInv(Probability, Mean, StandardDev : Single) : Single;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) or (StandardDev <= 0) then
+ RaiseStatError(stscStatBadParam);
+ Result := exp(Mean+StandardDev*NormSInv(Probability));
+end;
+
+function NormDist(X, Mean, StandardDev : Single;
+ Cumulative : Boolean) : Single;
+var
+ Z : Extended;
+begin
+ if (StandardDev <= 0) then
+ RaiseStatError(stscStatBadParam);
+ Z := (X-Mean)/StandardDev;
+ if (Cumulative) then
+ Result := NormSDist(Z)
+ else
+ Result := exp(-Z*Z/2.0)/(StandardDev*sqrt2pi);
+end;
+
+function NormInv(Probability, Mean, StandardDev : Single) : Single;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) or (StandardDev <= 0) then
+ RaiseStatError(stscStatBadParam);
+ Result := Mean+StandardDev*NormSInv(Probability);
+end;
+
+function Erfc(X : Single) : Single;
+var
+ t, z, ans : Double;
+begin
+ z := abs(X);
+ t := 1.0/(1.0+0.5*z);
+ ans := t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
+ t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
+ t*(-0.82215223+t*0.17087277)))))))));
+ if (X >= 0.0) then
+ Result := ans
+ else
+ Result := 2.0-ans;
+end;
+
+function NormSDist(Z : Single) : Single;
+const
+ sqrt2 = 1.41421356237310;
+begin
+ Result := 1.0-0.5*Erfc(Z/sqrt2);
+end;
+
+function NormSInv(Probability : Single) : Single;
+const
+ MAXIT = 100;
+ UNUSED = -1.11e30;
+ XACC = 3e-7;
+ FACTOR = 1.6;
+var
+ j : Integer;
+ ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) then
+ RaiseStatError(stscStatBadParam);
+ Result := 0.0;
+
+ {bracket the interval}
+ xl := -2.0;
+ xh := +2.0;
+ fl := NormSDist(xl)-Probability;
+ fh := NormSDist(xh)-Probability;
+ for j := 1 to MAXIT do begin
+ if (fl*fh < 0.0) then
+ {bracketed the root}
+ break;
+ if (abs(fl) < abs(fh)) then begin
+ xl := xl+FACTOR*(xl-xh);
+ fl := NormSDist(xl)-Probability;
+ end else begin
+ xh := xh+FACTOR*(xh-xl);
+ fh := NormSDist(xh)-Probability;
+ end;
+ end;
+ if (fl*fh >= 0.0) then
+ {couldn't bracket it, means Probability too close to 0.0 or 1.0}
+ RaiseStatError(stscStatNoConverge);
+
+ {Ridder's method of finding the root of NormSDist = Probability}
+ ans := UNUSED;
+
+ for j := 1 to MAXIT do begin
+ xm := 0.5*(xl+xh);
+ fm := NormSDist(xm)-Probability;
+ s := sqrt(fm*fm-fl*fh);
+ if (s = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+ if (fl >= fh) then
+ dsign := 1.0
+ else
+ dsign := -1.0;
+ xnew := xm+(xm-xl)*(dsign*fm/s);
+ if (abs(xnew-ans) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ ans := xnew;
+
+ fnew := NormSDist(ans)-Probability;
+ if (fnew = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+
+ {keep root bracketed on next iteration}
+ if (Sign(fm, fnew) <> fm) then begin
+ xl := xm;
+ fl := fm;
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fl, fnew) <> fl) then begin
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fh, fnew) <> fh) then begin
+ xl := ans;
+ fl := fnew;
+ end else
+ {shouldn't get here}
+ RaiseStatError(stscStatNoConverge);
+
+ if (abs(xh-xl) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ end;
+ RaiseStatError(stscStatNoConverge);
+end;
+
+function Poisson(X : Integer; Mean : Single; Cumulative : Boolean) : Single;
+begin
+ if (X < 0) or (Mean <= 0.0) then
+ RaiseStatError(stscStatBadParam);
+ if (Cumulative) then
+ Result := 1.0-GammP(X+1.0, Mean)
+ else
+ Result := IntPower(Mean, X)*exp(-Mean)/Factorial(X);
+end;
+
+function TDist(X : Single; DegreesFreedom : Integer;
+ TwoTails : Boolean) : Single;
+var
+ a : Double;
+begin
+ if (DegreesFreedom < 1) then
+ RaiseStatError(stscStatBadParam);
+
+ a := BetaDist(DegreesFreedom/(DegreesFreedom+X*X), DegreesFreedom/2.0,
+ 0.5, 0.0, 1.0);
+ if TwoTails then
+ Result := a
+ else
+ Result := 0.5*a;
+end;
+
+function TInv(Probability : Single; DegreesFreedom : Integer) : Single;
+const
+ MAXIT = 100;
+ UNUSED = -1.11e30;
+ XACC = 3e-7;
+ FACTOR = 1.6;
+var
+ j : Integer;
+ ans, fh, fl, fm, fnew, s, xh, xl, xm, xnew, dsign : Double;
+begin
+ if (Probability < 0.0) or (Probability > 1.0) or (DegreesFreedom < 1) then
+ RaiseStatError(stscStatBadParam);
+
+ Result := 0.0;
+ if (Probability = 1.0) then
+ exit;
+
+ {bracket the interval}
+ xl := 0.0;
+ xh := +2.0;
+ fl := TDist(xl, DegreesFreedom, true)-Probability;
+ fh := TDist(xh, DegreesFreedom, true)-Probability;
+ for j := 1 to MAXIT do begin
+ if (fl*fh < 0.0) then
+ {bracketed the root}
+ break;
+ if (abs(fl) < abs(fh)) then begin
+ xl := xl+FACTOR*(xl-xh);
+ fl := TDist(xl, DegreesFreedom, true)-Probability;
+ end else begin
+ xh := xh+FACTOR*(xh-xl);
+ fh := TDist(xh, DegreesFreedom, true)-Probability;
+ end;
+ end;
+ if (fl*fh >= 0.0) then
+ {couldn't bracket it, means Probability too close to 1.0}
+ RaiseStatError(stscStatNoConverge);
+
+ {Ridder's method of finding the root of TDist = Probability}
+ ans := UNUSED;
+
+ for j := 1 to MAXIT do begin
+ xm := 0.5*(xl+xh);
+ fm := TDist(xm, DegreesFreedom, true)-Probability;
+ s := sqrt(fm*fm-fl*fh);
+ if (s = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+ if (fl >= fh) then
+ dsign := 1.0
+ else
+ dsign := -1.0;
+ xnew := xm+(xm-xl)*(dsign*fm/s);
+ if (abs(xnew-ans) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ ans := xnew;
+
+ fnew := TDist(ans, DegreesFreedom, true)-Probability;
+ if (fnew = 0.0) then begin
+ Result := ans;
+ exit;
+ end;
+
+ {keep root bracketed on next iteration}
+ if (Sign(fm, fnew) <> fm) then begin
+ xl := xm;
+ fl := fm;
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fl, fnew) <> fl) then begin
+ xh := ans;
+ fh := fnew;
+ end else if (Sign(fh, fnew) <> fh) then begin
+ xl := ans;
+ fl := fnew;
+ end else
+ {shouldn't get here}
+ RaiseStatError(stscStatNoConverge);
+
+ if (abs(xh-xl) <= XACC) then begin
+ Result := ans;
+ exit;
+ end;
+ end;
+ RaiseStatError(stscStatNoConverge);
+end;
+
+procedure Initialize;
+ {-Fully initializes factorial lookup tables}
+var
+ i : Integer;
+begin
+ FactA[2] := 2.0;
+ for i := 3 to MFactA do
+ FactA[i] := i*FactA[i-1];
+ for i := 2 to MFactLna do
+ FactLnA[i] := GammaLn(i+1.0);
+end;
+
+initialization
+ Initialize;
+end.