From 560fd631fa99d78690a6177414cc448b89b2270d Mon Sep 17 00:00:00 2001 From: wp_xxyyzz Date: Wed, 17 Jan 2018 08:04:35 +0000 Subject: [PATCH] systools: Add units StRandom and StStat (and demo for StRandom). git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@6141 8e941d3f-bd1b-0410-a28a-d453659cc2b4 --- .../systools/examples/random/exrandom.lpi | 84 + .../systools/examples/random/exrandom.lpr | 43 + .../systools/examples/random/exrndu.lfm | 146 + .../systools/examples/random/exrndu.pas | 523 ++++ components/systools/laz_systools.lpk | 10 +- components/systools/laz_systools.pas | 2 +- components/systools/source/design/StReg.pas | 8 +- components/systools/source/run/strandom.pas | 735 ++++++ components/systools/source/run/ststat.pas | 2338 +++++++++++++++++ 9 files changed, 3885 insertions(+), 4 deletions(-) create mode 100644 components/systools/examples/random/exrandom.lpi create mode 100644 components/systools/examples/random/exrandom.lpr create mode 100644 components/systools/examples/random/exrndu.lfm create mode 100644 components/systools/examples/random/exrndu.pas create mode 100644 components/systools/source/run/strandom.pas create mode 100644 components/systools/source/run/ststat.pas 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 @@ + + + + + + + + + + + + + <Scaled Value="True"/> + <ResourceType Value="res"/> + <UseXPManifest Value="True"/> + <XPManifest> + <DpiAware Value="True"/> + </XPManifest> + <Icon Value="0"/> + </General> + <BuildModes Count="1"> + <Item1 Name="Default" Default="True"/> + </BuildModes> + <PublishOptions> + <Version Value="2"/> + </PublishOptions> + <RunParams> + <FormatVersion Value="2"/> + <Modes Count="0"/> + </RunParams> + <RequiredPackages Count="2"> + <Item1> + <PackageName Value="laz_systools"/> + </Item1> + <Item2> + <PackageName Value="LCL"/> + </Item2> + </RequiredPackages> + <Units Count="2"> + <Unit0> + <Filename Value="exrandom.lpr"/> + <IsPartOfProject Value="True"/> + </Unit0> + <Unit1> + <Filename Value="exrndu.pas"/> + <IsPartOfProject Value="True"/> + <ComponentName Value="Form1"/> + <HasResources Value="True"/> + <ResourceBaseClass Value="Form"/> + <UnitName Value="ExRndU"/> + </Unit1> + </Units> + </ProjectOptions> + <CompilerOptions> + <Version Value="11"/> + <PathDelim Value="\"/> + <Target> + <Filename Value="exrandom"/> + </Target> + <SearchPaths> + <UnitOutputDirectory Value="lib\$(TargetCPU)-$(TargetOS)"/> + </SearchPaths> + <Linking> + <Options> + <Win32> + <GraphicApplication Value="True"/> + </Win32> + </Options> + </Linking> + </CompilerOptions> + <Debugging> + <Exceptions Count="3"> + <Item1> + <Name Value="EAbort"/> + </Item1> + <Item2> + <Name Value="ECodetoolError"/> + </Item2> + <Item3> + <Name Value="EFOpenError"/> + </Item3> + </Exceptions> + </Debugging> +</CONFIG> 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 @@ <Description Value="Lazarus port of TurboPower SysTools (Delphi version at https://sourceforge.net/projects/tpsystools/) - runtime package"/> <License Value="MPL 1.1"/> <Version Major="4" Release="4"/> - <Files Count="21"> + <Files Count="23"> <Item1> <Filename Value="source\run\stbarc.pas"/> <UnitName Value="StBarC"/> @@ -101,6 +101,14 @@ <Filename Value="source\run\stmoney.pas"/> <UnitName Value="StMoney"/> </Item21> + <Item22> + <Filename Value="source\run\strandom.pas"/> + <UnitName Value="StRandom"/> + </Item22> + <Item23> + <Filename Value="source\run\ststat.pas"/> + <UnitName Value="StStat"/> + </Item23> </Files> <RequiredPkgs Count="2"> <Item1> 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.