pascal/lazarus-ccr
1
0
Fork 0
Code Issues Releases Activity

LazStats: Add and integrte new standalone units MatrixUnit and RegressionUnit. Move some code from MathUnit to MatrixUnit or RegressionUnit.

Browse Source 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7777 8e941d3f-bd1b-0410-a28a-d453659cc2b4
This commit is contained in:
wp_xxyyzz
2020-10-14 14:28:59 +00:00
parent 9013e01332
commit af6047d214
24 changed files with 975 additions and 408 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
applications/lazstats/source/LazStats.lpi
View File

@ -116,7 +116,7 @@
<PackageName Value="LCL"/>
</Item7>
</RequiredPackages>
<Units Count="181">
<Units Count="183">
<Unit0>
<Filename Value="LazStats.lpr"/>
<IsPartOfProject Value="True"/>
@ -1533,6 +1533,16 @@
<IsPartOfProject Value="True"/>
<UnitName Value="GridProcs"/>
</Unit180>
<Unit181>
<Filename Value="units\matrixunit.pas"/>
<IsPartOfProject Value="True"/>
<UnitName Value="MatrixUnit"/>
</Unit181>
<Unit182>
<Filename Value="units\regressionunit.pas"/>
<IsPartOfProject Value="True"/>
<UnitName Value="RegressionUnit"/>
</Unit182>
</Units>
</ProjectOptions>
<CompilerOptions>

2
applications/lazstats/source/LazStats.lpr
View File

@ -8,7 +8,7 @@ uses
{$ENDIF}{$ENDIF}
Interfaces, // this includes the LCL widgetset
Forms, tachartlazaruspkg, tachartprint, lhelpcontrolpkg, Globals, LicenseUnit,
OptionsUnit, MainDM, MainUnit, GridProcs;
OptionsUnit, MainDM, MainUnit, GridProcs, MatrixUnit, RegressionUnit;
{$R LazStats.res}

6
applications/lazstats/source/forms/analysis/descriptive/descriptiveunit.pas
View File

