From af6047d214f0cd79ab505bf4b90eed5cecc04013 Mon Sep 17 00:00:00 2001 From: wp_xxyyzz Date: Wed, 14 Oct 2020 14:28:59 +0000 Subject: [PATCH] LazStats: Add and integrte new standalone units MatrixUnit and RegressionUnit. Move some code from MathUnit to MatrixUnit or RegressionUnit. git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7777 8e941d3f-bd1b-0410-a28a-d453659cc2b4 --- applications/lazstats/source/LazStats.lpi | 12 +- applications/lazstats/source/LazStats.lpr | 2 +- .../analysis/descriptive/descriptiveunit.pas | 6 +- .../analysis/descriptive/multxvsyunit.pas | 2 +- .../forms/analysis/descriptive/plotxyunit.pas | 6 +- .../descriptive/resistancelineunit.pas | 16 +- .../analysis/descriptive/stemleafunit.pas | 4 +- .../analysis/descriptive/xvsmultyunit.pas | 2 +- .../measurement_programs/gradebookunit.pas | 2 +- .../measurement_programs/gradingunit.pas | 3 +- .../measurement_programs/testscoreunit.pas | 2 +- .../analysis/multiple_regression/wlsunit.pas | 129 ++-- .../analysis/multivariate/kmeansunit.pas | 2 +- .../multivariate/medianpolishunit.pas | 2 +- .../forms/analysis/multivariate/pathunit.pas | 2 +- .../analysis/multivariate/singlelinkunit.pas | 3 +- .../forms/analysis/nonparametric/sensunit.pas | 3 +- .../source/forms/misc/dictionaryunit.pas | 46 +- .../lazstats/source/units/anovatestsunit.pas | 2 +- .../lazstats/source/units/dataprocs.pas | 2 +- .../lazstats/source/units/mathunit.pas | 303 -------- .../lazstats/source/units/matrixlib.pas | 2 +- .../lazstats/source/units/matrixunit.pas | 724 ++++++++++++++++++ .../lazstats/source/units/regressionunit.pas | 106 +++ 24 files changed, 975 insertions(+), 408 deletions(-) create mode 100644 applications/lazstats/source/units/matrixunit.pas create mode 100644 applications/lazstats/source/units/regressionunit.pas diff --git a/applications/lazstats/source/LazStats.lpi b/applications/lazstats/source/LazStats.lpi index c4abe63b7..06442da67 100644 --- a/applications/lazstats/source/LazStats.lpi +++ b/applications/lazstats/source/LazStats.lpi @@ -116,7 +116,7 @@ - + @@ -1533,6 +1533,16 @@ + + + + + + + + + + diff --git a/applications/lazstats/source/LazStats.lpr b/applications/lazstats/source/LazStats.lpr index 42339e2dc..0fee7fa55 100644 --- a/applications/lazstats/source/LazStats.lpr +++ b/applications/lazstats/source/LazStats.lpr @@ -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} diff --git a/applications/lazstats/source/forms/analysis/descriptive/descriptiveunit.pas b/applications/lazstats/source/forms/analysis/descriptive/descriptiveunit.pas index 5a8ed33b6..54fd7463a 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/descriptiveunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/descriptiveunit.pas @@ -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); diff --git a/applications/lazstats/source/forms/analysis/descriptive/multxvsyunit.pas b/applications/lazstats/source/forms/analysis/descriptive/multxvsyunit.pas index 6a56bd1e9..753eb1e2e 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/multxvsyunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/multxvsyunit.pas @@ -67,7 +67,7 @@ implementation uses TATypes, - Math, Utils, MathUnit; + Math, Utils, MatrixUnit; { TMultXvsYFrm } diff --git a/applications/lazstats/source/forms/analysis/descriptive/plotxyunit.pas b/applications/lazstats/source/forms/analysis/descriptive/plotxyunit.pas index 946a5916b..d92f8c53d 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/plotxyunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/plotxyunit.pas @@ -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); diff --git a/applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas b/applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas index 05ee81a55..313dcdec0 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/resistancelineunit.pas @@ -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; diff --git a/applications/lazstats/source/forms/analysis/descriptive/stemleafunit.pas b/applications/lazstats/source/forms/analysis/descriptive/stemleafunit.pas index 45d547930..9c77021dd 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/stemleafunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/stemleafunit.pas @@ -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; diff --git a/applications/lazstats/source/forms/analysis/descriptive/xvsmultyunit.pas b/applications/lazstats/source/forms/analysis/descriptive/xvsmultyunit.pas index 6eecee104..77b33f657 100644 --- a/applications/lazstats/source/forms/analysis/descriptive/xvsmultyunit.pas +++ b/applications/lazstats/source/forms/analysis/descriptive/xvsmultyunit.pas @@ -63,7 +63,7 @@ implementation uses TAChartUtils, - MainUnit, MatrixLib, MathUnit, GridProcs, Utils; + MainUnit, MatrixLib, MatrixUnit, GridProcs, Utils; { TXvsMultYForm } diff --git a/applications/lazstats/source/forms/analysis/measurement_programs/gradebookunit.pas b/applications/lazstats/source/forms/analysis/measurement_programs/gradebookunit.pas index defede370..926bd4d37 100644 --- a/applications/lazstats/source/forms/analysis/measurement_programs/gradebookunit.pas +++ b/applications/lazstats/source/forms/analysis/measurement_programs/gradebookunit.pas @@ -105,7 +105,7 @@ implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit; { TGradebookFrm } diff --git a/applications/lazstats/source/forms/analysis/measurement_programs/gradingunit.pas b/applications/lazstats/source/forms/analysis/measurement_programs/gradingunit.pas index 8e5e38563..22d5795d4 100644 --- a/applications/lazstats/source/forms/analysis/measurement_programs/gradingunit.pas +++ b/applications/lazstats/source/forms/analysis/measurement_programs/gradingunit.pas @@ -68,11 +68,12 @@ type var GradingFrm: TGradingFrm; + implementation uses Math, - Utils, MathUnit, GradebookUnit; + Utils, MatrixUnit, GradebookUnit; { TGradingFrm } diff --git a/applications/lazstats/source/forms/analysis/measurement_programs/testscoreunit.pas b/applications/lazstats/source/forms/analysis/measurement_programs/testscoreunit.pas index 139becaac..edd940f4a 100644 --- a/applications/lazstats/source/forms/analysis/measurement_programs/testscoreunit.pas +++ b/applications/lazstats/source/forms/analysis/measurement_programs/testscoreunit.pas @@ -125,7 +125,7 @@ implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit, MathUnit; { TTestScoreFrm } diff --git a/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas b/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas index d1b65434f..743e7371c 100644 --- a/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas +++ b/applications/lazstats/source/forms/analysis/multiple_regression/wlsunit.pas @@ -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; diff --git a/applications/lazstats/source/forms/analysis/multivariate/kmeansunit.pas b/applications/lazstats/source/forms/analysis/multivariate/kmeansunit.pas index a871f7d95..4d2a948bc 100644 --- a/applications/lazstats/source/forms/analysis/multivariate/kmeansunit.pas +++ b/applications/lazstats/source/forms/analysis/multivariate/kmeansunit.pas @@ -86,7 +86,7 @@ implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit; { TKMeansFrm } diff --git a/applications/lazstats/source/forms/analysis/multivariate/medianpolishunit.pas b/applications/lazstats/source/forms/analysis/multivariate/medianpolishunit.pas index b2a56a73a..cd388153e 100644 --- a/applications/lazstats/source/forms/analysis/multivariate/medianpolishunit.pas +++ b/applications/lazstats/source/forms/analysis/multivariate/medianpolishunit.pas @@ -84,7 +84,7 @@ implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit; { TMedianPolishForm } diff --git a/applications/lazstats/source/forms/analysis/multivariate/pathunit.pas b/applications/lazstats/source/forms/analysis/multivariate/pathunit.pas index eedd6ea2a..b88f0cdf3 100644 --- a/applications/lazstats/source/forms/analysis/multivariate/pathunit.pas +++ b/applications/lazstats/source/forms/analysis/multivariate/pathunit.pas @@ -85,7 +85,7 @@ implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit, MathUnit; { TPathFrm } diff --git a/applications/lazstats/source/forms/analysis/multivariate/singlelinkunit.pas b/applications/lazstats/source/forms/analysis/multivariate/singlelinkunit.pas index 09cdea0b7..86760ac9b 100644 --- a/applications/lazstats/source/forms/analysis/multivariate/singlelinkunit.pas +++ b/applications/lazstats/source/forms/analysis/multivariate/singlelinkunit.pas @@ -58,11 +58,12 @@ type var SingleLinkFrm: TSingleLinkFrm; + implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit; { TSingleLinkFrm } diff --git a/applications/lazstats/source/forms/analysis/nonparametric/sensunit.pas b/applications/lazstats/source/forms/analysis/nonparametric/sensunit.pas index 59611d7ea..d7044001e 100644 --- a/applications/lazstats/source/forms/analysis/nonparametric/sensunit.pas +++ b/applications/lazstats/source/forms/analysis/nonparametric/sensunit.pas @@ -62,11 +62,12 @@ type var SensForm: TSensForm; + implementation uses Math, - Utils, MathUnit; + Utils, MatrixUnit; { TSensForm } diff --git a/applications/lazstats/source/forms/misc/dictionaryunit.pas b/applications/lazstats/source/forms/misc/dictionaryunit.pas index e1ef99d84..308cb2af6 100644 --- a/applications/lazstats/source/forms/misc/dictionaryunit.pas +++ b/applications/lazstats/source/forms/misc/dictionaryunit.pas @@ -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; //------------------------------------------------------------------- diff --git a/applications/lazstats/source/units/anovatestsunit.pas b/applications/lazstats/source/units/anovatestsunit.pas index a78d977c7..f959b2acb 100644 --- a/applications/lazstats/source/units/anovatestsunit.pas +++ b/applications/lazstats/source/units/anovatestsunit.pas @@ -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 } diff --git a/applications/lazstats/source/units/dataprocs.pas b/applications/lazstats/source/units/dataprocs.pas index e12e316fa..ae55ffb7e 100644 --- a/applications/lazstats/source/units/dataprocs.pas +++ b/applications/lazstats/source/units/dataprocs.pas @@ -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 diff --git a/applications/lazstats/source/units/mathunit.pas b/applications/lazstats/source/units/mathunit.pas index 51d9459a1..1d145999a 100644 --- a/applications/lazstats/source/units/mathunit.pas +++ b/applications/lazstats/source/units/mathunit.pas @@ -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(); diff --git a/applications/lazstats/source/units/matrixlib.pas b/applications/lazstats/source/units/matrixlib.pas index 8c4a1ff24..b0e98c478 100644 --- a/applications/lazstats/source/units/matrixlib.pas +++ b/applications/lazstats/source/units/matrixlib.pas @@ -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 diff --git a/applications/lazstats/source/units/matrixunit.pas b/applications/lazstats/source/units/matrixunit.pas new file mode 100644 index 000000000..e68a09dd9 --- /dev/null +++ b/applications/lazstats/source/units/matrixunit.pas @@ -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. + diff --git a/applications/lazstats/source/units/regressionunit.pas b/applications/lazstats/source/units/regressionunit.pas new file mode 100644 index 000000000..031c6f5eb --- /dev/null +++ b/applications/lazstats/source/units/regressionunit.pas @@ -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. +