@ -118,7 +118,7 @@ implementation
uses
Math,
Utils, MathUnit, GridProcs;
Utils, MatrixUnit, GridProcs;
{===============================================================================
@ -153,8 +153,8 @@ begin
SortOnX(FValues);
MathUnit.Calc_MaxMin(FValues, FMax, FMin);
MathUnit.Calc_MeanVarStdDev(FValues, FMean, FVariance, FStdDev);
VecMaxMin(FValues, FMax, FMin);
VecMeanVarStdDev(FValues, FMean, FVariance, FStdDev);
if FNumCases > 1 then
FStdErrorMean := sqrt(FVariance / FNumCases);

2
applications/lazstats/source/forms/analysis/descriptive/multxvsyunit.pas
View File

@ -67,7 +67,7 @@ implementation
uses
TATypes,
Math, Utils, MathUnit;
Math, Utils, MatrixUnit;
{ TMultXvsYFrm }

6
applications/lazstats/source/forms/analysis/descriptive/plotxyunit.pas
View File

@ -9,7 +9,7 @@ interface
uses
Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
StdCtrls, ExtCtrls, Buttons, ComCtrls, Grids,
MainUnit, Globals, MathUnit, BasicStatsReportAndChartFormUnit,
MainUnit, Globals, RegressionUnit, BasicStatsReportAndChartFormUnit,
ReportFrameUnit, ChartFrameUnit;
type
@ -71,7 +71,7 @@ implementation
uses
TAChartUtils, TAChartAxisUtils, TALegend, TASources, TACustomSeries, TASeries,
GridProcs, Utils;
GridProcs, Utils, MatrixUnit;
{ TPlotXYForm }
@ -145,7 +145,7 @@ begin
// Calculate regression
confBand := StrToFloat(ConfEdit.Text) / 100.0;
Calc_BivariateRegression(xValues, yValues, confBand, regressionRes);
BivariateRegression(xValues, yValues, confBand, regressionRes);
// Print the descriptive statistics to the output frame
WriteToReport(regressionRes);

16
applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas
View File

@ -13,7 +13,7 @@ uses
Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
StdCtrls, ExtCtrls, Buttons, Printers, ComCtrls, Grids,
MainUnit, Globals, DataProcs, DictionaryUnit,
MathUnit, BasicStatsReportAndChartFormUnit;
MatrixUnit, RegressionUnit, BasicStatsReportAndChartFormUnit;
type
@ -165,7 +165,7 @@ begin
// Calculate bivariate regression
confBand := StrToFloat(ConfEdit.Text) / 100;
Calc_BivariateRegression(xValues, yValues, confBand, regressionRes);
BivariateRegression(xValues, yValues, confBand, regressionRes);
// Do the resistant line analysis
ResistantLineAnalysis(xValues, yValues, xMedians, yMedians, grpSize);
@ -476,8 +476,8 @@ begin
xVector[i] := XData[i + offs];
yVector[i] := YData[i + offs];
end;
xMedians[0] := Calc_Median(xVector);
yMedians[0] := Calc_Median(yVector);
xMedians[0] := VecMedian(xVector);
yMedians[0] := VecMedian(yVector);
// Second group
SetLength(xVector, GrpSize[1]);
@ -488,8 +488,8 @@ begin
xVector[i] := XData[i + offs];
yVector[i] := YData[i + offs];
end;
xMedians[1] := Calc_Median(xVector);
yMedians[1] := Calc_Median(yVector);
xMedians[1] := VecMedian(xVector);
yMedians[1] := VecMedian(yVector);
// Third group
SetLength(xVector, GrpSize[2]);
@ -500,8 +500,8 @@ begin
xVector[i] := XData[i + offs];
yVector[i] := YData[i + offs];
end;
xMedians[2] := Calc_Median(xVector);
yMedians[2] := Calc_Median(yVector);
xMedians[2] := VecMedian(xVector);
yMedians[2] := VecMedian(yVector);
end;

4
applications/lazstats/source/forms/analysis/descriptive/stemleafunit.pas
View File

@ -48,7 +48,7 @@ implementation
{$R *.lfm}
uses
Utils, MathUnit, GridProcs;
Utils, MatrixUnit, GridProcs;
{ TStemLeafForm }
@ -132,7 +132,7 @@ begin
nCases := Length(values);
SetLength(valueString, nCases);
Calc_MaxMin(values, max, min);
VecMaxMin(values, max, min);
minSize := MaxInt;
maxSize := -Maxint;

2
applications/lazstats/source/forms/analysis/descriptive/xvsmultyunit.pas
View File

@ -63,7 +63,7 @@ implementation
uses
TAChartUtils,
MainUnit, MatrixLib, MathUnit, GridProcs, Utils;
MainUnit, MatrixLib, MatrixUnit, GridProcs, Utils;
{ TXvsMultYForm }

2
applications/lazstats/source/forms/analysis/measurement_programs/gradebookunit.pas
View File

@ -105,7 +105,7 @@ implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit;
{ TGradebookFrm }

3
applications/lazstats/source/forms/analysis/measurement_programs/gradingunit.pas
View File

@ -68,11 +68,12 @@ type
var
GradingFrm: TGradingFrm;
implementation
uses
Math,
Utils, MathUnit, GradebookUnit;
Utils, MatrixUnit, GradebookUnit;
{ TGradingFrm }

2
applications/lazstats/source/forms/analysis/measurement_programs/testscoreunit.pas
View File

@ -125,7 +125,7 @@ implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit, MathUnit;
{ TTestScoreFrm }

129
applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas
View File

@ -8,7 +8,7 @@ uses
Classes, SysUtils, FileUtil, Forms, Controls, Graphics, Dialogs,
StdCtrls, Buttons, ExtCtrls, ComCtrls,
Globals, MainUnit, DictionaryUnit, Matrixlib, DataProcs,
MathUnit, ReportFrameUnit, ChartFrameUnit, BasicStatsParamsFormUnit;
RegressionUnit, ReportFrameUnit, ChartFrameUnit, BasicStatsParamsFormUnit;
type
@ -114,7 +114,7 @@ implementation
uses
Math,
TAChartUtils, TAChartAxisUtils, TALegend, TASources, TACustomSeries,
Utils, GridProcs;
Utils, MatrixUnit, GridProcs;
{ TWLSFrm }
@ -215,7 +215,7 @@ var
i, j, noIndep, depCol, weightCol, oldDepCol, NCases, pos, col: integer;
IndepCols: IntDyneVec = nil;
RowLabels: StrDyneVec = nil;
X, Y: double;
X, weight: double;
Means: DblDyneVec = nil;
Variances: DblDyneVec = nil;
StdDevs: DblDyneVec = nil;
@ -224,7 +224,6 @@ var
BStdErrs: DblDyneVec = nil;
BtTests: DblDyneVec = nil;
tProbs: DblDyneVec = nil;
predicted: Double;
lReport: TStrings;
StdErrEst: Double = 0.0;
R2: Double = 0.0;
@ -259,7 +258,7 @@ begin
depCol := NoVariables;
Process_SquaredResidualsRegression(indepCols, noIndep, depCol, RowLabels,
BWeights, nCases, printDesc);
if WeightChk.Checked then
if SaveChk.Checked then //WeightChk.Checked then
AddPredictedStuffToGrid(indepCols, noIndep, BWeights);
// *** Display squared residuals for each independent variable ***
@ -273,13 +272,13 @@ begin
IndepCols[Noindep] := DepCol;
for i := 1 to NoCases do
begin
Y := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[NoVariables,i])); // weight
weight := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[NoVariables,i]));
for j := 0 to Noindep do
begin
pos := IndepCols[j];
X := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
X := X * Y;
OS3MainFrm.DataGrid.Cells[pos,i] := FloatToStr(X); // wp: DON'T OVERWRITE GRID CELLS
X := X * weight;
OS3MainFrm.DataGrid.Cells[pos, i] := FloatToStr(X); // wp: DON'T OVERWRITE GRID CELLS
end;
end;
@ -331,16 +330,16 @@ begin
indepCols[Noindep] := depCol; // wp: CALCULATION SHOULD NORMALIZE USER WEIGHTS HERE !!!
for i := 1 to NoCases do
begin
Y := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[weightCol,i])); // weight
weight := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[weightCol, i]));
for j := 0 to Noindep do
begin
pos := indepCols[j];
X := StrToFloat(OS3MainFrm.DataGrid.Cells[pos,i]);
X := X * Y;
OS3MainFrm.DataGrid.Cells[pos,i] := FloatToStr(X);
X := X * weight;
OS3MainFrm.DataGrid.Cells[pos, i] := FloatToStr(X);
end;
end;
if (Origin2Chk.Checked) then // get means of variables and subtract from the values
if Origin2Chk.Checked then // get means of variables and subtract from the values
begin
for j := 0 to Noindep do
begin
@ -351,7 +350,7 @@ begin
begin
if (DataProcs.ValidValue(i,pos)) then
begin
Means[j] := Means[j] + StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos,i]));
Means[j] := Means[j] + StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[pos, i]));
NCases := NCases + 1;
end;
end;
@ -565,7 +564,7 @@ begin
SortOnX(xPoints, yPoints);
// Regression
Calc_BivariateRegression(xPoints, yPoints, AConfLevel, regressionRes);
BivariateRegression(xPoints, yPoints, AConfLevel, regressionRes);
// Create tab with chart and report controls
CreateOrGetChartFrame(yCol, yLabel, memo, chartFrame);
@ -652,22 +651,33 @@ end;
procedure TWLSFrm.PredictIt(ColNoSelected: IntDyneVec; NoVars: integer;
Means, StdDevs, BetaWeights: DblDyneVec; StdErrEst: double; NoIndepVars: integer);
var
col, i, j, k, Index: integer;
col1, col2, col3, col4, col5, i, j, k, Index: integer;
predicted, zpredicted, z1, z2, resid, residsqr: double;
colName: String;
begin
col := NoVariables + 1;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.DictGrid.ColCount := 8;
DictionaryFrm.NewVar(col);
OS3MainFrm.DataGrid.Cells[col,0] := 'Pred.z';
DictionaryFrm.DictGrid.Cells[1,col] := 'Pred.z';
colName := 'Pred.z.';
col1 := GetVariableIndex(OS3MainFrm.DataGrid, colName);
if col1 = -1 then
begin
col1 := NoVariables + 1;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.DictGrid.ColCount := 8;
DictionaryFrm.NewVar(col1);
OS3MainFrm.DataGrid.Cells[col1, 0] := colName;
DictionaryFrm.DictGrid.Cells[1, col1] := colName;
end;
colName := 'z Resid.';
col2 := GetVariableIndex(OS3MainFrm.DataGrid, colName);
if col2 = -1 then
begin
col2 := NoVariables + 1;
DictionaryFrm.NewVar(col2);
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
OS3MainFrm.DataGrid.Cells[col2, 0] := colName;
DictionaryFrm.DictGrid.Cells[1, col2] := colName;
end;
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
OS3MainFrm.DataGrid.Cells[col,0] := 'z Resid.';
DictionaryFrm.DictGrid.Cells[1,col] := 'z Resid.';
OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 2;
for i := 1 to NoCases do
begin
zpredicted := 0.0;
@ -677,56 +687,73 @@ begin
z1 := (StrToFloat(OS3MainFrm.DataGrid.Cells[k,i]) - Means[j]) / StdDevs[j];
zpredicted := zpredicted + (z1 * BetaWeights[j]);
end;
OS3MainFrm.DataGrid.Cells[col-1,i] := Format('%.4f',[zpredicted]);
OS3MainFrm.DataGrid.Cells[col1, i] := Format('%.4f',[zpredicted]);
if StdDevs[NoVars-1] <> 0.0 then
begin
Index := ColNoSelected[NoVars-1];
z2 := StrToFloat(OS3MainFrm.DataGrid.Cells[Index,i]);
z2 := (z2 - Means[NoVars-1]) / StdDevs[NoVars-1]; // z score
OS3MainFrm.DataGrid.Cells[col,i] := Format('%.4f',[z2 - zpredicted]); // z residual
OS3MainFrm.DataGrid.Cells[col2, i] := Format('%.4f',[z2 - zpredicted]); // z residual
end;
end;
// Get raw predicted and residuals
col := NoVariables + 1;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Pred.Raw';
OS3MainFrm.DataGrid.Cells[col,0] := 'Pred.Raw';
colName := 'Pred.Raw';
col3 := GetVariableIndex(OS3MainFrm.DataGrid, colName);
if col3 = -1 then
begin
col3 := NoVariables + 1;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.NewVar(col3);
DictionaryFrm.DictGrid.Cells[1, col3] := colName;
OS3MainFrm.DataGrid.Cells[col3, 0] := colName;
OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 1;
end;
// calculate raw predicted scores and store in DataGrid at col
for i := 1 to NoCases do
begin // predicted raw obtained from previously predicted z score
predicted := StrToFloat(OS3MainFrm.DataGrid.Cells[col-2,i]) * StdDevs[NoVars-1] + Means[NoVars-1];
OS3MainFrm.DataGrid.Cells[col,i] := Format('%.3f',[predicted]);
predicted := StrToFloat(OS3MainFrm.DataGrid.Cells[col1, i]) * StdDevs[NoVars-1] + Means[NoVars-1];
OS3MainFrm.DataGrid.Cells[col3, i] := Format('%.3f',[predicted]);
end;
// Calculate residuals of predicted raw scores end;
col := NoVariables +1;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Raw Resid.';
OS3MainFrm.DataGrid.Cells[col,0] := 'Raw Resid.';
colName := 'Raw Resid.';
col4 := GetVariableIndex(OS3MainFrm.DataGrid, colName);
if col4 = -1 then
begin
col4 := NoVariables + 1;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.NewVar(col4);
DictionaryFrm.DictGrid.Cells[1, col4] := colName;
OS3MainFrm.DataGrid.Cells[col4, 0] := colName;
OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 1;
end;
for i := 1 to NoCases do
begin
Index := ColNoSelected[NoVars-1];
resid := StrToFloat(OS3MainFrm.DataGrid.Cells[col-1,i]) - StrToFloat(OS3MainFrm.DataGrid.Cells[Index,i]);
OS3MainFrm.DataGrid.Cells[col,i] := Format('%.3f',[resid]);
resid := StrToFloat(OS3MainFrm.DataGrid.Cells[col3, i]) - StrToFloat(OS3MainFrm.DataGrid.Cells[Index, i]);
OS3MainFrm.DataGrid.Cells[col4, i] := Format('%.3f', [resid]);
end;
// get square of raw residuals
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.DictGrid.Cells[1,col] := 'ResidSqr';
OS3MainFrm.DataGrid.Cells[col,0] := 'ResidSqr';
colName := 'Resid Sqr';
col5 := GetVariableIndex(OS3MainFrm.DataGrid, colName);
if col5 = -1 then
begin
col5 := NoVariables + 1;
DictionaryFrm.NewVar(col5);
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
DictionaryFrm.DictGrid.Cells[1, col5] := colName;
OS3MainFrm.DataGrid.Cells[col5, 0] := colName;
OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 1;
end;
for i := 1 to NoCases do
begin
residsqr := StrToFloat(OS3MainFrm.DataGrid.Cells[col-1,i]);
residsqr := StrToFloat(OS3MainFrm.DataGrid.Cells[col4,i]);
residsqr := residsqr * residsqr;
OS3MainFrm.DataGrid.Cells[col,i] := Format('%.3f',[residsqr]);
OS3MainFrm.DataGrid.Cells[col5,i] := Format('%.3f',[residsqr]);
end;
end;

2
applications/lazstats/source/forms/analysis/multivariate/kmeansunit.pas
View File

@ -86,7 +86,7 @@ implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit;
{ TKMeansFrm }

2
applications/lazstats/source/forms/analysis/multivariate/medianpolishunit.pas
View File

@ -84,7 +84,7 @@ implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit;
{ TMedianPolishForm }

2
applications/lazstats/source/forms/analysis/multivariate/pathunit.pas
View File

@ -85,7 +85,7 @@ implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit, MathUnit;
{ TPathFrm }

3
applications/lazstats/source/forms/analysis/multivariate/singlelinkunit.pas
View File

@ -58,11 +58,12 @@ type
var
SingleLinkFrm: TSingleLinkFrm;
implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit;
{ TSingleLinkFrm }

3
applications/lazstats/source/forms/analysis/nonparametric/sensunit.pas
View File

@ -62,11 +62,12 @@ type
var
SensForm: TSensForm;
implementation
uses
Math,
Utils, MathUnit;
Utils, MatrixUnit;
{ TSensForm }

46
applications/lazstats/source/forms/misc/dictionaryunit.pas
View File

@ -49,7 +49,7 @@ type
public
{ public declarations }
procedure DelRow(row : integer);
procedure NewVar(row : integer);
procedure NewVar(ARow: integer);
procedure PasteVar(row : integer);
procedure CopyVar(row : integer);
procedure Init;
@ -331,31 +331,31 @@ begin
end;
//-------------------------------------------------------------------
procedure TDictionaryFrm.NewVar(row : integer);
procedure TDictionaryFrm.NewVar(ARow: integer);
var
i, j : integer;
i, j: integer;
begin
DictGrid.RowCount := DictGrid.RowCount + 1; // add new row
NoVariables := NoVariables + 1;
if (row < NoVariables) AND (NoVariables > 1) then // move current rows down 1
begin
for i := NoVariables downto row + 1 do
begin
for j := 1 to 7 do DictGrid.Cells[j,i] := DictGrid.Cells[j,i-1];
VarDefined[i] := VarDefined[i-1];
end;
end;
// put default values in new variable
Defaults(Self,row);
VarDefined[row] := true;
DictGrid.RowCount := DictGrid.RowCount + 1; // add new row
NoVariables := NoVariables + 1;
if (ARow < NoVariables) and (NoVariables > 1) then // move current rows down 1
begin
for i := NoVariables downto ARow + 1 do
for j := 1 to 7 do DictGrid.Cells[j,i] := DictGrid.Cells[j,i-1];
VarDefined[i] := VarDefined[i-1];
end;
// add to grid if grid column does not exist
if OS3MainFrm.DataGrid.ColCount < row then
begin
OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 1;
OS3MainFrm.DataGrid.Cells[row,0] := DictGrid.Cells[1,row];
end;
ReturnBtnClick(Self);
// put default values in new variable
Defaults(Self, ARow);
VarDefined[ARow] := true;
// add to grid if grid column does not exist
if OS3MainFrm.DataGrid.ColCount < ARow then
begin
OS3MainFrm.DataGrid.ColCount := OS3MainFrm.DataGrid.ColCount + 1;
OS3MainFrm.DataGrid.Cells[ARow,0] := DictGrid.Cells[1, ARow];
end;
ReturnBtnClick(Self);
end;
//-------------------------------------------------------------------

2
applications/lazstats/source/units/anovatestsunit.pas
View File

@ -104,7 +104,7 @@ procedure HomogeneityTest(
implementation
uses
Utils, MathUnit;
Utils, MatrixUnit, MathUnit;
procedure Tukey(error_ms : double; { mean squared for residual }
error_df : double; { deg. freedom for residual }

2
applications/lazstats/source/units/dataprocs.pas
View File

@ -61,7 +61,7 @@ function StringsToInt(StrCol: integer; var NewCol: integer; Prompt: boolean): bo
implementation
uses
Utils, MathUnit, MainUnit;
Utils, MatrixUnit, MainUnit;
// NOTE: Do not call GridProcs.GoodRecord here because this old function may

303
applications/lazstats/source/units/mathunit.pas
View File

@ -3,7 +3,6 @@ unit MathUnit;
{ extract some math functions from functionslib for easier testing }
{$mode objfpc}{$H+}
{$modeswitch advancedrecords}
interface
@ -47,46 +46,11 @@ function FactorialLn(n: Integer): Double;
function PoissonPDF(n: integer; a: double): Double;
function PoissonCDF(n: Integer; a: double): Double;
procedure Calc_MaxMin(const AData: DblDyneVec; out AMax, AMin: Double);
procedure Calc_MeanStdDev(const AData: DblDyneVec; out AMean, AStdDev: Double);
procedure Calc_MeanVarStdDev(const AData: DblDyneVec; out AMean, AVariance, AStdDev: Double);
procedure Calc_MeanVarStdDevSS(const AData: DblDyneVec; out AMean, AVariance, AStdDev, ASumOfSquares: Double);
procedure Calc_SumSS(const AData: DblDyneVec; out Sum, SS: Double);
function Calc_Median(const AData: DblDyneVec): Double;
type
TBivariateRegressionResults = record
Slope, Intercept: Double;
XMean, YMean: Double;
XVariance, YVariance: Double;
XStdDev, YStdDev: Double;
StdErrorPredicted: Double;
R: Double;
t: Double;
SXX, SXY, SYY: Double;
Count, DF: Integer;
function ConfidenceLimits(x: Double; Upper: Boolean): Double;
end;
procedure Calc_BivariateRegression(const xData, yData: DblDyneVec; AConfLevel: Double;
out AResults: TBivariateRegressionResults);
procedure Exchange(var a, b: Double); overload;
procedure Exchange(var a, b: Integer); overload;
procedure Exchange(var a, b: String); overload;
procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil); // not 100% tested...
implementation
uses
Math;
// Utils;
// Calculates the error function
// /x
@ -600,273 +564,6 @@ end;
{===============================================================================
* Vector-based calculations
===============================================================================}
procedure Calc_MaxMin(const AData: DblDyneVec; out AMax, AMin: Double);
var
i: Integer;
begin
AMin := Infinity;
AMax := -Infinity;
for i := Low(AData) to High(AData) do
begin
if AData[i] < AMin then AMin := AData[i];
if AData[i] > AMax then AMax := AData[i];
end;
end;
procedure Calc_MeanStdDev(const AData: DblDyneVec; out AMean, AStdDev: Double);
var
variance: Double;
begin
Calc_MeanVarStdDev(AData, AMean, variance, AStdDev);
end;
procedure Calc_MeanVarStdDev(const AData: DblDyneVec; out AMean, AVariance, AStdDev: Double);
var
sum, ss: Double;
n: Integer;
begin
AMean := NaN;
AVariance := NaN;
AStdDev := NaN;
n := Length(AData);
if n = 0 then
exit;
Calc_SumSS(AData, sum, ss);
AMean := sum / n;
if n = 1 then
exit;
AVariance := ((ss - sqr(AMean)) / n) / (n - 1);
AStdDev := sqrt(AVariance);
end;
procedure Calc_MeanVarStdDevSS(const AData: DblDyneVec;
out AMean, AVariance, AStdDev, ASumOfSquares: Double);
var
sum: Double;
n: Integer;
begin
AMean := NaN;
AVariance := NaN;
AStdDev := NaN;
n := Length(AData);
if n = 0 then
exit;
Calc_SumSS(AData, sum, ASumOfSquares);
AMean := sum / n;
if n = 1 then
exit;
AVariance := (ASumOfSquares - sqr(sum) / n) / (n - 1);
AStdDev := sqrt(AVariance);
end;
procedure Calc_SumSS(const AData: DblDyneVec; out Sum, SS: Double);
var
i: Integer;
begin
Sum := 0;
SS := 0;
for i := Low(AData) to High(AData) do
begin
Sum := Sum + AData[i];
SS := SS + sqr(AData[i]);
end;
end;
function Calc_Median(const AData: DblDyneVec): Double;
var
N, midPt: integer;
begin
SortOnX(AData);
N := Length(AData);
midPt := N div 2;
if odd(N) then
Result := AData[midPt] // odd no. of values
else
Result := (AData[midPt-1] + AData[midPt]) / 2; // even no. of values
end;
// It is assumed that xData and yData contain at least 3 elements and
// have the same count of elements.
procedure Calc_BivariateRegression(const xData, yData: DblDyneVec; AConfLevel: Double;
out AResults: TBivariateRegressionResults);
var
i: Integer;
begin
with AResults do
begin
Count := Length(xData);
// Calculate means, variances, stddevs
Calc_MeanVarStdDevSS(xData, XMean, XVariance, XStdDev, SXX);
Calc_MeanVarStdDevSS(yData, YMean, YVariance, YStdDev, SYY);
SXY := 0;
for i := 0 to Count-1 do
SXY := SXY + xData[i] * yData[i];
R := (SXY - XMean * YMean * Count) / ((Count - 1) * XStdDev * YStdDev);
StdErrorPredicted := sqrt(1.0 - sqr(R)) * YStdDev * sqrt((Count - 1) / (Count - 2));
Slope := R * YStdDev / XStdDev;
Intercept := YMean - Slope * XMean;
DF := Count - 2;
t := InverseT(AConfLevel, DF);
end;
end;
function TBivariateRegressionResults.ConfidenceLimits(x: Double; Upper: Boolean): Double;
var
yPred, seData: Double;
begin
yPred := Intercept + Slope * x;
seData := StdErrorPredicted * sqrt(1.0 + 1/Count + sqr(x - XMean)/SXX);
if Upper then
Result := yPred + t*seData
else
Result := yPred - t*seData;
end;
procedure Exchange(var a, b: Double);
var
tmp: Double;
begin
tmp := a;
a := b;
b := tmp;
end;
procedure Exchange(var a, b: Integer);
var
tmp: Integer;
begin
tmp := a;
a := b;
b := tmp;
end;
procedure Exchange(var a, b: String);
var
tmp: String;
begin
tmp := a;
a := b;
b := tmp;
end;
procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
var
i, j, N: Integer;
begin
N := Length(X);
if (Y <> nil) and (N <> Length(Y)) then
raise Exception.Create('[SortOnX] Arrays must have the same length.');
if (Z <> nil) and (N <> Length(Z)) then
raise Exception.Create('[SortOnX] Arrays must have the same length.');
for i := 0 to N - 2 do
begin
for j := i + 1 to N - 1 do
begin
if X[i] > X[j] then //swap
begin
Exchange(X[i], X[j]);
if Y <> nil then
Exchange(Y[i], Y[j]);
if Z <> nil then
Exchange(Z[i], Z[j]);
end;
end;
end;
end;
// NOTE: The matrix Y is transposed relative to the typical usage in LazStats
procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
var
i, j, k, N, Ny: Integer;
begin
N := Length(X);
if N <> Length(Y[0]) then
raise Exception.Create('[SortOnX] Arrays X and Y (2nd index) must have the same length');
Ny := Length(Y);
for i := 0 to N-2 do
begin
for j := i+1 to N-1 do
if X[i] > X[j] then
begin
Exchange(X[i], X[j]);
for k := 0 to Ny-1 do
Exchange(Y[k, i], Y[k, j]);
end;
end;
end;
procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
procedure DoQuickSort(L, R: Integer);
var
I,J: Integer;
P: Integer;
begin
repeat
I := L;
J := R;
P := (L+R) div 2;
repeat
while CompareValue(X[P], X[I]) > 0 do inc(I);
while CompareValue(X[P], X[J]) < 0 do dec(J);
if I <= J then begin
if I <> J then begin
Exchange(X[I], X[J]);
if Y <> nil then
Exchange(Y[I], Y[J]);
if Z <> nil then
Exchange(Z[I], Z[J]);
end;
if P = I then
P := J
else if P = J then
P := I;
inc(I);
dec(J);
end;
until I > J;
if L < J then
DoQuickSort(L, J);
L := I;
until I >= R;
end;
begin
DoQuickSort(0, High(X));
end;
initialization
InitFactLn();

2
applications/lazstats/source/units/matrixlib.pas
View File

@ -147,7 +147,7 @@ implementation
uses
Math, StrUtils,
MathUnit;
MatrixUnit, MathUnit;
procedure GridDotProd(col1, col2: integer; out Product: double; var Ngood: integer);
// Get the cross-product of two vectors

724
applications/lazstats/source/units/matrixunit.pas Normal file
View File

@ -0,0 +1,724 @@
{ A general-purpose matrix library }
unit MatrixUnit;
{$mode objfpc}{$H+}
interface
uses
Classes, SysUtils,
Globals;
type
TDblMatrix = DblDyneMat;
TDblVector = DblDyneVec;
EMatrix = class(Exception);
// Vectors
procedure VecMaxMin(const AData: TDblVector;
out AMax, AMin: Double);
procedure VecMeanStdDev(const AData: TDblVector;
out AMean, AStdDev: Double);
procedure VecMeanVarStdDev(const AData: TDblVector;
out AMean, AVariance, AStdDev: Double);
procedure VecMeanVarStdDevSS(const AData: TDblVector;
out AMean, AVariance, AStdDev, ASumOfSquares: Double);
procedure VecSumSS(const AData: TDblVector;
out Sum, SS: Double);
function VecMedian(const AData: TDblVector): Double;
// Matrices
operator + (A, B: TDblMatrix): TDblMatrix;
operator - (A, B: TDblMatrix): TDblMatrix;
operator * (A, B: TDblMatrix): TDblMatrix;
{ NOTE: Indices follow math convention:
- 1st index is the row index, i.e. runs vertically
- 2nd index is the col index, i.e. runs horizontally
All indices are 0-based. }
function MatAppendColVector(A: TDblMatrix; v: TDblVector): TDblMatrix;
procedure MatCheck(A: TDblMatrix);
procedure MatCheckSquare(A: TDblMatrix; out n: Integer);
function MatCopy(A: TDblMatrix): TDblMatrix;
function MatDeterminant(A: TDblMatrix): double;
function MatEqualSize(A, B: TDblMatrix): Boolean;
procedure MatExchangeElements(A: TDblMatrix; i1, j1, i2, j2: Integer);
procedure MatExchangeRows(A: TDblMatrix; i1, i2: Integer);
function MatInverse(A: TDblMatrix): TDblMatrix;
function MatIsSquare(A: TDblMatrix; out n: Integer): Boolean;
function MatNumCols(A: TDblMatrix): Integer;
function MatNumRows(A: TDblMatrix): Integer;
procedure MatSize(A: TDblMatrix; out n, m: Integer);
function MatTranspose(A: TDblMatrix): TDblMatrix;
// Sorting
procedure Exchange(var a, b: Double); overload;
procedure Exchange(var a, b: Integer); overload;
procedure Exchange(var a, b: String); overload;
procedure SortOnX(X: TDblVector; Y: TDblVector = nil; Z: TDblVector = nil);
procedure SortOnX(X: TDblVector; Y: TDblMatrix);
procedure QuickSortOnX(X: TDblVector; Y: TDblVector = nil; Z: TDblVector = nil); // not 100% tested...
type
TLUSolver = class
private
FLU: TDblMatrix; // LU-decomposed matrix
FIndex: array of integer; // records permutations
d: Double; // odd or even row permuations
procedure LUDecompose;
public
constructor Create(A: TDblMatrix);
// function CreateL: TDblMatrix;
// function CreateU: TDblMatrix;
procedure Solve(b: TDblVector);
property LU: TDblMatrix read FLU;
end;
implementation
uses
Math;
{===============================================================================
Statistical vector operations
===============================================================================}
procedure VecMaxMin(const AData: TDblVector; out AMax, AMin: Double);
var
i: Integer;
begin
AMin := Infinity;
AMax := -Infinity;
for i := Low(AData) to High(AData) do
begin
if AData[i] < AMin then AMin := AData[i];
if AData[i] > AMax then AMax := AData[i];
end;
end;
procedure VecMeanStdDev(const AData: TDblVector; out AMean, AStdDev: Double);
var
variance: Double;
begin
VecMeanVarStdDev(AData, AMean, variance, AStdDev);
end;
procedure VecMeanVarStdDev(const AData: TDblVector;
out AMean, AVariance, AStdDev: Double);
var
sum, ss: Double;
n: Integer;
begin
AMean := NaN;
AVariance := NaN;
AStdDev := NaN;
n := Length(AData);
if n = 0 then
exit;
VecSumSS(AData, sum, ss);
AMean := sum / n;
if n = 1 then
exit;
AVariance := ((ss - sqr(AMean)) / n) / (n - 1);
AStdDev := sqrt(AVariance);
end;
procedure VecMeanVarStdDevSS(const AData: TDblVector;
out AMean, AVariance, AStdDev, ASumOfSquares: Double);
var
sum: Double;
n: Integer;
begin
AMean := NaN;
AVariance := NaN;
AStdDev := NaN;
n := Length(AData);
if n = 0 then
exit;
VecSumSS(AData, sum, ASumOfSquares);
AMean := sum / n;
if n = 1 then
exit;
AVariance := (ASumOfSquares - sqr(sum) / n) / (n - 1);
AStdDev := sqrt(AVariance);
end;
procedure VecSumSS(const AData: TDblVector; out Sum, SS: Double);
var
i: Integer;
begin
Sum := 0;
SS := 0;
for i := Low(AData) to High(AData) do
begin
Sum := Sum + AData[i];
SS := SS + sqr(AData[i]);
end;
end;
function VecMedian(const AData: TDblVector): Double;
var
N, midPt: integer;
begin
SortOnX(AData);
N := Length(AData);
midPt := N div 2;
if odd(N) then
Result := AData[midPt] // odd no. of values
else
Result := (AData[midPt-1] + AData[midPt]) / 2; // even no. of values
end;
{===============================================================================
Matrix operators
===============================================================================}
operator + (A, B: TDblMatrix): TDblMatrix;
var
n, m, i, j: Integer;
begin
n := Length(A);
m := Length(A[0]);
if (n <> Length(B)) or (m <> Length(B[0])) then
raise EMatrix.Create('Matrix subtraction: dimension error');
SetLength(Result, n,m);
for i := 0 to n-1 do
for j := 0 to m-1 do
Result[i, j] := A[i, j] + B[i, j];
end;
operator - (A, B: TDblMatrix): TDblMatrix;
var
n, m, i, j: Integer;
begin
n := Length(A);
m := Length(A[0]);
if (n <> Length(B)) or (m <> Length(B[0])) then
raise EMatrix.Create('Matrix subtraction: dimension error');
SetLength(Result, n,m);
for i := 0 to n-1 do
for j := 0 to m-1 do
Result[i, j] := A[i, j] - B[i, j];
end;
operator * (A, B: TDblMatrix): TDblMatrix;
var
na, ma, nb, mb, i, j, k: Integer;
sum: Double;
begin
MatSize(A, na, ma);
MatSize(B, nb, mb);
if ma <> nb then
raise EMatrix.Create('Matrix product: dimension error');
SetLength(Result, na,mb);
for i := 0 to na-1 do
for j := 0 to mb-1 do begin
sum := 0;
for k := 0 to mA - 1 do
sum := sum + A[i, k] * B[k, j];
Result[i, j] := sum;
end;
end;
{ Adds the elements of the vector v to the rows of the matrix A, i.e. the number
of columns increases by 1 }
function MatAppendColVector(A: TDblMatrix; v: TDblVector): TDblMatrix;
var
i, j, n, m, nv: Integer;
begin
MatSize(A, n, m);
nv := Length(v);
if n <> nv then
raise EMatrix.Create('Dimension error.');
SetLength(Result, n, m+1);
for i := 0 to n-1 do
begin
for j := 0 to m-1 do
Result[i, j] := A[i, j];
Result[i, m] := v[i];
end;
end;
{ Checks whether the matrix A has been initialized (i.e. by calling SetLength).
Raises an Exception otherwise. }
procedure MatCheck(A: TDblMatrix);
begin
if (A = nil) and (Length(A) = 0) then
raise EMatrix.Create('Matrix not created.');
end;
procedure MatCheckSquare(A: TDblMatrix; out n: Integer);
begin
if not MatIsSquare(A, n) then
raise EMatrix.Create('Matrix is not square.');
end;
function MatCopy(A: TDblMatrix): TDblMatrix;
var
n, m, i, j: Integer;
begin
MatSize(A, n,m);
SetLength(Result, n, m);
for i := 0 to n-1 do
for j := 0 to m-1 do
Result[i, j] := A[i, j];
end;
function MatDeterminant(A: TDblMatrix): double;
var
i, n: integer;
S: TLUSolver;
d: Double;
begin
Result := NaN;
MatCheck(A);
MatCheckSquare(A, n);
S := TLUSolver.Create(A);
try
d := S.d;
for i := 0 to n-1 do
d := d * S.LU[i, i];
Result := d;
finally
S.Free;
end;
end;
function MatEqualSize(A, B: TDblMatrix): Boolean;
var
na, ma, nb, mb: Integer;
begin
MatSize(A, na, ma);
MatSize(B, nb, mb);
Result := (na = ma) and (nb = mb);
end;
{ Exchanges the elements (i1,j1) and (i2,j2) of the matrix A }
procedure MatExchangeElements(A: TDblMatrix; i1,j1, i2,j2: Integer);
var
tmp: Double;
begin
tmp := A[i1, j1];
A[i1, j1] := A[i2, j2];
A[i2, j2] := tmp;
end;
{ Exchanges the rows i1 and i2 of matrix A }
procedure MatExchangeRows(A: TDblMatrix; i1, i2: Integer);
var
j, n, m: integer;
begin
MatSize(A, n,m);
for j := 0 to m-1 do
MatExchangeElements(A, i1,j, i2,j);
end;
function MatInverse(A: TDblMatrix): TDblMatrix;
var
i, j, n: integer;
S: TLUSolver;
v: TDblVector;
begin
MatCheck(A);
MatCheckSquare(A, n);
SetLength(Result, n,n);
SetLength(v, n);
S := TLUSolver.Create(A);
try
for j := 0 to n-1 do
begin
for i := 0 to n-1 do
v[i] := 0.0;
v[j] := 1.0;
S.Solve(v);
for i := 0 to n-1 do Result[i, j] := v[i];
end;
finally
S.Free;
end;
end;
function MatIsSquare(A: TDblMatrix; out n: Integer): Boolean;
var
m: Integer;
begin
MatSize(A, n, m);
Result := (n = m);
end;
function MatNumCols(A: TDblMatrix): Integer;
var
n: Integer;
begin
MatSize(A, n, Result);
end;
function MatNumRows(A: TDblMatrix): Integer;
var
m: Integer;
begin
MatSize(A, Result, m);
end;
procedure MatSize(A: TDblMatrix; out n, m: Integer);
begin
n := Length(A);
m := Length(A[0]);
end;
function MatTranspose(A: TDblMatrix): TDblMatrix;
var
n, i, j: Integer;
begin
MatCheck(A);
MatCheckSquare(A, n);
SetLength(Result, n, n);
for i := 0 to n-1 do
for j := 0 to n-1 do
Result[j, i] := A[i, j];
end;
{===============================================================================
TLUSolver
--------------------------------------------------------------------------------
Solves the linear system of equations, A x = b by means of LU decomposition.
Matrix A is passed when the solver is created and it is decomposed in
upper and lower triangular matrices, A = L * U where L has non-zero elements
only below, and U only above the diagonal.
The equation system is is solved by means of the insertion method (Solve).
Ref:
Press et al, Numerical Recipies in Pascal
===============================================================================}
constructor TLUSolver.Create(A: TDblMatrix);
var
n: Integer;
begin
MatCheckSquare(A, n);
FLU := MatCopy(A);
LUDecompose;
end;
(*
{ A has been decomposed into the product of two matrices, L*U. L and U are
saved in a minimal way in the single matrix FLU.
CreateL creates a new matrix and copies the elements of L to the correct
position in the result. The product of the two matrices created by CreateL
and CreateU, is A. }
function TLUSolver.CreateL: TDblMatrix;
var
n, m, i, j: integer;
begin
MatrixSize(FLU, n,m);
SetLength(Result, n, m);
for i := 0 to n-1 do begin
for j := 0 to i do begin
if (i > j) then
Result[i, j] := FLU[i,j]
else
if (i = j) then
Result[i, j] := 1.0;
end;
end;
for i := 0 to n - 1 do MatExchangeRows(Result, i, FIndex[i]);
end;
{ A has been decomposed into the product of two matrices, L*U. L and U are
saved in a minimal way in the single matrix FLU.
CreateU creates a new matrix and copies the elements of U to the correct
position in the result. The product of the two matrices created by CreateL
and CreateU, is A. }
function TLUSolver.CreateU: TDblMatrix;
var
n, m, i, j: integer;
begin
MatSize(FLU, n,m);
SetLength(Result, n, m);
for i := 0 to n-1 do
for j := i to m-1 do
if (i <= j) then Result[i,j] := FLU[i,j];
end;
*)
{ Executes the LU decompositon. Main procedure of the solver. }
procedure TLUSolver.LUDecompose;
const
tiny = 1.0e-20;
var
k, j, imax, i: integer;
sum, dum, big: float;
vv: TDblVector;
n : integer;
begin
n := MatNumRows(FLU);
SetLength(vv, n);
SetLength(FIndex, n);
d := 1.0;
for i := 0 to n-1 do // Scaling information
begin
big := 0.0;
for j := 0 to n-1 do
begin
dum := abs(FLU[i, j]);
if dum > big then big := dum;
end;
if big = 0.0 then
raise EMatrix.Create('Matrix is singular.');
vv[i] := 1.0 / big;
end;
for j := 0 to n-1 do begin
for i := 0 to j-1 do begin
sum := FLU[i, j];
for k := 0 to i-1 do
sum := sum - FLU[i, k] * FLU[k, j];
FLU[i, j] := sum;
end;
big := 0.0;
for i := j to n-1 do begin
sum := FLU[i,j];
for k := 0 to j-1 do
sum := sum - FLU[i, k] * FLU[k, j];
FLU[i, j] := sum;
dum := vv[i] * abs(sum);
if dum > big then begin
big := dum;
imax := i;
end;
end;
if j <> imax then begin
for k := 0 to n-1 do begin
dum := FLU[imax, k];
FLU[imax, k] := FLU[j, k];
FLU[j, k] := dum;
end;
d := -d;
vv[imax] := vv[j];
end;
FIndex[j] := imax;
if FLU[j, j]=0.0 then FLU[j, j] := tiny;
if j <> n-1 then
begin
dum := 1.0 / FLU[j, j];
for i := succ(j) to n-1 do
FLU[i, j] := FLU[i, j] * dum;
end;
end;
end;
{ Solves the equation system A*x = b for x and returns the solution in b.
A already had been LU-decomposed when the solver was created. This means
that the method can be reused again and again for different b vectors.
NOTE: the input b is overwritten by the calculation! }
procedure TLUSolver.Solve(b: TDblVector);
var
nB, n, m: integer;
i, ii, ip, j: integer;
sum: Double;
begin
nB := Length(b);
MatSize(FLU, n,m);
if nB <> n then
raise EMatrix.Create('TLUSolver: Dimension error.');
ii := -1;
for i := 0 to n-1 do
begin
ip := FIndex[i];
sum := b[ip];
b[ip] := b[i];
if ii <> -1 then
for j := ii to pred(i) do
sum := sum - FLU[i, j] * b[j]
else
if (sum <> 0) then
ii := i;
b[i] := sum;
end;
for i := n-1 downto 0 do
begin
sum := b[i];
for j := succ(i) to n-1 do
sum := sum - FLU[i, j] * b[j];
b[i] := sum / FLU[i,i];
end;
end;
{===============================================================================
Sorting
===============================================================================}
procedure Exchange(var a, b: Double);
var
tmp: Double;
begin
tmp := a;
a := b;
b := tmp;
end;
procedure Exchange(var a, b: Integer);
var
tmp: Integer;
begin
tmp := a;
a := b;
b := tmp;
end;
procedure Exchange(var a, b: String);
var
tmp: String;
begin
tmp := a;
a := b;
b := tmp;
end;
procedure SortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
var
i, j, N: Integer;
begin
N := Length(X);
if (Y <> nil) and (N <> Length(Y)) then
raise Exception.Create('[SortOnX] Arrays must have the same length.');
if (Z <> nil) and (N <> Length(Z)) then
raise Exception.Create('[SortOnX] Arrays must have the same length.');
for i := 0 to N - 2 do
begin
for j := i + 1 to N - 1 do
begin
if X[i] > X[j] then //swap
begin
Exchange(X[i], X[j]);
if Y <> nil then
Exchange(Y[i], Y[j]);
if Z <> nil then
Exchange(Z[i], Z[j]);
end;
end;
end;
end;
// NOTE: The matrix Y is transposed relative to the typical usage in LazStats
procedure SortOnX(X: DblDyneVec; Y: DblDyneMat);
var
i, j, k, N, Ny: Integer;
begin
N := Length(X);
if N <> Length(Y[0]) then
raise Exception.Create('[SortOnX] Arrays X and Y (2nd index) must have the same length');
Ny := Length(Y);
for i := 0 to N-2 do
begin
for j := i+1 to N-1 do
if X[i] > X[j] then
begin
Exchange(X[i], X[j]);
for k := 0 to Ny-1 do
Exchange(Y[k, i], Y[k, j]);
end;
end;
end;
procedure QuickSortOnX(X: DblDyneVec; Y: DblDyneVec = nil; Z: DblDyneVec = nil);
procedure DoQuickSort(L, R: Integer);
var
I,J: Integer;
P: Integer;
begin
repeat
I := L;
J := R;
P := (L+R) div 2;
repeat
while CompareValue(X[P], X[I]) > 0 do inc(I);
while CompareValue(X[P], X[J]) < 0 do dec(J);
if I <= J then begin
if I <> J then begin
Exchange(X[I], X[J]);
if Y <> nil then
Exchange(Y[I], Y[J]);
if Z <> nil then
Exchange(Z[I], Z[J]);
end;
if P = I then
P := J
else if P = J then
P := I;
inc(I);
dec(J);
end;
until I > J;
if L < J then
DoQuickSort(L, J);
L := I;
until I >= R;
end;
begin
DoQuickSort(0, High(X));
end;
end.

106
applications/lazstats/source/units/regressionunit.pas Normal file
View File

@ -0,0 +1,106 @@
unit RegressionUnit;
{$mode objfpc}{$H+}
{$modeswitch advancedrecords}
interface
uses
Classes, SysUtils,
MatrixUnit;
type
TBivariateRegressionResults = record
Slope, Intercept: Double;
XMean, YMean: Double;
XVariance, YVariance: Double;
XStdDev, YStdDev: Double;
StdErrorPredicted: Double;
R: Double;
t: Double;
SXX, SXY, SYY: Double;
Count, DF: Integer;
function ConfidenceLimits(x: Double; Upper: Boolean): Double;
end;
procedure BivariateRegression(const xData, yData: TDblVector;
AConfLevel: Double; out AResults: TBivariateRegressionResults);
type
TMultipleRegressionResults = record
end;
procedure MultipleRegression(const xData: TDblMatrix; const yData: TDblVector;
AConfLevel: Double; out AResults: TMultipleRegressionResults);
implementation
uses
MathUnit;
{===============================================================================
Bivariate regression
Calculates a bivariate regression y = a x + b
x is a vector of independent observations, y of dependent observations.
Both vectors must have the same length.
It is assumed that xData and yData contain at least 3 elements and
have the same count of elements.
===============================================================================}
procedure BivariateRegression(const xData, yData: TDblVector;
AConfLevel: Double; out AResults: TBivariateRegressionResults);
var
i: Integer;
begin
with AResults do
begin
Count := Length(xData);
// Calculate means, variances, stddevs
VecMeanVarStdDevSS(xData, XMean, XVariance, XStdDev, SXX);
VecMeanVarStdDevSS(yData, YMean, YVariance, YStdDev, SYY);
SXY := 0;
for i := 0 to Count-1 do
SXY := SXY + xData[i] * yData[i];
R := (SXY - XMean * YMean * Count) / ((Count - 1) * XStdDev * YStdDev);
StdErrorPredicted := sqrt(1.0 - sqr(R)) * YStdDev * sqrt((Count - 1) / (Count - 2));
Slope := R * YStdDev / XStdDev;
Intercept := YMean - Slope * XMean;
DF := Count - 2;
t := InverseT(AConfLevel, DF);
end;
end;
function TBivariateRegressionResults.ConfidenceLimits(x: Double; Upper: Boolean): Double;
var
yPred, seData: Double;
begin
yPred := Intercept + Slope * x;
seData := StdErrorPredicted * sqrt(1.0 + 1/Count + sqr(x - XMean)/SXX);
if Upper then
Result := yPred + t*seData
else
Result := yPred - t*seData;
end;
{===============================================================================
Multiple regression
===============================================================================}
procedure MultipleRegression(const xData: TDblMatrix; const yData: TDblVector;
AConfLevel: Double; out AResults: TMultipleRegressionResults);
begin
end;
end.