From 44d2ab491830e4e5eb05b2403947eed8c9ec7633 Mon Sep 17 00:00:00 2001 From: wp_xxyyzz Date: Tue, 15 Dec 2020 22:54:06 +0000 Subject: [PATCH] LazStats: Massive refactoring of AutoCorUnit regarding usability. git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7930 8e941d3f-bd1b-0410-a28a-d453659cc2b4 --- applications/lazstats/source/LazStats.lpi | 593 ++++--- applications/lazstats/source/LazStats.lpr | 3 +- .../analysis/correlation/autocorunit.lfm | 195 +-- .../analysis/correlation/autocorunit.pas | 1487 +++++++---------- .../analysis/correlation/expsmoothunit.pas | 8 +- .../analysis/correlation/smoothingunit.lfm | 491 ++++++ .../analysis/correlation/smoothingunit.pas | 1048 ++++++++++++ .../lazstats/source/units/mathunit.pas | 2 + .../lazstats/source/units/matrixunit.pas | 18 +- .../lazstats/source/units/regressionunit.pas | 9 +- 10 files changed, 2558 insertions(+), 1296 deletions(-) create mode 100644 applications/lazstats/source/forms/analysis/correlation/smoothingunit.lfm create mode 100644 applications/lazstats/source/forms/analysis/correlation/smoothingunit.pas diff --git a/applications/lazstats/source/LazStats.lpi b/applications/lazstats/source/LazStats.lpi index 51ea93554..9b36c775c 100644 --- a/applications/lazstats/source/LazStats.lpi +++ b/applications/lazstats/source/LazStats.lpi @@ -404,1145 +404,1145 @@ - - - - - - - - - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + + + + + + + + + @@ -1568,11 +1568,6 @@ - - - - - diff --git a/applications/lazstats/source/LazStats.lpr b/applications/lazstats/source/LazStats.lpr index 0fee7fa55..476c07935 100644 --- a/applications/lazstats/source/LazStats.lpr +++ b/applications/lazstats/source/LazStats.lpr @@ -8,7 +8,8 @@ uses {$ENDIF}{$ENDIF} Interfaces, // this includes the LCL widgetset Forms, tachartlazaruspkg, tachartprint, lhelpcontrolpkg, Globals, LicenseUnit, - OptionsUnit, MainDM, MainUnit, GridProcs, MatrixUnit, RegressionUnit; + OptionsUnit, MainDM, MainUnit, GridProcs, MatrixUnit, RegressionUnit, +SmoothingUnit; {$R LazStats.res} diff --git a/applications/lazstats/source/forms/analysis/correlation/autocorunit.lfm b/applications/lazstats/source/forms/analysis/correlation/autocorunit.lfm index 8380808fb..0aa625268 100644 --- a/applications/lazstats/source/forms/analysis/correlation/autocorunit.lfm +++ b/applications/lazstats/source/forms/analysis/correlation/autocorunit.lfm @@ -1,7 +1,7 @@ inherited AutoCorrForm: TAutoCorrForm - Left = 767 + Left = 446 Height = 512 - Top = 157 + Top = 176 Width = 991 HelpType = htKeyword HelpKeyword = 'html/Autocorrelation.htm' @@ -10,34 +10,34 @@ inherited AutoCorrForm: TAutoCorrForm ClientWidth = 991 inherited ParamsPanel: TPanel Height = 496 - Width = 448 + Width = 400 ClientHeight = 496 - ClientWidth = 448 + ClientWidth = 400 TabOrder = 1 inherited CloseBtn: TButton - Left = 393 + Left = 345 Top = 471 TabOrder = 14 end inherited ComputeBtn: TButton - Left = 309 + Left = 261 Top = 471 TabOrder = 13 end inherited ResetBtn: TButton - Left = 247 + Left = 199 Top = 471 TabOrder = 12 end inherited HelpBtn: TButton Tag = 104 - Left = 188 + Left = 140 Top = 471 TabOrder = 11 end inherited ButtonBevel: TBevel Top = 455 - Width = 448 + Width = 400 end object GroupBox1: TGroupBox[5] AnchorSideLeft.Control = ParamsPanel @@ -90,20 +90,19 @@ inherited AutoCorrForm: TAutoCorrForm Left = 142 Height = 68 Top = 0 - Width = 306 - Anchors = [akTop, akLeft, akRight] + Width = 252 AutoSize = True BorderSpacing.Left = 16 Caption = 'Include Cases:' ClientHeight = 48 - ClientWidth = 302 + ClientWidth = 248 TabOrder = 1 object Label1: TLabel AnchorSideLeft.Control = FromCaseEdit AnchorSideLeft.Side = asrBottom AnchorSideTop.Control = OnlyCasesBtn AnchorSideTop.Side = asrCenter - Left = 189 + Left = 178 Height = 15 Top = 23 Width = 12 @@ -144,7 +143,7 @@ inherited AutoCorrForm: TAutoCorrForm Left = 128 Height = 23 Top = 19 - Width = 57 + Width = 46 Alignment = taRightJustify BorderSpacing.Left = 4 TabOrder = 2 @@ -155,10 +154,10 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideLeft.Side = asrBottom AnchorSideTop.Control = OnlyCasesBtn AnchorSideTop.Side = asrCenter - Left = 205 + Left = 194 Height = 23 Top = 19 - Width = 58 + Width = 46 Alignment = taRightJustify BorderSpacing.Left = 4 BorderSpacing.Right = 8 @@ -185,9 +184,9 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideRight.Control = InBtn AnchorSideBottom.Control = GroupBox4 Left = 0 - Height = 169 + Height = 262 Top = 97 - Width = 205 + Width = 155 Anchors = [akTop, akLeft, akRight, akBottom] BorderSpacing.Top = 2 BorderSpacing.Right = 6 @@ -198,13 +197,15 @@ inherited AutoCorrForm: TAutoCorrForm TabOrder = 2 end object InBtn: TBitBtn[9] - AnchorSideLeft.Control = ParamsPanel AnchorSideLeft.Side = asrCenter AnchorSideTop.Control = VarList - Left = 211 + AnchorSideRight.Control = DepVarEdit + Left = 161 Height = 26 Top = 97 Width = 26 + Anchors = [akTop, akRight] + BorderSpacing.Right = 6 Images = MainDataModule.ImageList ImageIndex = 1 OnClick = InBtnClick @@ -212,15 +213,17 @@ inherited AutoCorrForm: TAutoCorrForm TabOrder = 3 end object OutBtn: TBitBtn[10] - AnchorSideLeft.Control = ParamsPanel AnchorSideLeft.Side = asrCenter AnchorSideTop.Control = InBtn AnchorSideTop.Side = asrBottom - Left = 211 + AnchorSideRight.Control = DepVarEdit + Left = 161 Height = 26 Top = 127 Width = 26 + Anchors = [akTop, akRight] BorderSpacing.Top = 4 + BorderSpacing.Right = 6 Images = MainDataModule.ImageList ImageIndex = 0 OnClick = OutBtnClick @@ -230,7 +233,7 @@ inherited AutoCorrForm: TAutoCorrForm object Label4: TLabel[11] AnchorSideLeft.Control = DepVarEdit AnchorSideBottom.Control = DepVarEdit - Left = 243 + Left = 193 Height = 15 Top = 101 Width = 88 @@ -240,18 +243,16 @@ inherited AutoCorrForm: TAutoCorrForm ParentColor = False end object DepVarEdit: TEdit[12] - AnchorSideLeft.Control = InBtn - AnchorSideLeft.Side = asrBottom + AnchorSideLeft.Control = GroupBox5 AnchorSideRight.Control = ParamsPanel AnchorSideRight.Side = asrBottom AnchorSideBottom.Control = OutBtn AnchorSideBottom.Side = asrBottom - Left = 243 + Left = 193 Height = 23 Top = 118 - Width = 205 + Width = 207 Anchors = [akLeft, akRight, akBottom] - BorderSpacing.Left = 6 BorderSpacing.Bottom = 12 ReadOnly = True TabOrder = 5 @@ -261,7 +262,7 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideTop.Control = AlphaEdit AnchorSideTop.Side = asrCenter AnchorSideRight.Control = AlphaEdit - Left = 312 + Left = 264 Height = 15 Top = 161 Width = 67 @@ -274,7 +275,7 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideTop.Control = MaxLagEdit AnchorSideTop.Side = asrCenter AnchorSideRight.Control = MaxLagEdit - Left = 296 + Left = 248 Height = 15 Top = 188 Width = 83 @@ -288,7 +289,7 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideTop.Side = asrBottom AnchorSideRight.Control = ParamsPanel AnchorSideRight.Side = asrBottom - Left = 387 + Left = 339 Height = 23 Top = 157 Width = 61 @@ -303,7 +304,7 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideTop.Side = asrBottom AnchorSideRight.Control = ParamsPanel AnchorSideRight.Side = asrBottom - Left = 387 + Left = 339 Height = 23 Top = 184 Width = 61 @@ -319,23 +320,23 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideTop.Side = asrBottom AnchorSideRight.Control = ParamsPanel AnchorSideRight.Side = asrBottom - Left = 221 + Left = 193 Height = 51 Top = 223 - Width = 227 + Width = 207 Anchors = [akTop, akLeft, akRight] AutoSize = True BorderSpacing.Top = 16 Caption = 'Projection Option' ClientHeight = 31 - ClientWidth = 223 + ClientWidth = 203 TabOrder = 8 object Label2: TLabel AnchorSideLeft.Control = ProjPtsEdit AnchorSideLeft.Side = asrBottom AnchorSideTop.Control = ProjPtsEdit AnchorSideTop.Side = asrCenter - Left = 135 + Left = 131 Height = 15 Top = 6 Width = 33 @@ -361,12 +362,12 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideLeft.Control = ProjectChk AnchorSideLeft.Side = asrBottom AnchorSideTop.Control = GroupBox3 - Left = 85 + Left = 81 Height = 23 Top = 2 Width = 42 Alignment = taRightJustify - BorderSpacing.Left = 16 + BorderSpacing.Left = 12 BorderSpacing.Top = 2 BorderSpacing.Bottom = 6 TabOrder = 1 @@ -378,75 +379,38 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideRight.Side = asrBottom AnchorSideBottom.Control = ButtonBevel Left = 0 - Height = 177 - Top = 278 - Width = 205 - Anchors = [akLeft, akRight, akBottom] + Height = 84 + Top = 371 + Width = 177 + Anchors = [akLeft, akBottom] AutoSize = True Caption = 'Data Smoothing:' ChildSizing.LeftRightSpacing = 12 ChildSizing.TopBottomSpacing = 6 ChildSizing.VerticalSpacing = 2 ChildSizing.Layout = cclLeftToRightThenTopToBottom - ClientHeight = 157 - ClientWidth = 201 + ChildSizing.ControlsPerLine = 1 + ClientHeight = 64 + ClientWidth = 173 TabOrder = 9 - object MeanChk: TCheckBox + object SmoothingParamsBtn: TButton Left = 12 - Height = 19 + Height = 25 Top = 6 - Width = 185 - Caption = 'Center on Mean' - TabOrder = 0 - end - object DifferenceChk: TCheckBox - Left = 12 - Height = 19 - Top = 27 - Width = 185 - Caption = 'Difference Smoothing' + Width = 149 + AutoSize = True + Caption = 'Parameters...' + OnClick = SmoothingParamsBtnClick TabOrder = 1 end - object MoveAvgChk: TCheckBox + object AutoRegSmoothChk: TCheckBox Left = 12 Height = 19 - Top = 48 - Width = 185 - Caption = 'Moving Average Smooth' - TabOrder = 2 - end - object ExpSmoothChk: TCheckBox - Left = 12 - Height = 19 - Top = 69 - Width = 185 - Caption = 'Exponentially Smooth' - TabOrder = 3 - end - object FourierSmoothChk: TCheckBox - Left = 12 - Height = 19 - Top = 90 - Width = 185 - Caption = 'Fourier Filter Smooth' - TabOrder = 4 - end - object PolyChk: TCheckBox - Left = 12 - Height = 19 - Top = 111 - Width = 185 - Caption = 'Polynomial Regression Smooth' - TabOrder = 5 - end - object MRegSmoothChk: TCheckBox - Left = 12 - Height = 19 - Top = 132 - Width = 185 - Caption = 'Multiple Regression Smooth' - OnChange = MRegSmoothChkChange - TabOrder = 6 + Top = 39 + Width = 149 + BorderSpacing.Top = 8 + Caption = 'Auto-regression smooth' + TabOrder = 0 end end object GroupBox5: TGroupBox[19] @@ -457,19 +421,21 @@ inherited AutoCorrForm: TAutoCorrForm AnchorSideRight.Side = asrBottom AnchorSideBottom.Control = GroupBox4 AnchorSideBottom.Side = asrBottom - Left = 221 - Height = 177 - Top = 278 - Width = 227 - Anchors = [akTop, akLeft, akRight, akBottom] + Left = 193 + Height = 135 + Top = 320 + Width = 207 + Anchors = [akLeft, akRight, akBottom] + AutoSize = True BorderSpacing.Left = 16 + BorderSpacing.Top = 16 Caption = 'Analysis / Output Options' ChildSizing.LeftRightSpacing = 12 ChildSizing.TopBottomSpacing = 6 ChildSizing.VerticalSpacing = 2 ChildSizing.Layout = cclLeftToRightThenTopToBottom - ClientHeight = 157 - ClientWidth = 223 + ClientHeight = 115 + ClientWidth = 203 TabOrder = 10 object PlotChk: TCheckBox Left = 12 @@ -477,6 +443,8 @@ inherited AutoCorrForm: TAutoCorrForm Top = 6 Width = 148 Caption = 'Plot correlogram' + Checked = True + State = cbChecked TabOrder = 0 end object StatsChk: TCheckBox @@ -511,28 +479,19 @@ inherited AutoCorrForm: TAutoCorrForm Caption = 'Yule-Walker coefficients' TabOrder = 4 end - object ResidChk: TCheckBox - Left = 12 - Height = 19 - Top = 111 - Width = 148 - Caption = 'Residual plot' - Enabled = False - TabOrder = 5 - end end end inherited ParamsSplitter: TSplitter - Left = 460 + Left = 412 Height = 512 end inherited PageControl: TPageControl - Left = 469 + Left = 421 Height = 496 - Width = 514 - ActivePage = YuleWalkerPage + Width = 562 + ActivePage = AutoRegPage MultiLine = True - TabIndex = 4 + TabIndex = 5 TabOrder = 3 Options = [nboMultiLine] inherited ChartPage: TTabSheet @@ -540,7 +499,7 @@ inherited AutoCorrForm: TAutoCorrForm TabVisible = False end object CorrelationMatrixPage: TTabSheet[2] - Caption = 'Correlation matrix' + Caption = 'Correlation Matrix' TabVisible = False end object PartialsPage: TTabSheet[3] @@ -551,6 +510,10 @@ inherited AutoCorrForm: TAutoCorrForm Caption = 'Yule-Walker' TabVisible = False end + object AutoRegPage: TTabSheet[5] + Caption = 'Auto-Regr. Smooth' + TabVisible = False + end end object Panel2: TPanel[3] AnchorSideRight.Control = Owner diff --git a/applications/lazstats/source/forms/analysis/correlation/autocorunit.pas b/applications/lazstats/source/forms/analysis/correlation/autocorunit.pas index e28b38044..4a51ba4a1 100644 --- a/applications/lazstats/source/forms/analysis/correlation/autocorunit.pas +++ b/applications/lazstats/source/forms/analysis/correlation/autocorunit.pas @@ -2,17 +2,25 @@ // - original file not found, but could be reconstructed from graphs. // --> there is no EXACT agreement of numbers with the pdf file "autocorrelation.pdf". +// Code looks as if conversion from C++ has not yet been completed. +// Fourier smooth differs from results obtained by OpenStat. But also the +// unrefactored Lazstats differs from OpenStat. +// +// Possible copyright issue: +// Fourier smooth uses code from Numerical Recipes + unit AutoCorUnit; {$mode objfpc}{$H+} +//{$modeswitch advancedrecords} interface uses Classes, SysUtils, Forms, Controls, Graphics, Dialogs, StdCtrls, ExtCtrls, - Buttons, ComCtrls, MainUnit, FunctionsLib, Globals, GraphLib, DataProcs, - MatrixLib, PointsUnit, DifferenceUnit, ReportFrameUnit, - BasicStatsReportAndChartFormUnit; + Buttons, ComCtrls, MainUnit, FunctionsLib, Globals, //GraphLib, + MatrixLib, ReportFrameUnit, ChartFrameUnit, + BasicStatsReportAndChartFormUnit, SmoothingUnit; type @@ -20,6 +28,7 @@ type TAutoCorrForm = class(TBasicStatsReportAndChartForm) AlphaEdit: TEdit; + SmoothingParamsBtn: TButton; Panel2: TPanel; PlotChk: TCheckBox; StatsChk: TCheckBox; @@ -27,9 +36,9 @@ type PartialsChk: TCheckBox; PartialsPage: TTabSheet; CorrelationMatrixPage: TTabSheet; + AutoRegPage: TTabSheet; YuleWalkerPage: TTabSheet; YuleWalkerChk: TCheckBox; - ResidChk: TCheckBox; GroupBox5: TGroupBox; MaxLagEdit: TEdit; InBtn: TBitBtn; @@ -40,13 +49,7 @@ type Label3: TLabel; Label4: TLabel; VarList: TListBox; - MRegSmoothChk: TCheckBox; - PolyChk: TCheckBox; - FourierSmoothChk: TCheckBox; - ExpSmoothChk: TCheckBox; - MoveAvgChk: TCheckBox; - DifferenceChk: TCheckBox; - MeanChk: TCheckBox; + AutoRegSmoothChk: TCheckBox; GroupBox4: TGroupBox; Label2: TLabel; ProjPtsEdit: TEdit; @@ -61,9 +64,9 @@ type AllCasesBtn: TRadioButton; OnlyCasesBtn: TRadioButton; RowBtn: TRadioButton; + procedure SmoothingParamsBtnClick(Sender: TObject); procedure ColBtnClick(Sender: TObject); procedure InBtnClick(Sender: TObject); - procedure MRegSmoothChkChange(Sender: TObject); procedure OutBtnClick(Sender: TObject); procedure RowBtnClick(Sender: TObject); procedure VarListDblClick(Sender: TObject); @@ -75,24 +78,15 @@ type FCorrelationMatrixReportFrame: TReportFrame; FPartialsReportFrame: TReportFrame; FYuleWalkerReportFrame: TReportFrame; + FAutoRegChartFrame: TChartFrame; + FParams: TSmoothingParams; - FExpSmoothAlpha: Double; - FMovAvgRawWeights: DblDyneVec; - FMovAvgOrder: Integer; - FPolyOrder: Integer; procedure AutoPlot(const ACorrs, APartialCorrs: DblDyneVec; AUpperLimit, ALowerLimit: Double); + procedure AutoRegPlot(const APoints, ASmoothedPoints: DblDyneVec); function CalcMean(const pts: DblDyneVec; NoPts: Integer): Double; - procedure ExponentialSmooth(var Pts: DblDyneVec; NoPts: Integer); - procedure Four1(var data: DblDyneVec; nn: LongWord; isign: integer); - procedure Fourier(var data: DblDyneVec; n, npts: integer); - procedure FourierSmooth(var Pts: DblDyneVec; NoPts: Integer); procedure GetPts(var Pts: DblDyneVec; var ANumPts: Integer; DepVar: Integer); - procedure MovingAverage(var Pts: DblDyneVec; NoPts: Integer); - procedure PlotDifferencesForLag(var Pts: DblDyneVec; NoPts: Integer); - procedure PolyFit(const pts: DblDyneVec; var avg: DblDyneVec; NoPts, Order: integer); - procedure PolynomialSmooth(var Pts: DblDyneVec; NoPts: Integer); - procedure RealFT(var data: DblDyneVec; n: LongWord; isign: integer); + function Prepare(out ADepVar: Integer; out APts: DblDyneVec): Boolean; procedure RemoveMeans(var Pts: DblDyneVec; NoPts: Integer; AMean: Double); protected @@ -109,6 +103,7 @@ type var AutoCorrForm: TAutoCorrForm; + implementation {$R *.lfm} @@ -116,8 +111,8 @@ implementation uses Math, TALegend, - GridProcs, MathUnit, MatrixUnit, ChartFrameUnit, - MoveAvgUnit, AutoPlotUnit, PolynomialUnit, ExpSmoothUnit, FFTUnit; + Utils, GridProcs, MathUnit, MatrixUnit; + { TAutoCorrForm } @@ -139,6 +134,10 @@ begin FYuleWalkerReportFrame.Parent := YuleWalkerPage; FYuleWalkerReportFrame.Align := alClient; + FAutoRegChartFrame := TChartFrame.Create(self); + FAutoRegChartFrame.Parent := AutoRegPage; + FAutoRegChartFrame.Align := alClient; + FCreated := true; Reset; end; @@ -150,11 +149,10 @@ begin ParamsPanel.Constraints.MinWidth := MaxValue([ GroupBox2.Left + GroupBox2.Width, - 4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left{, - GroupBox5.Width + GroupBox4.BorderSpacing.Left + GroupBox4.Width} + 4*CloseBtn.Width + 3*CloseBtn.BorderSpacing.Left ]); ParamsPanel.Constraints.MinHeight := GroupBox3.Top + GroupBox3.Height + - VarList.BorderSpacing.Bottom + GroupBox4.Height + + GroupBox5.BorderSpacing.Top + GroupBox5.Height + ButtonBevel.Height + CloseBtn.BorderSpacing.Top + CloseBtn.Height; end; @@ -191,6 +189,36 @@ begin end; + procedure TAutoCorrForm.AutoRegPlot(const APoints, ASmoothedPoints: DblDyneVec); + var + x: DblDyneVec; + i: Integer; + begin + FAutoRegChartFrame.SetTitle('Autoregression Smoothing'); + + SetLength(x, Length(APoints)); + for i := 0 to Length(APoints)-1 do x[i] := i; + + FAutoRegChartFrame.PlotXY(ptLinesAndSymbols, x, APoints, nil, nil, 'Original', DATA_COLORS[0]); + FAutoRegChartFrame.PlotXY(ptLinesAndSymbols, x, ASmoothedPoints, nil, nil, 'Smoothed', DATA_COLORS[1]); + + FAutoRegChartFrame.Chart.Legend.Alignment := laBottomCenter; + FAutoRegChartFrame.Chart.Legend.ColumnCount := 3; + (* + + avg[0] := pts[0]; + PointsFrm.pts := pts; + PointsFrm.avg := avg; + PointsFrm.NoCases := NoPts + noproj; + PointsFrm.LabelOne := 'Original'; + PointsFrm.LabelTwo := 'Smoothed'; + PointsFrm.Title := 'Autoregressive Smoothed'; + PointsFrm.Caption := 'Autoregression Smoothing'; + PointsFrm.ShowModal; + *) + end; + + procedure TAutoCorrForm.Reset; var i: integer; @@ -200,14 +228,18 @@ begin if not FCreated then exit; + FParams.CopyFrom(DEFAULT_SMOOTHING_PARAMS); + ReportPage.Caption := 'Reports'; ChartPage.TabVisible := false; CorrelationMatrixPage.TabVisible := false; PartialsPage.TabVisible := false; YuleWalkerPage.TabVisible := false; + AutoRegPage.TabVisible := false; FCorrelationMatrixReportFrame.Clear; FPartialsReportFrame.Clear; FYuleWalkerReportFrame.Clear; + FAutoRegChartFrame.Clear; VarList.Clear; if ColBtn.Checked then @@ -223,19 +255,12 @@ begin DepVarEdit.Text := ''; MaxLagEdit.Text := '30'; + PlotChk.Checked := true; StatsChk.Checked := false; RmatChk.Checked := false; PartialsChk.Checked := false; - PlotChk.Checked := false; - ResidChk.Checked := false; - DifferenceChk.Checked := false; - PolyChk.Checked := false; - MeanChk.Checked := false; - MoveAvgChk.Checked := false; - ExpSmoothChk.Checked := false; - FourierSmoothChk.Checked := false; YuleWalkerChk.Checked := false; - ResidChk.Enabled := MRegSmoothChk.Checked; + AutoRegSmoothChk.Checked := false; FromCaseEdit.Text := ''; ToCaseEdit.Text := ''; @@ -243,11 +268,6 @@ begin AlphaEdit.Text := FormatFloat('0.00', DEFAULT_ALPHA_LEVEL); ProjPtsEdit.Text := ''; - FMovAvgOrder := 1; - FMovAvgRawWeights := nil; - FExpSmoothAlpha := 0.99; - FPolyOrder := 1; - UpdateBtnStates; end; @@ -262,11 +282,368 @@ begin DepVarEdit.Clear; GroupBox2.Caption := 'Include Columns:'; AllCasesBtn.Caption := 'All Variables'; - OnlyCasesBtn.Caption := 'Only Columns From:'; + OnlyCasesBtn.Caption := 'Only Cols From:'; UpdateBtnStates; end; +procedure TAutoCorrForm.Compute; +var + ColNoSelected: IntDyneVec = nil; + RowLabels: StrDyneVec = nil; + ColLabels: StrDyneVec = nil; + pts: DblDyneVec = nil; + avg: DblDyneVec = nil; + betas: DblDyneVec = nil; + residual: DblDyneVec = nil; + rxy: DblDyneVec = nil; + a: DblDyneMat = nil; + correlations: DblDyneMat = nil; + PartCors: DblDyneVec = nil; + noPts: Integer = 0; + DepVar, noSelected, noProj, nPoints, count, lag, maxLag, nCors: Integer; + i, j, k: Integer; + X, Y, mx, my, vx, vy, sx, sy, r, sampTrans, stdErr, zConf, UCL, LCL: Double; + mean, covZero, confidence, constant, uplimit, lowlimit, YHat, varResid: Double; + title: String; + lReport: TStrings; + +begin + if not Prepare(DepVar, pts) then + exit; + + noPts := Length(pts); + noSelected := 1; + SetLength(ColNoSelected, noSelected); + ColNoSelected[0] := DepVar; + + if ProjectChk.Checked then + noProj := StrToInt(ProjPtsEdit.Text) + else + noProj := 0; + + // Get the maximum lag value + maxLag := StrToInt(MaxLagEdit.Text); + if maxlag > NoPts div 2 then + maxlag := NoPts div 2; + if StrToInt(MaxLagEdit.Text) > maxlag then + MaxLagEdit.Text := IntToStr(maxlag); + nPoints := maxlag + 2; + + // Remove mean from all observations + mean := CalcMean(pts, noPts); + if FParams.CenterOnMean then + RemoveMeans(pts, noPts, mean); + + // Execute smoothing operations + if FParams.DifferenceParams.Active then + DifferenceSmooth(pts, FParams.DifferenceParams.Order, FParams.DifferenceParams.MaxLag); + if FParams.MovingAvgParams.Active then + MovingAvgSmooth(pts, FParams.MovingAvgParams.Weights, FParams.MovingAvgParams.Order, noProj); + if FParams.ExponentialParams.Active then + ExponentialSmooth(pts, FParams.ExponentialParams.Alpha, noProj); + if FParams.FourierParams.Active then + FourierSmooth(pts, FParams.FourierParams.Strength, noProj); + if FParams.PolynomialParams.Active then + PolynomialSmooth(pts, FParams.PolynomialParams.Order); + + // Get mean and variance of (transformed) points + mean := 0.0; + covZero := 0.0; + for i := 0 to NoPts - 1 do + mean := mean + pts[i]; + mean := mean / NoPts; + for i := 0 to NoPts-1 do + if FParams.CenterOnMean then + covZero := covZero + sqr(pts[i]) + else + covZero := covZero + sqr(pts[i] - mean); + covZero := covZero / NoPts; + + confidence := StrToFloat(AlphaEdit.Text); + nCors := 0; + + if maxlag > NoPts - 3 then + begin + maxlag := NoPts - 3; + MaxLagEdit.Text := IntToStr(maxlag); + end; + + // Allocate space for covariance and correlation matrices, etc. + SetLength(correlations, nPoints+1, nPoints+1); + SetLength(PartCors, npoints); + SetLength(a, npoints, npoints); + SetLength(betas, npoints); + SetLength(rxy, npoints); + + SetLength(avg, noPts+noproj); // wp: was +10 + SetLength(residual, noPts+noproj); // wp: was +10 + + // Prepare label arrays + SetLength(RowLabels, npoints); + SetLength(ColLabels, npoints); + for i := 0 to npoints-1 do + begin + RowLabels[i] := 'Lag ' + IntToStr(i); + ColLabels[i] := RowLabels[i]; + end; + + lReport := TStringList.Create; + try + if StatsChk.Checked then + begin + lReport.Add('Overall mean: %8.3f', [mean]); + lReport.Add('Variance: %8.3f', [covzero]); + lReport.Add(''); + + lReport.Add(' Lag Rxy MeanX MeanY Std.Dev.X Std.Dev.Y Cases LCL UCL '); + lReport.Add('----- ------------ ------------ ------------ ------------ ------------ ------- ------------ ------------'); + end; + + if covzero = 0 then + lReport.Add('ERROR: zero variance detected.') + else + begin + correlations[0, 0] := 1.0; + + for lag := 0 to maxlag do + begin + r := 0.0; + vx := 0.0; + vy := 0.0; + mx := 0.0; + my := 0.0; + count := 0; + for i := 1 to (NoPts - lag) do + begin + X := pts[i-1]; + Y := pts[i-1+lag]; + if FParams.CenterOnMean then + r := r + X * Y + else + r := r + (X - mean) * (Y - mean); + vx := vx + sqr(X); + vy := vy + sqr(Y); + mx := mx + X; + my := my + Y; + count := count + 1; + end; + r := r / NoPts; + vx := vx - sqr(mx) / count; + vx := vx / (count - 1); + sx := sqrt(vx); + + vy := vy - sqr(my) / Count; + vy := vy / (count - 1); + sy := sqrt(vy); + + mx := mx / count; + my := my / count; + r := r / covzero; + if (abs(r) < 1.0) then + sampTrans := ln((1.0 + r) / (1.0 - r)) / 2.0; + stdErr := sqrt(1.0 / (NoPts - 3.0)); + zconf := abs(InverseZ(confidence * 0.5)); + if (abs(r) < 1.0) then + begin + //z := samptrans / StdErr; + UCL := sampTrans + zconf * StdErr; + LCL := sampTrans - zconf * StdErr; + UCL := (exp(2.0 * UCL) - 1.0) / (exp(2.0 * UCL) + 1.0); + LCL := (exp(2.0 * LCL) - 1.0) / (exp(2.0 * LCL) + 1.0); + end else + begin + UCL := 1.0; + LCL := 1.0; + end; + if StatsChk.Checked then + lReport.Add('%5d %12.4f %12.4f %12.4f %12.4f %12.4f %7d %12.4f %12.4f', + [lag, r, mx, my, sx, sy, count, LCL, UCL] + ); + ncors := ncors + 1; + correlations[0, lag] := r; + correlations[lag, 0] := r; + end; // next lag + end; + + if StatsChk.Checked then + begin + ReportPage.caption := 'Statistics'; + FStatsReportFrame.DisplayReport(lReport); + lReport.Clear; + end; + ReportPage.TabVisible := StatsChk.Checked; + + // build correlation matrix + for i := 0 to maxlag do + correlations[i, i] := 1.0; + for i := 1 to maxlag do + for j := i+1 to maxlag do + begin + correlations[i, j] := correlations[0, j-i]; + correlations[j, i] := correlations[i, j]; + end; + + // Print the correlation matrix if elected + if RmatChk.Checked then + begin + title := 'Matrix of Lagged Variable: ' + DepVarEdit.Text; + MatPrint(correlations, maxLag+1, maxLag+1, title, RowLabels, ColLabels, NoPts, lReport); + FCorrelationMatrixReportFrame.DisplayReport(lReport); + lReport.Clear; + end; + CorrelationMatrixPage.TabVisible := RmatChk.Checked; + + // Calculate partial correlations + PartCors[0] := 1.0; + for i := 1 to maxlag do // start at lag 1 + begin + for j := 1 to i do + begin + for k := 1 to i do + a[j-1, k-1] := correlations[j, k]; + rxy[j-1] := correlations[0,j]; + end; + + SVDinverse(a, i); + + // get betas as product of inverse times vector + for j := 1 to i do + begin + betas[j-1] := 0.0; + for k := 1 to i do + betas[j-1] := betas[j-1] + (a[j-1,k-1] * rxy[k-1]); + end; + + // get regression constant + // Note - since variance of Y and each X is the same, B = beta for an X + Constant := 0; + if not FParams.CenterOnMean then + begin + for j := 1 to i do + Constant := Constant + betas[j-1] * Mean; + Constant := Mean - Constant; + end; + + // calculate predicted value and residual + // Note - the dependent variable predicted is the next value in the + // time series using each of the previous time period values + YHat := 0.0; + for j := 0 to i-1 do + YHat := YHat + (betas[j] * pts[j]); + YHat := YHat + Constant; + avg[i] := YHat; + residual[i] := pts[i] - YHat; + + // print betas if elected + if YuleWalkerChk.Checked then + begin + if lReport.Count > 0 then + begin + lReport.Add(''); + if i = 1 then lReport.Add(DIVIDER) else lReport.Add(DIVIDER_SMALL); + lReport.Add(''); + end; + title := 'Yule-Walker Coefficients for lag ' + IntToStr(i); + DynVectorPrint(betas, i, title, ColLabels, NoPts, lReport); + lReport.Add('Constant: %10.3f', [Constant]); + end; + + PartCors[i] := betas[i-1]; + end; // next i (lag from 1 to maxlag) + + // Print Yule-Walker coefficients + if YuleWalkerChk.Checked then + begin + FYuleWalkerReportFrame.DisplayReport(lReport); + lReport.Clear; + end; + YuleWalkerPage.TabVisible := YuleWalkerChk.Checked; + + // Print partial correlations if elected + if PartialsChk.Checked then + begin + title := 'Partial Correlation Coefficients'; + DynVectorPrint(PartCors, maxlag, title, ColLabels, NoPts, lReport); + FPartialsReportFrame.DisplayReport(lReport); + lReport.Clear; + end; + PartialsPage.TabVisible := PartialsChk.Checked; + + // plot correlations if elected + if PlotChk.Checked then + begin + uplimit := 1.96 * (1.0 / sqrt(count)); + lowlimit := -1.96 * (1.0 / sqrt(count)); + + SetLength(rXY, maxLag+1); + for i := 0 to maxlag do + rxy[i] := correlations[0, i]; + + AutoPlot(rxy, partCors, upLimit, lowLimit); + end; + ChartPage.TabVisible := PlotChk.Checked; + + if AutoRegSmoothChk.Checked then + begin + // calculate predicted values and residuals for remaining points + // Note - the dependent variable predicted is the next value in + // the time series using each of the previous time period values + // as predictors + for i := maxlag to (NoPts + noproj - 1) do + begin + Yhat := 0.0; + for j := 0 to maxlag do + Yhat := Yhat + (betas[j] * pts[i-maxlag+j]); + Yhat := Yhat + Constant; + avg[i] := Yhat; + residual[i] := pts[i] - Yhat; + end; + + AutoRegPlot(pts, avg); + (* + // plot points smoothed by autoregression + avg[0] := pts[0]; + PointsFrm.pts := pts; + PointsFrm.avg := avg; + PointsFrm.NoCases := NoPts + noproj; + PointsFrm.LabelOne := 'Original'; + PointsFrm.LabelTwo := 'Smoothed'; + PointsFrm.Title := 'Autoregressive Smoothed'; + PointsFrm.Caption := 'Autoregression Smoothing'; + PointsFrm.ShowModal; *) + + // plot residuals if elected + (* + if ResidChk.Checked then + begin + varResid := 0.0; + residual[0] := 0.0; + for i := 1 to maxlag do + varResid := varResid + sqr(residual[i]); + varResid := varResid / maxlag; + StdErr := sqrt(varResid); + + // plot the residuals + PointsFrm.pts := pts; + PointsFrm.avg := residual; + PointsFrm.NoCases := NoPts; + PointsFrm.MsgEdit.Text := 'Std. Err. Residuals: ' + FloatToStr(StdErr); + PointsFrm.LabelOne := 'Original'; + PointsFrm.LabelTwo := 'Residuals'; + PointsFrm.Title := 'Residuals from Autoregression Smoothing'; + PointsFrm.Caption := 'Autoregressive Residuals'; + PointsFrm.ShowModal; + end; + *) + end; + AutoRegPage.TabVisible := AutoRegSmoothChk.Checked; + + finally + lReport.Free; + end; +end; + (* procedure TAutoCorrForm.Compute; var X, Y, count, covzero, mean: double; @@ -294,9 +671,14 @@ var confidence: double; lReport: TStrings; begin - NoSelected := 1; - SetLength(ColNoSelected, NoSelected); + if not Prepare(DepVar, pts) then + exit; + noPts := Length(pts); + noSelected := 1; + SetLength(ColNoSelected, noSelected); + ColNoSelected[0] := DepVar; + { if ColBtn.Checked then begin // Get column of the selected variable @@ -323,12 +705,12 @@ begin end; end; ColNoSelected[0] := DepVar; - + } if ProjectChk.Checked then noProj := StrToInt(ProjPtsEdit.Text) else noProj := 0; - + { // Get points to analyze GetPts(pts, noPts, DepVar); if noPts = 0 then @@ -336,6 +718,7 @@ begin MessageDlg('No data points found.', mtError, [mbOk], 0); exit; end; + } count := noPts; SetLength(pts, noPts+noproj+10); // why +10? @@ -368,55 +751,58 @@ begin ColLabels[i] := RowLabels[i]; end; + // Calculate mean of all values + mean := CalcMean(pts, noPts); + + // Remove mean from all observations if elected + if MeanChk.Checked then + RemoveMeans(pts, noPts, mean); + + correlations[0, 0] := 1.0; + + // Get differences for lag specified + if DifferenceSmoothChk.Checked then + PlotDifferencesForLag(Pts, NoPts); + + // Get moving average if checked + if MovingAvgSmoothChk.Checked then + MovingAverage(Pts, NoPts); + + // Do exponential smoothing if requested + if ExpSmoothChk.Checked then + ExponentialSmooth(pts, NoPts); + + // Fast Fourier smoothing, if requested + if FourierSmoothChk.Checked then + FourierSmooth(pts, NoPts); + + // Get polynomial regression fit if elected + if PolySmoothChk.Checked then + PolynomialSmooth(pts, NoPts); + + // get mean and variance of (transformed?) points + mean := 0.0; + covZero := 0.0; + for i := 0 to NoPts - 1 do + mean := mean + pts[i]; + mean := mean / NoPts; + for i := 0 to NoPts-1 do + if MeanChk.Checked then + covZero := covZero + sqr(pts[i]) + else + covZero := covZero + sqr(pts[i] - mean); + covZero := covZero / NoPts; + + confidence := StrToFloat(AlphaEdit.Text); + ncors := 0; + if maxlag > NoPts - 3 then + begin + maxlag := NoPts - 3; + MaxLagEdit.Text := IntToStr(maxlag); + end; + lReport := TStringList.Create; try - - uplimit := 0.0; - lowlimit := 0.0; - covzero := 0.0; - - // Calculate mean of all values - mean := CalcMean(pts, noPts); - - // Remove mean from all observations if elected - if MeanChk.Checked then - RemoveMeans(pts, noPts, mean); - - correlations[0, 0] := 1.0; - - // Get differences for lag specified - if DifferenceChk.Checked then - PlotDifferencesForLag(Pts, NoPts); - - // Get moving average if checked - if MoveAvgChk.Checked then - MovingAverage(Pts, NoPts); - - // Do exponential smoothing if requested - if ExpSmoothChk.Checked then - ExponentialSmooth(pts, NoPts); - - // Fast Fourier smoothing, if requested - if FourierSmoothChk.Checked then - FourierSmooth(pts, NoPts); - - // Get polynomial regression fit if elected - if PolyChk.Checked then - PolynomialSmooth(pts, NoPts); - - // get mean and variance of (transformed?) points - mean := 0.0; - for i := 0 to NoPts - 1 do - mean := mean + pts[i]; - mean := mean / NoPts; - for i := 0 to NoPts-1 do - if MeanChk.Checked then - covzero := covzero + sqr(pts[i]) - else - covzero := covzero + sqr(pts[i] - mean); - covzero := covzero / NoPts; - - // get correlations for each lag 0 to maxlag if StatsChk.Checked then begin lReport.Add('Overall mean: %8.3f', [mean]); @@ -427,74 +813,71 @@ begin lReport.Add('----- ------------ ------------ ------------ ------------ ------------ ------- ------------ ------------'); end; - confidence := StrToFloat(AlphaEdit.Text); - ncors := 0; - if maxlag > NoPts - 3 then + if covzero = 0 then + lReport.Add('ERROR: zero variance detected.') + else begin - maxlag := NoPts - 3; - maxlagedit.Text := IntToStr(maxlag); + for lag := 0 to maxlag do + begin + r := 0.0; + vx := 0.0; + vy := 0.0; + mx := 0.0; + my := 0.0; + Count := 0.0; + // lagvalue[lag] := lag; + for i := 1 to (NoPts - lag) do + begin + X := pts[i-1]; + Y := pts[i-1+lag]; + if MeanChk.Checked then + r := r + X * Y + else + r := r + (X - mean) * (Y - mean); + vx := vx + sqr(X); + vy := vy + sqr(Y); + mx := mx + X; + my := my + Y; + Count := Count + 1.0; + end; + r := r / NoPts; + vx := vx - sqr(mx) / Count; + vx := vx / (Count - 1.0); + sx := sqrt(vx); + + vy := vy - sqr(my) / Count; + vy := vy / (Count - 1.0); + sy := sqrt(vy); + + mx := mx / Count; + my := my / Count; + r := r / covzero; + if (abs(r) < 1.0) then + samptrans := ln((1.0 + r) / (1.0 - r)) / 2.0; + stdErr := sqrt(1.0 / (NoPts - 3.0)); + zconf := abs(InverseZ(confidence * 0.5)); + if (abs(r) < 1.0) then + begin + //z := samptrans / StdErr; + UCL := samptrans + zconf * StdErr; + LCL := samptrans - zconf * StdErr; + UCL := (exp(2.0 * UCL) - 1.0) / (exp(2.0 * UCL) + 1.0); + LCL := (exp(2.0 * LCL) - 1.0) / (exp(2.0 * LCL) + 1.0); + end else + begin + UCL := 1.0; + LCL := 1.0; + end; + if StatsChk.Checked then + lReport.Add('%5d %12.4f %12.4f %12.4f %12.4f %12.4f %7.0f %12.4f %12.4f', + [lag, r, mx, my, sx, sy, Count, LCL, UCL] + ); + ncors := ncors + 1; + correlations[0,lag] := r; + correlations[lag,0] := r; + end; // next lag end; - for lag := 0 to maxlag do - begin - r := 0.0; - vx := 0.0; - vy := 0.0; - mx := 0.0; - my := 0.0; - Count := 0.0; -// lagvalue[lag] := lag; - for i := 1 to (NoPts - lag) do - begin - X := pts[i-1]; - Y := pts[i-1+lag]; - if MeanChk.Checked then - r := r + X * Y - else - r := r + (X - mean) * (Y - mean); - vx := vx + sqr(X); - vy := vy + sqr(Y); - mx := mx + X; - my := my + Y; - Count := Count + 1.0; - end; - r := r / NoPts; - vx := vx - sqr(mx) / Count; - vx := vx / (Count - 1.0); - sx := sqrt(vx); - - vy := vy - sqr(my) / Count; - vy := vy / (Count - 1.0); - sy := sqrt(vy); - - mx := mx / Count; - my := my / Count; - r := r / covzero; - if (abs(r) < 1.0) then - samptrans := ln((1.0 + r) / (1.0 - r)) / 2.0; - stdErr := sqrt(1.0 / (NoPts - 3.0)); - zconf := abs(InverseZ(confidence * 0.5)); - if (abs(r) < 1.0) then - begin - //z := samptrans / StdErr; - UCL := samptrans + zconf * StdErr; - LCL := samptrans - zconf * StdErr; - UCL := (exp(2.0 * UCL) - 1.0) / (exp(2.0 * UCL) + 1.0); - LCL := (exp(2.0 * LCL) - 1.0) / (exp(2.0 * LCL) + 1.0); - end else - begin - UCL := 1.0; - LCL := 1.0; - end; - if StatsChk.Checked then - lReport.Add('%5d %12.4f %12.4f %12.4f %12.4f %12.4f %7.0f %12.4f %12.4f', - [lag, r, mx, my, sx, sy, Count, LCL, UCL] - ); - ncors := ncors + 1; - correlations[0,lag] := r; - correlations[lag,0] := r; - end; // next lag - if StatsChk.Checked then begin ReportPage.caption := 'Statistics'; @@ -626,7 +1009,7 @@ begin end; ChartPage.TabVisible := PlotChk.Checked; - if MRegSmoothChk.Checked then + if AutoRegSmoothChk.Checked then begin // calculate predicted values and residuals for remaining points // Note - the dependent variable predicted is the next value in @@ -643,8 +1026,6 @@ begin end; // plot points smoothed by autoregression - if PointsFrm = nil then - Application.CreateForm(TPointsFrm, PointsFrm); avg[0] := pts[0]; PointsFrm.pts := pts; PointsFrm.avg := avg; @@ -682,7 +1063,7 @@ begin lReport.Free; end; end; - + *) procedure TAutoCorrForm.ColBtnClick(Sender: TObject); var i: integer; @@ -697,6 +1078,7 @@ begin UpdateBtnStates; end; + procedure TAutoCorrForm.InBtnClick(Sender: TObject); var index: integer; @@ -710,10 +1092,6 @@ begin end; end; -procedure TAutoCorrForm.MRegSmoothChkChange(Sender: TObject); -begin - ResidChk.Enabled := MRegSmoothChk.Checked; -end; procedure TAutoCorrForm.OutBtnClick(Sender: TObject); begin @@ -725,300 +1103,6 @@ begin end; end; -procedure TAutoCorrForm.Four1(var data: DblDyneVec; nn: longword; isign: integer); -var - n, mmax, m, j, istep, i: longword; - wtemp, wr, wpr, wpi, wi, theta: double; - tempr, tempi: double; -begin - n := 2 * nn; - j := 1; - i := 1; - while i < n do - begin - if (j > i) then - begin - tempr := data[j]; - tempi := data[j+1]; - data[j] := data[i]; - data[j+1] := data[i+1]; - data[i] := tempr; - data[i+1] := tempi; - end; - m := n div 2; - while (m >= 2) and (j > m) do - begin - j := j - m; - m := m div 2; - end; - j := j + m; - i := i + 2; - end; - mmax := 2; - while (n > mmax) do - begin - istep := 2 * mmax; - theta := isign * TWO_PI / mmax; - wtemp := sin(0.5 * theta); - wpr := -2.0 * wtemp * wtemp; - wpi := sin(theta); - wr := 1.0; - wi := 0.0; - m := 1; - while m < mmax do - begin - i := m; - while i <= n do - begin - j := i + mmax; - tempr := wr * data[j] - wi * data[j+1]; - tempi := wr * data[j+1] + wi * data[j]; - data[j] := data[i] - tempr; - data[j+1] := data[i+1] - tempi; - data[i] := data[i] + tempr; - data[i+1] := data[i+1] + tempi; - i := i + istep; - end; - wtemp := wr; - wr := wr * wpr - wi * wpi + wr; - wi := wi * wpr + wtemp * wpi + wi; - m := m + 2; - end; - mmax := istep; - end; -end; - -procedure TAutoCorrForm.RealFt(var data: DblDyneVec; n: longword; isign: integer); -var - i,i1,i2,i3,i4,np3 : integer; // was: longword; - c1,c2,h1r,h1i,h2r,h2i : double; - wr,wi,wpr,wpi,wtemp,theta : double; -begin - c1 := 0.5; - theta := PI / ( n div 2); - if (isign = 1) then - begin - c2 := -0.5; - Four1(data, n div 2, 1); - end else - begin - c2 := 0.5; - theta := -theta; - end; - wtemp := sin(0.5 * theta); - wpr := -2.0 * wtemp * wtemp; - wpi := sin(theta); - wr := 1.0 + wpr; - wi := wpi; - np3 := n + 3; - for i := 2 to n div 2 do - begin - i1 := i + i - 1; - i2 := 1 + i1; - i3 := np3 - i2; - i4 := 1 + i3; - h1r := c1 * (data[i1] + data[i3]); - h1i := c1 * (data[i2] - data[i4]); - h2r := -c2 * (data[i2] + data[i4]); - h2i := c2 * (data[i1] - data[i3]); - data[i1] := h1r + wr * h2r - wi * h2i; - data[i2] := h1i + wr * h2i + wi * h2r; - data[i3] := h1r - wr * h2r + wi * h2i; - data[i4] := -h1i + wr * h2i + wi * h2r; - wtemp := wr; - wr := wtemp * wpr - wi * wpi + wr; - wi := wi * wpr + wtemp * wpi + wi; - end; - - if (isign = 1) then - begin - h1r := data[1]; - data[1] := h1r + data[2]; - data[2] := h1r - data[2]; - end else - begin - h1r := data[1]; - data[1] := c1 * (h1r + data[2]); - data[2] := c1 * (h1r - data[2]); - Four1(data,n div 2,-1); - end; -end; - -procedure TAutoCorrForm.Fourier(var data: DblDyneVec; n, npts: integer); -var - nmin, m, mo2, i, k, j: integer; - yn, y1, rn1, fac, cnst: double; - y: DblDyneVec; -begin - m := 2; - nmin := n + 2*npts; - while (m < nmin) do m := m * 2; - - cnst := sqr(npts / m); - - SetLength(y, m+1); - for i := 0 to n - 1 do - y[i+1] := data[i]; - y1 := y[1]; - yn := y[n]; - rn1 := 1.0 / (n - 1); - for j := 1 to n do - y[j] := y[j] + (-rn1 * (y1 * (n - j) + y1 * (j - 1))); - for j := n+1 to m do - y[j] := 0.0; - mo2 := m div 2; - - RealFT(y, mo2, 1); - - y[1] := y[1] / mo2; - fac := 1.0; - for j := 1 to mo2 - 1 do - begin - k := 2 * j + 1; - if (fac <> 0) then - begin - fac := (1.0 - cnst * j * j) / mo2; - if (fac < 0.0) then fac := 0.0; - y[k] := fac * y[k]; - y[k + 1] := fac * y[k + 1]; - end - else - y[k + 1] := 0.0; - y[k] := 0.0; - end; - - fac := (1.0 - 0.25 * npts * npts) / mo2; - if (fac < 0.0) then fac := 0.0; - y[2] := y[2] * fac; - - RealFT(y,mo2,-1); - - for j := 1 to n do - y[j] := y[j] + rn1 * (y1 * (n - j) + yn * (j - 1)); - for j := 0 to n - 1 do - data[j] := y[j+1]; - - y := nil; -end; - -procedure TAutoCorrForm.PolyFit(const pts: DblDyneVec; var avg: DblDyneVec; - NoPts, Order: integer); -var - X: DblDyneMat; - XY: DblDyneVec; - XTX: DblDyneMat; - Beta: DblDyneVec; - t, Yhat: double; - i, j, k: integer; - RowLabels, ColLabels: StrDyneVec; -begin - SetLength(X, NoPts, order+1); - SetLength(XTX, order+2, order+2); - SetLength(XY, order+1); - SetLength(Beta, order+1); - SetLength(RowLabels, NoPts+1); - SetLength(ColLabels, NoPts+1); - - // initialize - for i := 0 to NoPts - 1 do - for j := 0 to order do - X[i,j] := 0.0; - - for i := 0 to order do - begin - XY[i] := 0.0; - Beta[i] := 0.0; - for j := 0 to order do - XTX[i,j] := 0.0; - end; - - for i := 0 to NoPts - 1 do - for j := 0 to order do - begin - t := i+1; - X[i,j] := Power(t,j); - end; - - // print the X matrix as a check - for i := 0 to NoPts - 1 do - RowLabels[i] := 'Case' + IntToStr(i+1); - - for i := 0 to order+1 do - ColLabels[i] := 'Order' + IntToStr(i); - - { - OutputFrm.RichEdit.Clear; - Title := 'X Matrix'; - DynMatPrint(X,NoPts,order+1,Title,RowLabels,ColLabels,NoPts); - OutputFrm.ShowModal; -} - // Get X transpose times X - for j := 0 to order do - begin - for k := 0 to order do - begin - XTX[j,k] := 0.0; - for i := 0 to NoPts - 1 do - XTX[j,k] := XTX[j,k] + (X[i,j] * X[i,k]); - end; - end; -{ - // print the XTX matrix - OutputFrm.RichEdit.Clear; - Title := 'XTX Matrix (Offset by 1)'; - DynMatPrint(XTX,order+2,order+2,Title,ColLabels,ColLabels,NoPts); - OutputFrm.ShowModal; -} - // Get X transpose Y - for j := 0 to order do - for i := 0 to NoPts - 1 do - XY[j] := XY[j] + (X[i,j] * pts[i]); -{ - // print the XY vector - OutputFrm.RichEdit.Clear; - Title := 'XY vector'; - DynVectorPrint(XY,order+1,Title,ColLabels,NoPts); - OutputFrm.ShowModal; -} - // get inverse of XTX - SVDinverse(XTX, order+1); -{ - // print the inverse matrix - OutputFrm.RichEdit.Clear; - Title := 'XTX Inverse Matrix'; - DynMatPrint(XTX,order+2,order+2,Title,ColLabels,ColLabels,NoPts); - OutputFrm.ShowModal; -} - // get betas - for j := 0 to order do - for k := 0 to order do - Beta[j] := Beta[j] + (XTX[j,k] * XY[k]); - - { - // print the betas - OutputFrm.RichEdit.Clear; - Title := 'Betas vector'; - DynVectorPrint(Beta,order+1,Title,ColLabels,NoPts); - OutputFrm.ShowModal; -} - // get predicted values - for i := 0 to NoPts - 1 do - begin - Yhat := 0.0; - t := i; - for j := 0 to order do - Yhat := Yhat + (Beta[j] * Power(t,j)); - avg[i] := Yhat; - end; - - //cleanup - ColLabels := nil; - RowLabels := nil; - Beta := nil; - XY := nil; - XTX := nil; - X := nil; -end; procedure TAutoCorrForm.UpdateBtnStates; begin @@ -1030,13 +1114,13 @@ begin FCorrelationMatrixReportFrame.UpdateBtnStates; FPartialsReportFrame.UpdateBtnStates; FYuleWalkerReportFrame.UpdateBtnStates; + FAutoRegChartFrame.UpdateBtnStates; InBtn.Enabled := (VarList.ItemIndex > -1) and (DepVarEdit.Text = ''); OutBtn.Enabled := (DepVarEdit.Text <> ''); end; - { Routines called from Compute } function TAutoCorrForm.CalcMean(const Pts: DblDyneVec; NoPts: Integer): Double; @@ -1049,164 +1133,6 @@ begin Result := Result / NoPts; end; -procedure TAutoCorrForm.ExponentialSmooth(var Pts: DblDyneVec; NoPts: Integer); -var - F: TExpSmoothFrm; - noProj: Integer; - avg: DblDyneVec = nil; - residual: DblDyneVec = nil; - i: Integer; - varRes: Double; - stdErr: Double; -begin - F := TExpSmoothFrm.Create(nil); - try - F.Alpha := FExpSmoothAlpha; - if F.ShowModal <> mrOK then - exit; - - FExpSmoothAlpha := F.Alpha; - - if ProjectChk.Checked then - noProj := StrToInt(ProjPtsEdit.Text) - else - noProj := 0; - - SetLength(avg, NoPts + noProj); - avg[0] := pts[0]; // set first value = to observed - for i := 1 to NoPts - 1 do // case pointer - avg[i] := FExpSmoothAlpha * pts[i] + (1.0 - FExpSmoothAlpha) * avg[i-1]; - - if ProjectChk.Checked then - begin - for i := 0 to noproj-1 do - begin - avg[NoPts + i] := FExpSmoothAlpha * pts[NoPts + i - 1]; - avg[NoPts + i] := avg[NoPts + i] + (1.0 - FExpSmoothAlpha) * avg[NoPts + i - 1]; - pts[NoPts + i] := avg[NoPts + i]; - end; - end; - - // plot the points - PointsFrm.pts := pts; - PointsFrm.avg := avg; - PointsFrm.NoCases := NoPts + noproj; - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Smoothed'; - PointsFrm.Title := 'Exponential Smoothed'; - PointsFrm.Caption := 'Exponential Smoothing'; - PointsFrm.ShowModal; - - if ResidChk.Checked then // calculate and plot residuals; - begin - SetLength(residual, NoPts); - varRes := 0.0; - for i := 0 to NoPts - 1 do - begin - residual[i] := pts[i] - avg[i]; - varRes := varRes + sqr(residual[i]); - end; - varRes := varRes / NoPts; - StdErr := sqrt(varRes); - - // plot the residuals - PointsFrm.pts := pts; - PointsFrm.avg := residual; - PointsFrm.NoCases := NoPts; - PointsFrm.MsgEdit.Text := 'Std. Err. Residuals: ' + FloatToStr(StdErr); - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Residuals'; - PointsFrm.Title := 'Residuals from Exponential Smoothing'; - PointsFrm.Caption := 'Exponential Residuals'; - PointsFrm.ShowModal; - end; - - // replace original points with smoothed values - for i := 0 to (NoPts + noproj - 1) do - pts[i] := avg[i]; - - finally - F.Free; - end; -end; - - -procedure TAutoCorrForm.FourierSmooth(var Pts: DblDyneVec; NoPts: Integer); -var - F: TFFTFrm; - avg: DblDyneVec; - residual: DblDyneVec; - NValues: Integer; - i: Integer; - noProj: Integer; - varRes: double; - stdErr: Double; -begin - if ProjectChk.Checked then - noProj := StrToInt(ProjPtsEdit.Text) - else - noProj := 0; - NValues := NoPts + noProj + 1; - - F := TFFtFrm.Create(nil); - try - F.NptsEdit.Text := IntToStr(NValues); - if F.ShowModal <> mrOK then - exit; - - SetLength(avg, NoPts + noProj); - if ProjectChk.Checked then - begin - for i := 0 to noproj - 1 do - begin - avg[i] := pts[NoPts-1-noproj+i]; - pts[i] := avg[i]; - end; - end; - - Fourier(avg, NValues, NValues); - - PointsFrm.pts := pts; - PointsFrm.avg := avg; - PointsFrm.NoCases := NoPts + noproj; - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Smoothed'; - PointsFrm.Title := 'Fourier Smoothed'; - PointsFrm.Caption := 'Fourier Smoothing'; - PointsFrm.ShowModal; - - // calculate and plot residuals; - if ResidChk.Checked then - begin - SetLength(residual, NoPts); - varRes := 0.0; - for i := 0 to NoPts - 1 do - begin - residual[i] := pts[i] - avg[i]; - varRes := varRes + sqr(residual[i]); - end; - varRes := varRes / NoPts; - stdErr := sqrt(varRes); - - // plot the residuals - PointsFrm.pts := pts; - PointsFrm.avg := residual; - PointsFrm.NoCases := NoPts; - PointsFrm.MsgEdit.Text := 'Std. Err. Residuals: ' + FloatToStr(StdErr); - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Residuals'; - PointsFrm.Title := 'Residuals from Fourier Smoothing'; - PointsFrm.Caption := 'Fourier Residuals'; - PointsFrm.ShowModal; - end; - - // replace original points with smoothed values - for i := 0 to (NoPts + noproj - 1) do pts[i] := avg[i]; - - finally - F.Free; - end; -end; procedure TAutoCorrForm.GetPts(var Pts: DblDyneVec; var ANumPts: Integer; DepVar: Integer); @@ -1263,259 +1189,56 @@ begin end; -procedure TAutoCorrForm.MovingAverage(var Pts: DblDyneVec; NoPts: Integer); +function TAutoCorrform.Prepare(out ADepVar: Integer; out APts: DblDyneVec): Boolean; var - F: TMoveAvgFrm; - i, j: Integer; - avg: DblDyneVec = nil; - residual: DblDyneVec = nil; - noProj: Integer; - varRes: Double; - stdErr: Double; -begin - F := TMoveAvgFrm.Create(nil); - try - F.Order := FMovAvgOrder; - F.RawWeights := FMovAvgRawWeights; - if F.ShowModal <> mrOK then - exit; - - FMovAvgRawWeights := F.RawWeights; - FMovAvgOrder := F.Order; - if FMovAvgOrder <= 0 then - exit; - - if ProjectChk.Checked then - noProj := StrToInt(ProjPtsEdit.Text) - else - noProj := 0; - - // Calculate moving average - SetLength(avg, NoPts + NoProj); - for i := F.Order to NoPts - F.Order - 1 do - begin - avg[i] := Pts[i] * F.Weights[0]; // middle value - for j := 1 to F.Order do - avg[i] := avg[i] + (Pts[i-j] + Pts[i+j]) * F.Weights[j]; - end; - - // fill in unestimable averages with original points - for i := 0 to F.Order - 1 do // left values - begin - avg[i] := pts[i] * F.Weights[0]; - for j := 1 to F.Order do - avg[i] := avg[i] + pts[i+j] * 2.0 * F.Weights[j]; - end; - for i := NoPts - F.Order to NoPts - 1 do //right values - begin - avg[i] := pts[i] * F.Weights[0]; - for j := 1 to F.Order do - avg[i] := avg[i] + (pts[i-j] * 2.0 * F.Weights[j]); - end; - - if ProjectChk.Checked then - begin - noProj := StrToInt(ProjPtsEdit.Text); - for i := 0 to NoProj-1 do - begin - avg[NoPts + i] := avg[NoPts - 1]; - pts[NoPts + i] := pts[NoPts - 1]; - end; - end; - - // plot the original points and the smoothed average - PointsFrm.pts := pts; - PointsFrm.avg := avg; - PointsFrm.NoCases := NoPts + noProj; - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Smoothed'; - PointsFrm.MsgEdit.Text := 'No. points: '+ IntToStr(NoPts); - PointsFrm.Title := 'Moving Average Smoothed'; - PointsFrm.Caption := 'Moving Average Smoothing'; - PointsFrm.ShowModal; - - // calculate and plot residuals; - if ResidChk.Checked then - begin - SetLength(residual, NoPts); - varRes := 0.0; - for i := 0 to NoPts - 1 do - begin - residual[i] := pts[i] - avg[i]; - varRes := varRes + sqr(residual[i]); - end; - varRes := varRes / NoPts; - stdErr := sqrt(varRes); - - // plot the residuals - PointsFrm.pts := pts; - PointsFrm.avg := residual; - PointsFrm.NoCases := NoPts; - PointsFrm.MsgEdit.Text := 'Std. Err. Residuals: ' + FloatToStr(StdErr); - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Residuals'; - PointsFrm.Title := 'Residuals from Moving Average'; - PointsFrm.Caption := 'Moving Average Residuals'; - PointsFrm.ShowModal; - end; - - // replace original points with smoothed values - for i := 0 to (NoPts + noproj - 1) do - pts[i] := avg[i]; - - finally - F.Free; - end; -end; - -procedure TAutoCorrForm.PlotDifferencesForLag(var Pts: DblDyneVec; NoPts: Integer); -var - F: TDifferenceFrm; - lag: Integer; - avg: DblDyneVec = nil; - residual: DblDyneVec = nil; - varRes: Double; - stdErr: Double; - i, j: Integer; -begin - F := TDifferenceFrm.Create(nil); - try - if F.ShowModal <> mrOK then - exit; - - lag := StrToInt(MaxLagEdit.Text); - - SetLength(avg, NoPts); - for i := 0 to NoPts-1 do - avg[i] := pts[i]; - - for j := 1 to StrToInt(F.OrderEdit.Text) do - for i := NoPts downto lag do - avg[i] := avg[i] - avg[i - lag]; // wp: what with avg[0]..avg[lag] ? - - // plot the original and differenced values - PointsFrm.pts := pts; - PointsFrm.avg := avg; - PointsFrm.NoCases := NoPts; - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Differenced'; - PointsFrm.MsgEdit.Text := 'No. points = ' + IntToStr(NoCases); - PointsFrm.Title := 'Differencing Smoothed'; - PointsFrm.Caption := 'Difference Smoothing'; - PointsFrm.ShowModal; - - // calculate and plot residuals; - if ResidChk.Checked then - begin - SetLength(residual, NoPts); - varRes := 0.0; - for i := 0 to NoPts - 1 do - begin - residual[i] := pts[i] - avg[i]; - varRes := varRes + sqr(residual[i]); - end; - varRes := varRes / NoPts; - stdErr := sqrt(varRes); - - // plot the residuals - PointsFrm.pts := pts; - PointsFrm.avg := residual; - PointsFrm.NoCases := NoPts; - PointsFrm.MsgEdit.Text := 'Std. Err. Residuals: ' + FloatToStr(StdErr); - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Residuals'; - PointsFrm.Title := 'Residuals from Difference Smoothing'; - PointsFrm.Caption := 'Difference Residuals'; - PointsFrm.ShowModal; - end; - - // replace original points by smoothed values - for i := 0 to NoPts - 1 do - pts[i] := avg[i]; - finally - F.Free; - end; -end; - -procedure TAutoCorrForm.PolynomialSmooth(var Pts: DblDyneVec; NoPts: Integer); -var - F: TPolynomialFrm; - noProj: Integer; + noPts: Integer; i: Integer; - avg: DblDyneVec = nil; - residual: DblDyneVec = nil; - varRes: Double; - stdErr: Double; begin - F := TPolynomialFrm.Create(nil); - try - F.Order := FPolyOrder; - if F.ShowModal <> mrOK then + Result := false; + + if ColBtn.Checked then + begin + // Get column of the selected variable + ADepVar := GetVariableIndex(OS3MainFrm.DataGrid, DepVarEdit.Text); + if (ADepVar = -1)then + begin + ErrorMsg('No variable selected to analyze.'); exit; - - FPolyOrder := F.Order; - - if ProjectChk.Checked then - noProj := StrToInt(ProjPtsEdit.Text) - else - noProj := 0; - - SetLength(avg, NoPts + noProj); - if ProjectChk.Checked then - begin - for i := 0 to noproj - 1 do - begin - avg[i] := pts[NoPts-1-noproj+i]; - pts[i] := avg[i]; - end; end; - - PolyFit(pts, avg, NoPts + noproj, FPolyOrder); - - // plot original and smoothed data - PointsFrm.pts := pts; - PointsFrm.avg := avg; - PointsFrm.NoCases := NoPts + noproj; - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Smoothed'; - PointsFrm.Title := 'Polynomial Smoothed'; - PointsFrm.Caption := 'Polynomial Smoothing'; - PointsFrm.ShowModal; - - // plot residuals if checked - if ResidChk.Checked then + end + else + begin + // Get row of the selected case + ADepVar := -1; + for i := 1 to NoCases do begin - SetLength(residual, NoPts); - varRes := 0.0; - for i := 0 to NoPts - 1 do - begin - residual[i] := pts[i] - avg[i]; - varRes := varRes + sqr(residual[i]); - end; - varRes := varRes / NoPts; - stdErr := sqrt(varRes); - - // plot the residuals - PointsFrm.pts := pts; - PointsFrm.avg := residual; - PointsFrm.NoCases := NoPts; - PointsFrm.MsgEdit.Text := 'Std. Err. Residuals: ' + FloatToStr(StdErr); - PointsFrm.LabelOne := 'Original'; - PointsFrm.LabelTwo := 'Residuals'; - PointsFrm.Title := 'Residuals from Polynomial Smoothing'; - PointsFrm.Caption := 'Polynomial Residuals'; - PointsFrm.ShowModal; +// if not GoodRecord(OS3MainFrm.DataGrid, i, ColNoSelected) then continue; + if (OS3MainFrm.DataGrid.Cells[0,i] = DepVarEdit.Text) then ADepVar := i; + end; + if (ADepVar = -1)then + begin + ErrorMsg('No variable selected to analyze.'); + exit; end; - - // replace original points with smoothed values - for i := 0 to (NoPts + noproj - 1) do - pts[i] := avg[i]; - - finally - F.Free; end; + (* + if ProjectChk.Checked then + noProj := StrToInt(ProjPtsEdit.Text) + else + noProj := 0; *) + + // Get points to analyze + GetPts(APts, noPts, ADepVar); + if noPts = 0 then + begin + ErrorMsg('No data points found.'); + exit; + end; + + Result := true; end; + procedure TAutoCorrForm.RemoveMeans(var Pts: DblDyneVec; NoPts: Integer; AMean: Double); var i: Integer; @@ -1525,6 +1248,40 @@ begin end; +procedure TAutoCorrForm.SmoothingParamsBtnClick(Sender: TObject); +var + msg: String; + C: TWinControl; + depVar: Integer; + pts: DblDyneVec; + res: TModalResult; +begin + if not Validate(msg, C) then + begin + C.SetFocus; + ErrorMsg(msg); + exit; + end; + + Prepare(DepVar, pts); + + if SmoothingForm = nil then + SmoothingForm := TSmoothingForm.Create(Application); + + SmoothingForm.Points := pts; + + if ProjectChk.Checked then + SmoothingForm.NumProj := StrToInt(ProjPtsEdit.Text) + else + SmoothingForm.NumProj := 0; + SmoothingForm.Params := FParams; + + res := SmoothingForm.ShowModal; + if res = mrOK then + FParams := SmoothingForm.Params; +end; + + function TAutoCorrForm.Validate(out AMsg: String; out AControl: TWinControl): Boolean; var x: Double; diff --git a/applications/lazstats/source/forms/analysis/correlation/expsmoothunit.pas b/applications/lazstats/source/forms/analysis/correlation/expsmoothunit.pas index cea7f4ff5..8644ba68d 100644 --- a/applications/lazstats/source/forms/analysis/correlation/expsmoothunit.pas +++ b/applications/lazstats/source/forms/analysis/correlation/expsmoothunit.pas @@ -5,7 +5,7 @@ unit ExpSmoothUnit; interface uses - Classes, SysUtils, FileUtil, LResources, Forms, Controls, Graphics, Dialogs, + Classes, SysUtils, Forms, Controls, Graphics, Dialogs, StdCtrls, ExtCtrls; type @@ -36,8 +36,11 @@ type var ExpSmoothFrm: TExpSmoothFrm; + implementation +{$R *.lfm} + uses Math; @@ -80,8 +83,5 @@ begin end; -initialization - {$I expsmoothunit.lrs} - end. diff --git a/applications/lazstats/source/forms/analysis/correlation/smoothingunit.lfm b/applications/lazstats/source/forms/analysis/correlation/smoothingunit.lfm new file mode 100644 index 000000000..8f578b506 --- /dev/null +++ b/applications/lazstats/source/forms/analysis/correlation/smoothingunit.lfm @@ -0,0 +1,491 @@ +inherited SmoothingForm: TSmoothingForm + Left = 603 + Height = 463 + Top = 292 + Width = 823 + Caption = 'SmoothingForm' + ClientHeight = 463 + ClientWidth = 823 + inherited ParamsPanel: TPanel + Height = 447 + Width = 296 + ClientHeight = 447 + ClientWidth = 296 + object OKBtn: TButton[0] + AnchorSideRight.Control = ParamsPanel + AnchorSideRight.Side = asrBottom + AnchorSideBottom.Control = ParamsPanel + AnchorSideBottom.Side = asrBottom + Left = 254 + Height = 25 + Top = 422 + Width = 42 + Anchors = [akRight, akBottom] + AutoSize = True + BorderSpacing.Left = 8 + BorderSpacing.Top = 8 + Caption = 'OK' + ModalResult = 1 + TabOrder = 9 + end + inherited CloseBtn: TButton[1] + AnchorSideRight.Control = HelpBtn + AnchorSideRight.Side = asrTop + Left = -88 + Top = 422 + ModalResult = 11 + Visible = False + end + inherited ComputeBtn: TButton[2] + AnchorSideRight.Control = OKBtn + Left = 88 + Top = 422 + Anchors = [akTop] + Visible = False + end + inherited ResetBtn: TButton[3] + Left = 26 + Top = 422 + Visible = False + end + inherited HelpBtn: TButton[4] + Left = -33 + Top = 422 + Visible = False + end + inherited ButtonBevel: TBevel[5] + AnchorSideBottom.Control = OKBtn + Top = 406 + Width = 296 + end + object DataSmoothingChkLB: TCheckListBox[6] + AnchorSideLeft.Control = ParamsPanel + AnchorSideTop.Control = Label2 + AnchorSideTop.Side = asrBottom + AnchorSideRight.Control = ParamsPanel + AnchorSideRight.Side = asrBottom + AnchorSideBottom.Control = Splitter1 + Left = 0 + Height = 129 + Top = 64 + Width = 296 + Anchors = [akTop, akLeft, akRight, akBottom] + BorderSpacing.Top = 2 + BorderSpacing.Bottom = 8 + Items.Strings = ( + 'Difference Smoothing' + 'Moving Average Smoothing' + 'Exponential Smoothing' + 'Fourier Filter Smoothing' + 'Polynomial Regression Smooth' + 'Multiple Regression Smoothing' + ) + ItemHeight = 17 + OnSelectionChange = DataSmoothingChkLBSelectionChange + TabOrder = 4 + Data = { + 06000000000000000000 + } + end + object Label2: TLabel[7] + AnchorSideLeft.Control = ParamsPanel + AnchorSideTop.Control = Bevel1 + AnchorSideTop.Side = asrBottom + Left = 0 + Height = 15 + Top = 47 + Width = 136 + BorderSpacing.Top = 12 + Caption = 'Data Smoothing Methods' + ParentColor = False + end + object MeanChk: TCheckBox[8] + AnchorSideLeft.Control = ParamsPanel + AnchorSideTop.Control = ParamsPanel + Left = 0 + Height = 19 + Top = 0 + Width = 105 + BorderSpacing.Bottom = 8 + Caption = 'Center on Mean' + OnChange = MeanChkChange + TabOrder = 5 + end + object ResidChk: TCheckBox[9] + AnchorSideLeft.Control = MeanChk + AnchorSideLeft.Side = asrBottom + AnchorSideTop.Control = ParamsPanel + Left = 121 + Height = 19 + Top = 0 + Width = 88 + BorderSpacing.Left = 16 + Caption = 'Residual plot' + OnChange = ResidChkChange + TabOrder = 6 + end + object Bevel1: TBevel[10] + AnchorSideLeft.Control = ParamsPanel + AnchorSideTop.Control = MeanChk + AnchorSideTop.Side = asrBottom + AnchorSideRight.Control = ParamsPanel + AnchorSideRight.Side = asrBottom + Left = 0 + Height = 8 + Top = 27 + Width = 296 + Anchors = [akTop, akLeft, akRight] + Shape = bsBottomLine + end + object NotebookGroup: TGroupBox[11] + AnchorSideLeft.Control = ParamsPanel + AnchorSideTop.Control = Splitter1 + AnchorSideTop.Side = asrBottom + AnchorSideRight.Control = ParamsPanel + AnchorSideRight.Side = asrBottom + AnchorSideBottom.Control = ButtonBevel + Left = 0 + Height = 188 + Top = 218 + Width = 296 + Anchors = [akTop, akLeft, akRight, akBottom] + BorderSpacing.Top = 12 + Caption = 'NotebookGroup' + ClientHeight = 168 + ClientWidth = 292 + TabOrder = 7 + object Notebook: TNotebook + AnchorSideLeft.Control = NotebookGroup + AnchorSideTop.Control = NotebookGroup + AnchorSideRight.Control = NotebookGroup + AnchorSideRight.Side = asrBottom + AnchorSideBottom.Control = NotebookGroup + AnchorSideBottom.Side = asrBottom + Left = 12 + Height = 156 + Top = 4 + Width = 268 + PageIndex = 0 + Anchors = [akTop, akLeft, akRight, akBottom] + BorderSpacing.Left = 12 + BorderSpacing.Top = 4 + BorderSpacing.Right = 12 + BorderSpacing.Bottom = 8 + TabOrder = 0 + object DifferencePage: TPage + object Label8: TLabel + AnchorSideTop.Side = asrCenter + Left = 13 + Height = 15 + Top = 4 + Width = 108 + Anchors = [akTop, akRight] + BorderSpacing.Right = 8 + Caption = 'Difference for lag of:' + ParentColor = False + end + object Label9: TLabel + AnchorSideLeft.Control = DifferencePage + AnchorSideTop.Side = asrCenter + Left = 0 + Height = 15 + Top = 31 + Width = 121 + BorderSpacing.Right = 8 + Caption = 'No. of times to repeat: ' + ParentColor = False + end + object LagEdit: TSpinEdit + AnchorSideLeft.Control = Label8 + AnchorSideLeft.Side = asrBottom + AnchorSideTop.Control = DifferencePage + Left = 129 + Height = 23 + Top = 0 + Width = 66 + Alignment = taRightJustify + MaxValue = 65535 + OnChange = LagEditChange + TabOrder = 0 + end + object DiffOrderEdit: TSpinEdit + AnchorSideLeft.Control = Label9 + AnchorSideLeft.Side = asrBottom + Left = 129 + Height = 23 + Top = 28 + Width = 66 + Alignment = taRightJustify + MaxValue = 65535 + MinValue = 1 + OnChange = DiffOrderEditChange + TabOrder = 1 + Value = 1 + end + end + object MovingAvgPage: TPage + object Label3: TLabel + AnchorSideLeft.Control = MovingAvgPage + AnchorSideTop.Control = WeightOrder + AnchorSideTop.Side = asrCenter + Left = 0 + Height = 15 + Top = 4 + Width = 36 + BorderSpacing.Right = 8 + Caption = 'Order: ' + ParentColor = False + end + object WeightsGrid: TStringGrid + AnchorSideLeft.Control = MovingAvgPage + AnchorSideTop.Side = asrBottom + AnchorSideRight.Control = MovingAvgPage + AnchorSideRight.Side = asrBottom + AnchorSideBottom.Control = MovingAvgPage + AnchorSideBottom.Side = asrBottom + Left = 0 + Height = 124 + Top = 32 + Width = 260 + Anchors = [akTop, akLeft, akRight, akBottom] + AutoAdvance = aaDown + AutoFillColumns = True + BorderSpacing.Top = 8 + BorderSpacing.Right = 8 + ColCount = 3 + Options = [goFixedVertLine, goFixedHorzLine, goVertLine, goHorzLine, goRangeSelect, goDrawFocusSelected, goEditing, goThumbTracking, goSmoothScroll, goFixedColSizing] + RowCount = 2 + TabOrder = 0 + OnEditingDone = WeightsGridEditingDone + OnPrepareCanvas = WeightsGridPrepareCanvas + OnSelectEditor = WeightsGridSelectEditor + ColWidths = ( + 78 + 91 + 91 + ) + Cells = ( + 3 + 0 + 0 + 'Weight #' + 1 + 0 + 'Weight Value' + 2 + 0 + 'Normalized' + ) + end + object WeightOrder: TSpinEdit + AnchorSideLeft.Control = Label3 + AnchorSideLeft.Side = asrBottom + AnchorSideTop.Control = MovingAvgPage + Left = 44 + Height = 23 + Top = 0 + Width = 68 + Alignment = taRightJustify + OnChange = WeightOrderChange + OnEditingDone = WeightOrderEditingDone + TabOrder = 1 + end + end + object ExpPage: TPage + object Label5: TLabel + AnchorSideLeft.Side = asrBottom + AnchorSideTop.Control = AlphaEdit + AnchorSideTop.Side = asrCenter + AnchorSideRight.Control = AlphaEdit + Left = 56 + Height = 15 + Top = 12 + Width = 31 + Anchors = [akTop, akRight] + BorderSpacing.Left = 16 + BorderSpacing.Right = 8 + Caption = 'Alpha' + ParentColor = False + end + object AlphaScroll: TScrollBar + AnchorSideLeft.Control = ExpPage + AnchorSideTop.Side = asrBottom + AnchorSideRight.Control = ExpPage + AnchorSideRight.Side = asrBottom + Left = 16 + Height = 23 + Top = 39 + Width = 236 + Anchors = [akTop, akLeft, akRight] + BorderSpacing.Left = 16 + BorderSpacing.Top = 12 + BorderSpacing.Right = 16 + PageSize = 0 + Position = 99 + TabOrder = 0 + OnChange = AlphaScrollChange + end + object Label6: TLabel + AnchorSideLeft.Control = AlphaScroll + AnchorSideTop.Control = AlphaScroll + AnchorSideTop.Side = asrBottom + Left = 24 + Height = 15 + Top = 66 + Width = 15 + BorderSpacing.Left = 8 + BorderSpacing.Top = 4 + Caption = '0.0' + ParentColor = False + end + object Label7: TLabel + AnchorSideTop.Control = AlphaScroll + AnchorSideTop.Side = asrBottom + AnchorSideRight.Control = AlphaScroll + AnchorSideRight.Side = asrBottom + Left = 229 + Height = 15 + Top = 66 + Width = 15 + Anchors = [akTop, akRight] + BorderSpacing.Top = 4 + BorderSpacing.Right = 8 + Caption = '1.0' + ParentColor = False + end + object AlphaEdit: TFloatSpinEdit + AnchorSideLeft.Control = ExpPage + AnchorSideLeft.Side = asrCenter + Left = 95 + Height = 23 + Top = 8 + Width = 78 + Alignment = taRightJustify + Increment = 0.01 + MaxValue = 1 + OnChange = AlphaEditChange + TabOrder = 1 + Value = 0.01 + end + end + object FourierPage: TPage + object Label4: TLabel + AnchorSideLeft.Control = FourierPage + AnchorSideTop.Control = StrengthEdit + AnchorSideTop.Side = asrCenter + Left = 0 + Height = 15 + Top = 16 + Width = 76 + BorderSpacing.Right = 8 + Caption = 'Filter strength:' + ParentColor = False + end + object StrengthEdit: TSpinEdit + AnchorSideLeft.Control = Label4 + AnchorSideLeft.Side = asrBottom + AnchorSideTop.Control = FourierPage + AnchorSideRight.Side = asrBottom + Left = 84 + Height = 23 + Top = 12 + Width = 68 + Alignment = taRightJustify + Anchors = [akTop, akLeft, akRight] + BorderSpacing.Top = 12 + OnChange = StrengthEditChange + OnEditingDone = StrengthEditEditingDone + TabOrder = 0 + end + end + object PolyPage: TPage + AnchorSideTop.Side = asrBottom + object Label1: TLabel + AnchorSideLeft.Control = PolyPage + AnchorSideTop.Side = asrCenter + Left = 0 + Height = 15 + Top = 4 + Width = 97 + BorderSpacing.Right = 8 + Caption = 'Polynomial order :' + ParentColor = False + end + object PolyOrderEdit: TSpinEdit + AnchorSideLeft.Control = Label1 + AnchorSideLeft.Side = asrBottom + Left = 105 + Height = 23 + Top = 0 + Width = 50 + Alignment = taRightJustify + MinValue = 1 + OnChange = PolyOrderEditChange + TabOrder = 0 + Value = 1 + end + end + object MRegPage: TPage + end + end + end + object Splitter1: TSplitter[12] + AnchorSideLeft.Control = ParamsPanel + AnchorSideRight.Control = ParamsPanel + AnchorSideRight.Side = asrBottom + Cursor = crVSplit + Left = 0 + Height = 5 + Top = 201 + Width = 296 + Align = alNone + ResizeAnchor = akBottom + end + object CancelBtn: TButton[13] + AnchorSideTop.Control = OKBtn + AnchorSideRight.Control = OKBtn + Left = 184 + Height = 25 + Top = 422 + Width = 62 + Anchors = [akTop, akRight] + AutoSize = True + BorderSpacing.Left = 8 + BorderSpacing.Right = 8 + Caption = 'Cancel' + ModalResult = 2 + TabOrder = 10 + end + end + inherited ParamsSplitter: TSplitter + Left = 308 + Height = 463 + end + object ButtonBevel2: TBevel[2] + AnchorSideLeft.Control = ParamsPanel + AnchorSideRight.Control = ParamsPanel + AnchorSideRight.Side = asrBottom + AnchorSideBottom.Control = CloseBtn + Left = 8 + Height = 8 + Top = 406 + Width = 296 + Anchors = [akLeft, akRight, akBottom] + Shape = bsBottomLine + end + object OriginalDataChartSource: TUserDefinedChartSource[3] + OnGetChartDataItem = OriginalDataChartSourceGetChartDataItem + Left = 429 + Top = 68 + end + object SmoothedDataChartSource: TUserDefinedChartSource[4] + OnGetChartDataItem = SmoothedDataChartSourceGetChartDataItem + Left = 430 + Top = 122 + end + object ResidualsChartSource: TUserDefinedChartSource[5] + OnGetChartDataItem = ResidualsChartSourceGetChartDataItem + Left = 433 + Top = 182 + end +end diff --git a/applications/lazstats/source/forms/analysis/correlation/smoothingunit.pas b/applications/lazstats/source/forms/analysis/correlation/smoothingunit.pas new file mode 100644 index 000000000..60ce023ab --- /dev/null +++ b/applications/lazstats/source/forms/analysis/correlation/smoothingunit.pas @@ -0,0 +1,1048 @@ +unit SmoothingUnit; + +{$mode objfpc}{$H+} +{$modeswitch advancedrecords} + +interface + +uses + Classes, SysUtils, Forms, Controls, Graphics, Dialogs, StdCtrls, ExtCtrls, + CheckLst, Spin, Grids, TASeries, TASources, Globals, BasicStatsChartFormUnit, + TACustomSource; + +type + TSmoothingDiffParams = record + Active: Boolean; + MaxLag: Integer; + Order: Integer; + end; + + TSmoothingMovingAvgParams = record + Active: Boolean; + RawWeights: DblDyneVec; + Weights: DblDyneVec; + procedure NormalizeWeights; + function Order: Integer; + end; + + TSmoothingExpParams = record + Active: Boolean; + Alpha: Double; + end; + + TSmoothingFourierParams = record + Active: Boolean; + Strength: Integer; + end; + + TSmoothingPolynomialParams = record + Active: Boolean; + Order: Integer; + end; + + TSmoothingParams = record + CenterOnMean: Boolean; + DifferenceParams: TSmoothingDiffParams; + MovingAvgParams: TSmoothingMovingAvgParams; + ExponentialParams: TSmoothingExpParams; + FourierParams: TSmoothingFourierParams; + PolynomialParams: TSmoothingPolynomialParams; + procedure CopyFrom(AParams: TSmoothingParams); + end; + + { TSmoothingForm } + + TSmoothingForm = class(TBasicStatsChartForm) + AlphaScroll: TScrollBar; + Bevel1: TBevel; + CancelBtn: TButton; + OKBtn: TButton; + DataSmoothingChkLB: TCheckListBox; + DifferencePage: TPage; + AlphaEdit: TFloatSpinEdit; + Label3: TLabel; + Label4: TLabel; + Label5: TLabel; + Label6: TLabel; + Label7: TLabel; + Label8: TLabel; + Label9: TLabel; + NotebookGroup: TGroupBox; + Label2: TLabel; + Label1: TLabel; + MovingAvgPage: TPage; + ExpPage: TPage; + FourierPage: TPage; + MRegPage: TPage; + PolyPage: TPage; + ButtonBevel2: TBevel; + MeanChk: TCheckBox; + Notebook: TNotebook; + ResidChk: TCheckBox; + OriginalDataChartSource: TUserDefinedChartSource; + SmoothedDataChartSource: TUserDefinedChartSource; + ResidualsChartSource: TUserDefinedChartSource; + PolyOrderEdit: TSpinEdit; + LagEdit: TSpinEdit; + DiffOrderEdit: TSpinEdit; + Splitter1: TSplitter; + StrengthEdit: TSpinEdit; + WeightOrder: TSpinEdit; + WeightsGrid: TStringGrid; + procedure AlphaEditChange(Sender: TObject); + procedure AlphaScrollChange(Sender: TObject); + procedure DataSmoothingChkLBSelectionChange(Sender: TObject; {%H-}User: boolean); + procedure DiffOrderEditChange(Sender: TObject); + procedure LagEditChange(Sender: TObject); + procedure MeanChkChange(Sender: TObject); + procedure StrengthEditChange(Sender: TObject); + procedure StrengthEditEditingDone(Sender: TObject); + procedure PolyOrderEditChange(Sender: TObject); + procedure ResidChkChange(Sender: TObject); + procedure ResidualsChartSourceGetChartDataItem({%H-}ASource: TUserDefinedChartSource; + AIndex: Integer; var AItem: TChartDataItem); + procedure SmoothedDataChartSourceGetChartDataItem({%H-}ASource: TUserDefinedChartSource; + AIndex: Integer; var AItem: TChartDataItem); + procedure WeightOrderChange(Sender: TObject); + procedure WeightOrderEditingDone(Sender: TObject); + procedure OriginalDataChartSourceGetChartDataItem({%H-}ASource: TUserDefinedChartSource; + AIndex: Integer; var AItem: TChartDataItem); + procedure WeightsGridEditingDone(Sender: TObject); + procedure WeightsGridPrepareCanvas({%H-}Sender: TObject; {%H-}aCol, {%H-}aRow: Integer; + {%H-}aState: TGridDrawState); + procedure WeightsGridSelectEditor(Sender: TObject; aCol, {%H-}aRow: Integer; + var Editor: TWinControl); + + private + FParams: TSmoothingParams; + FNumProj: Integer; + FOrigPoints: DblDyneVec; + FSmoothedPoints: DblDyneVec; + FMean: Double; + FDataSeries: TLineSeries; + FSmoothedSeries: TLineSeries; + FResidualSeries: TLineSeries; + procedure SetParams(const AParams: TSmoothingParams); + procedure SetPoints(const APoints: DblDyneVec); + procedure SetRawWeights(const AWeights: DblDyneVec); + procedure WeightsToGrid; + + protected + procedure AdjustConstraints; override; + procedure Compute; override; + function Validate(out AMsg: String; out AControl: TWinControl): boolean; override; + + public + constructor Create(AOwner: TComponent); override; + + property Params: TSmoothingParams read FParams write SetParams; + property Points: DblDyneVec read FOrigPoints write SetPoints; + property NumProj: integer read FNumProj write FNumProj; + + end; + +var + SmoothingForm: TSmoothingForm; + + +const + // Default to "no smoothing" + DEFAULT_SMOOTHING_PARAMS: TSmoothingParams = ( + CenterOnMean: false; + DifferenceParams: (Active: false; MaxLag:1; Order:1); + MovingAvgParams: (Active: false; RawWeights: nil; Weights: nil); + ExponentialParams: (Active: false; Alpha: 1.00); + FourierParams: (Active: false; Strength: 0); + PolynomialParams: (Active: false; Order: 2) + ); + +procedure DifferenceSmooth(var AData: DblDyneVec; AOrder, AMaxLag: Integer); +procedure ExponentialSmooth(var AData: DblDyneVec; Alpha: Double; ANumProj: Integer); +procedure MovingAvgSmooth(var AData: DblDyneVec; AWeights: DblDyneVec; + AOrder, ANumProj: Integer); +procedure FourierSmooth(var AData: DblDyneVec; AStrength, ANumProj: integer); +procedure PolynomialSmooth(var AData: DblDyneVec; APolyOrder: Integer); + + +implementation + +{$R *.lfm} + +uses + Math, + Utils, MathUnit, MatrixUnit, ChartFrameUnit; + + +// Normalize all values so that their sum is 1. +// Note that the weights are considered to be symmetric around index 0. +procedure TSmoothingMovingAvgParams.NormalizeWeights; +var + sum: Double; + i: Integer; +begin + if (RawWeights = nil) or (Length(RawWeights) = 0) then + begin + SetLength(RawWeights, 1); + RawWeights[0] := 1.0; + end; + + sum := RawWeights[0]; + for i := 1 to High(RawWeights) do + sum := sum + RawWeights[i] * 2; + + SetLength(Weights, Length(RawWeights)); + if sum = 0 then + begin + Weights[0] := 1.0; + for i := 1 to High(Weights) do + Weights[i] := 0; + end else + for i := 0 to High(RawWeights) do + Weights[i] := RawWeights[i] / sum; +end; + +function TSmoothingMovingAvgParams.Order: Integer; +begin + Result := Length(RawWeights) - 1; +end; + + +procedure TSmoothingParams.CopyFrom(AParams: TSmoothingParams); +begin + CenterOnMean := AParams.CenterOnMean; + DifferenceParams := AParams.DifferenceParams; + MovingAvgParams := AParams.MovingAvgParams; + ExponentialParams := AParams.ExponentialParams; + FourierParams := AParams.FourierParams; + PolynomialParams := AParams.PolynomialParams; +end; + + +{ Smoothing by difference calculation } + +procedure DifferenceSmooth(var AData: DblDyneVec; AOrder, AMaxLag: Integer); +var + i, j, lag, numPts: Integer; + avg: DblDyneVec = nil; +begin + numPts := Length(AData); + if numPts = 0 then + exit; + + lag := AMaxLag; + + avg := VecCopy(AData); + + for j := 1 to AOrder do + for i := numPts-1 downto lag do + avg[i] := avg[i] - avg[i - lag]; + + // replace original points by smoothed values + AData := VecCopy(avg); +end; + + + +{ Exponential smoothing } + +procedure ExponentialSmooth(var AData: DblDyneVec; Alpha: Double; ANumProj: Integer); +var + avg: DblDyneVec = nil; + i, j, numPts: Integer; +begin + numPts := Length(AData); + if numPts = 0 then + exit; + + SetLength(avg, numPts + ANumProj); + avg[0] := AData[0]; // set first value = to observed + for i := 1 to numPts - 1 do // case pointer + avg[i] := Alpha * AData[i] + (1.0 - Alpha) * avg[i-1]; + + if ANumProj > 0 then + begin + for i := 0 to ANumProj-1 do + begin + j := numPts + i; + avg[j] := Alpha * AData[j-1] + (1 - Alpha) * avg[j-1]; + end; + SetLength(AData, numPts + ANumProj); + end; + + // replace original points by smoothed values + for i := 0 to (numPts + ANumProj - 1) do + AData[i] := avg[i]; +end; + + +{ Smoothing by moving average calculation } + +procedure MovingAvgSmooth(var AData: DblDyneVec; AWeights: DblDyneVec; + AOrder, ANumProj: Integer); +var + i, j: Integer; + avg: DblDyneVec = nil; + w: Double; + numPts: Integer; +begin + numPts := Length(AData); + if numPts = 0 then + exit; + + // Calculate moving average + + SetLength(avg, numPts + ANumProj); + for i := AOrder to numPts - AOrder - 1 do + begin + w := AWeights[0]; + avg[i] := AData[i] * w; // center value + for j := 1 to AOrder do + begin + w := AWeights[j]; + avg[i] := avg[i] + (AData[i-j] + AData[i+j]) * w; + end; + end; + + // Fill in unestimable averages with original points + + // Left values + for i := 0 to AOrder - 1 do + begin + w := AWeights[0]; + avg[i] := AData[i] * w; + for j := 1 to AOrder do + begin + w := AWeights[j]; + avg[i] := avg[i] + AData[i+j] * 2.0 * w; + end; + end; + + // Right values + for i := numPts - AOrder to numPts - 1 do + begin + w := AWeights[0]; + avg[i] := AData[i] * w; + for j := 1 to AOrder do + begin + w := AWeights[j]; + avg[i] := avg[i] + (AData[i-j] * 2.0 * w); + end; + end; + + if ANumProj > 0 then + begin + for i := 0 to ANumProj-1 do + avg[numPts + i] := avg[numPts - 1]; + SetLength(AData, numPts + ANumProj); + end; + + // replace original points with smoothed values + for i := 0 to numPts + ANumProj - 1 do + AData[i] := avg[i]; +end; + + +{ Smoothing by fitting a polynomial of given order } + +procedure PolynomialSmooth(var AData: DblDyneVec; APolyOrder: Integer); +var + X: DblDyneMat = nil; + XT: DblDyneMat = nil; + XTX: DblDyneMat = nil; + XTXinv: DblDyneMat = nil; + Beta: DblDyneVec = nil; + XY: DblDyneVec = nil; + numPts: Integer; + i, j: Integer; + t: Double; +begin + numPts := Length(AData); + if numPts = 0 then + exit; + + // Get polynomial input matrix + SetLength(X, numPts, APolyOrder + 1); + for i := 0 to numPts-1 do + begin + t := i + 1; + for j := 0 to APolyOrder do + begin + X[i, j] := IntPower(t, j); + end; + end; + + // Transposed matrix of X + XT := MatTransposed(X); + + // Calculate product XTransposed times X + XTX := XT * X; + + // Calulate inverse marix of XTX + XTXinv := MatInverse(XTX); + + // Calculate product xTransposed times Y; + XY := XT * AData; + + // Get fit coefficients (beta) + Beta := XTXInv * XY; + + // Get predicted values and return as APoints + AData := X * Beta; +end; + +{ Smoothing by Fourier transformation } + +{ Basic FFT routine, ref: Numerical Recipes + AData is an array 1..2*nn } +procedure Four1(var AData: TDblVector; nn,isign: integer); +var + ii, jj, n, mmax, m, j, istep, i: integer; + wtemp, wr, wpr, wpi, wi, theta: double; + tempr, tempi: double; +begin + n := 2*nn; + j := 1; + for ii := 1 to nn do + begin + i := 2*ii - 1; + if (j > i) then + begin + Exchange(AData[j], AData[i]); + Exchange(AData[j+1], AData[i+1]); + end; + m := n div 2; + while ((m >= 2) and (j > m)) do + begin + j := j - m; + m := m DIV 2 + end; + j := j+m + end; + mmax := 2; + while (n > mmax) do + begin + istep := 2*mmax; + theta := TWO_PI / (isign*mmax); + wpr := -2.0*sqr(sin(0.5*theta)); + wpi := sin(theta); + wr := 1.0; + wi := 0.0; + for ii := 1 to mmax div 2 do + begin + m := 2*ii - 1; + for jj := 0 to ((n-m) div istep) do + begin + i := m + jj*istep; + j := i + mmax; + tempr := wr*AData[j] - wi*AData[j+1]; + tempi := wr*AData[j+1] + wi*AData[j]; + AData[j] := AData[i] - tempr; + AData[j+1] := AData[i+1] - tempi; + AData[i] := AData[i] + tempr; + AData[i+1] := AData[i+1] + tempi + end; + wtemp := wr; + wr := wr*wpr - wi*wpi + wr; + wi := wi*wpr + wtemp*wpi + wi; + end; + mmax := istep; + end; +end; + +{ FFT of an array with real values, ref: Numerical Recipes. + AData is an array 1..2*n. n must be power of 2 + isign = +1 ---> forward transform (time to frequency) + isign = -1 ---> backward transform (frequency to time) + Input array is overwriten by procedure. Note that the calling routine must + divide it by n in the forward transformation case. } +procedure RealFT(var AData: TDblVector; n, isign: integer); +var + wr, wi, wpr, wpi, wtemp, theta: double; + i, i1, i2, i3, i4: integer; + c1, c2, h1r, h1i, h2r, h2i, wrs, wis: double; +begin + theta := TWO_PI / (2*n); + c1 := 0.5; + if (isign = 1) then + begin + c2 := -0.5; + Four1(AData, n, 1); + end else + begin + c2 := 0.5; + theta := -theta; + end; + wpr := -2.0*sqr(sin(0.5*theta)); + wpi := sin(theta); + wr := 1.0 + wpr; + wi := wpi; + for i := 2 to (n div 2) do // case i=1 done separately below + begin + i1 := i + i - 1; + i2 := i1 + 1; + i3 := n + n + 3 - i2; + i4 := i3 + 1; + wrs := wr; //sngl(wr); + wis := wi; //sngl(wi); + h1r := c1 * (AData[i1] + AData[i3]); + h1i := c1 * (AData[i2] - AData[i4]); + h2r := -c2 * (AData[i2] + AData[i4]); + h2i := c2 * (AData[i1] - AData[i3]); + AData[i1] := h1r + wrs*h2r - wis*h2i; + AData[i2] := h1i + wrs*h2i + wis*h2r; + AData[i3] := h1r - wrs*h2r + wis*h2i; + AData[i4] := -h1i + wrs*h2i + wis*h2r; + wtemp := wr; + wr := wr*wpr - wi*wpi + wr; + wi := wi*wpr + wtemp*wpi + wi + end; + if (isign = 1) then + begin + h1r := AData[1]; + AData[1] := h1r + AData[2]; + AData[2] := h1r - AData[2] + end else + begin + h1r := AData[1]; + AData[1] := c1 * (h1r + AData[2]); + AData[2] := c1 * (h1r - AData[2]); + Four1(AData, n, -1) + end; +end; + + +// NOTE: AData is a 0-based array (opposed to the other Fourier routines } +procedure DoFourierSmooth(var AData: DblDyneVec; n, AStrength: integer); +var + mmin, m, mo2, i, k, j: integer; + fac, cnst: double; + y: DblDyneVec = nil; +begin + m := 2; + mmin := n; // + 2*npts; + while (m < mmin) do m := m * 2; + mo2 := m div 2; // "m over 2" + + // y is a 1-based copy of the AData array, as required by the aux routines. + // y must be over-dimensioned by factor 2 + SetLength(y, m+1); + for i := 0 to n - 1 do + y[i+1] := AData[i]; + + { wp: why this? + y1 := y[1]; + yn := y[n]; + rn1 := 1.0 / (n - 1); + for j := 1 to n do + y[j] := y[j] + (-rn1 * (y1 * (n - j) + y1 * (j - 1))); + for j := n+1 to m do + y[j] := 0.0; + } + + // Transformation to frequency domain + RealFT(y, mo2, 1); + + // Filtering: Multiplies frequency spectrum by a parabola + cnst := sqr(AStrength / m); + y[1] := y[1] / mo2; + fac := 1.0; + for j := 1 to mo2 - 1 do + begin + k := 2 * j + 1; + if (fac <> 0) then + begin + fac := (1.0 - cnst * j * j) / mo2; + if (fac < 0.0) then fac := 0.0; + y[k] := fac * y[k]; + y[k + 1] := fac * y[k + 1]; + end + else begin + y[k + 1] := 0.0; + y[k] := 0.0; + end; + end; + + fac := (1.0 - sqr(AStrength/2)) / mo2; + if (fac < 0.0) then fac := 0.0; + y[2] := y[2] * fac; + + // Back-transformation to time domain + RealFT(y, mo2, -1); + + { wp: seems to invert the transformation above + for j := 1 to n do + y[j] := y[j] + rn1 * (y1 * (n - j) + yn * (j - 1)); + } + for j := 0 to n - 1 do + AData[j] := y[j+1]; + + y := nil; +end; + + +procedure FourierSmooth(var AData: DblDyneVec; AStrength, ANumProj: Integer); +var + n, numPts: Integer; + avg: DblDyneVec = nil; + i: Integer; +begin + numPts := Length(AData); + if numPts = 0 then + exit; + + n := numPts + ANumProj; + + SetLength(avg, n); + for i := 0 to numPts-1 do + avg[i] := AData[i]; + + if ANumProj > 0 then + for i := 0 to ANumProj-1 do + begin + avg[i] := AData[numPts - 1 - ANumProj + i]; + AData[i] := avg[i]; + end; + + DoFourierSmooth(avg, n, AStrength); + + // Replace original points with smoothed values + AData := VecCopy(avg); +end; + + +{ TSmoothingForm } + +constructor TSmoothingForm.Create(AOwner: TComponent); +begin + inherited; + + SetParams(DEFAULT_SMOOTHING_PARAMS); + + FChartFrame.SetXTitle('Cases'); + + FDataSeries := TLineSeries(FChartFrame.PlotXY(ptLinesAndSymbols, nil, nil, nil, nil, 'Original', DATA_COLORS[0])); + FDataSeries.Source := OriginalDataChartSource; + + FSmoothedSeries := TLineSeries(FChartFrame.PlotXY(ptLinesAndSymbols, nil, nil, nil, nil, 'Smoothed', DATA_COLORS[1])); + FSmoothedSeries.Source := SmoothedDataChartSource; + + FResidualSeries := TLineSeries(FChartFrame.PlotXY(ptLinesAndSymbols, nil, nil, nil, nil, 'Residual', DATA_COLORS[2])); + FResidualSeries.Source := ResidualsChartSource; + FResidualSeries.Active := false; +end; + + +procedure TSmoothingForm.AdjustConstraints; +var + w: Integer; +begin + w := Max(OKBtn.Width, CancelBtn.Width); + OKBtn.Constraints.MinWidth := w; + CancelBtn.Constraints.MinWidth := w; + + ParamsPanel.Constraints.MinWidth := Max( + 2*OKBtn.Width + OKBtn.BorderSpacing.Left, + ResidChk.Left + ResidChk.Width + ); +end; + + +procedure TSmoothingForm.AlphaEditChange(Sender: TObject); +var + x: Double; +begin + if not TryStrToFloat(AlphaEdit.Text, x) then + exit; + Compute; +end; + + +procedure TSmoothingForm.AlphaScrollChange(Sender: TObject); +begin + AlphaEdit.Value := AlphaScroll.Position / 100; +end; + + +procedure TSmoothingForm.Compute; +var + msg: String; + C: TWinControl; +begin + if (FSmoothedSeries = nil) or (FResidualSeries = nil) then + exit; + + if not Validate(msg, C) then + begin + C.SetFocus; + ErrorMsg(msg); + exit; + end; + + FParams.CenterOnMean := MeanChk.Checked; + FParams.DifferenceParams.Active := DataSmoothingChkLB.Checked[0]; + FParams.MovingAvgParams.Active := DataSmoothingChkLB.Checked[1]; + FParams.ExponentialParams.Active := DataSmoothingChkLB.Checked[2]; + FParams.FourierParams.Active := DataSmoothingChkLB.Checked[3]; + FParams.PolynomialParams.Active := DataSmoothingChkLB.Checked[4]; + + FSmoothedPoints := VecCopy(FOrigPoints); + + if FParams.DifferenceParams.Active then + begin + FParams.DifferenceParams.Order := DiffOrderEdit.Value; + FParams.DifferenceParams.MaxLag := LagEdit.Value; + DifferenceSmooth(FSmoothedPoints, DiffOrderEdit.Value, LagEdit.Value); + end; + + if FParams.MovingAvgParams.Active then + begin + FParams.MovingAvgParams.NormalizeWeights; + MovingAvgSmooth(FSmoothedPoints, FParams.MovingAvgParams.Weights, FParams.MovingAvgParams.Order, FNumProj); + end; + + if FParams.ExponentialParams.Active then + begin + FParams.ExponentialParams.Alpha := AlphaEdit.Value; + ExponentialSmooth(FSmoothedPoints, AlphaEdit.Value, FNumProj); + end; + + if FParams.FourierParams.Active then + begin + FParams.FourierParams.Strength := StrengthEdit.Value; + FourierSmooth(FSmoothedPoints, StrengthEdit.Value, FNumProj); + end; + + if FParams.PolynomialParams.Active then + begin + FParams.PolynomialParams.Order := PolyOrderEdit.Value; + PolynomialSmooth(FSmoothedPoints, FParams.PolynomialParams.Order); + end; + + FSmoothedSeries.Active := (FSmoothedPoints <> nil) and (Length(FSmoothedPoints) > 0); + SmoothedDataChartSource.PointsNumber := Length(FSmoothedPoints); + SmoothedDataChartSource.Reset; + + FResidualSeries.Active := ResidChk.Checked and FSmoothedSeries.Active; + if FResidualSeries.Active then + begin + ResidualsChartSource.PointsNumber := Length(FSmoothedPoints); + ResidualsChartSource.Reset; + end; +end; + + +procedure TSmoothingForm.DataSmoothingChkLBSelectionChange(Sender: TObject; + User: boolean); +var + index: Integer; +begin + index := DataSmoothingChkLB.ItemIndex; + Notebook.PageIndex := index; + NotebookGroup.Caption := DataSmoothingChkLB.Items[index] + ' Parameters'; + Notebook.Enabled := DataSmoothingChkLB.Checked[index]; + + Compute; +end; + + +procedure TSmoothingForm.DiffOrderEditChange(Sender: TObject); +begin + Compute; +end; + + +procedure TSmoothingForm.LagEditChange(Sender: TObject); +begin + Compute; +end; + + +procedure TSmoothingForm.MeanChkChange(Sender: TObject); +var + i: Integer; +begin + FParams.CenterOnMean := MeanChk.Checked; + if (Length(FOrigPoints) > 0) then + begin + if FParams.CenterOnMean then + for i := 0 to High(FOrigPoints) do + FOrigPoints[i] := FOrigPoints[i] - FMean + else + for i := 0 to High(FOrigPoints) do + FOrigPoints[i] := FOrigPoints[i] + FMean; + end; + OriginalDataChartSource.Reset; + Compute; +end; + + +procedure TSmoothingForm.PolyOrderEditChange(Sender: TObject); +begin + Compute; +end; + + +procedure TSmoothingForm.ResidChkChange(Sender: TObject); +begin + Compute; +end; + + +procedure TSmoothingForm.ResidualsChartSourceGetChartDataItem( + ASource: TUserDefinedChartSource; AIndex: Integer; var AItem: TChartDataItem); +begin + AItem.X := AIndex + 1; + AItem.Y := FSmoothedPoints[AIndex] - FOrigPoints[AIndex]; +end; + + +procedure TSmoothingForm.SmoothedDataChartSourceGetChartDataItem( + ASource: TUserDefinedChartSource; AIndex: Integer; var AItem: TChartDataItem); +begin + AItem.X := AIndex + 1; + AItem.Y := FSmoothedPoints[AIndex]; +end; + + +procedure TSmoothingForm.OriginalDataChartSourceGetChartDataItem( + ASource: TUserDefinedChartSource; AIndex: Integer; var AItem: TChartDataItem); +begin + AItem.X := AIndex + 1; + AItem.Y := FOrigPoints[AIndex]; +end; + + +procedure TSmoothingForm.SetParams(const AParams: TSmoothingParams); +begin + FParams := AParams; + + MeanChk.Checked := AParams.CenterOnMean; + + SetRawWeights(AParams.MovingAvgParams.RawWeights); + LagEdit.Value := FParams.DifferenceParams.MaxLag; + DiffOrderEdit.Value := FParams.DifferenceParams.Order; + AlphaEdit.Value := FParams.ExponentialParams.Alpha; + StrengthEdit.Value := FParams.FourierParams.Strength; + PolyOrderEdit.Value := FParams.PolynomialParams.Order; +end; + + +procedure TSmoothingForm.SetPoints(const APoints: DblDyneVec); +var + i: Integer; +begin + SetLength(FOrigPoints, Length(APoints)); + for i := 0 to High(APoints) do + FOrigPoints[i] := APoints[i]; + + FMean := VecMean(FOrigPoints); + if FParams.CenterOnMean then + for i := 0 to High(FOrigPoints) do + FOrigPoints[i] := FOrigPoints[i] - FMean; + + OriginaldataChartSource.PointsNumber := Length(FOrigPoints); + OriginalDataChartSource.Reset; + + Compute; +end; + + +procedure TSmoothingForm.SetRawWeights(const AWeights: DblDyneVec); +begin + FParams.MovingAvgParams.RawWeights := VecCopy(AWeights); + FParams.MovingAvgParams.NormalizeWeights; + WeightOrder.OnChange := nil; + WeightOrder.Value := Length(FParams.MovingAvgParams.RawWeights) - 1; + WeightOrder.OnChange := @WeightOrderChange; + WeightsToGrid; +end; + + +procedure TSmoothingForm.StrengthEditChange(Sender: TObject); +begin + Compute; +end; + + +procedure TSmoothingForm.StrengthEditEditingDone(Sender: TObject); +begin + Compute; +end; + + +function TSmoothingForm.Validate(out AMsg: String; + out AControl: TWinControl): boolean; +var + n, i: Integer; + x: Double; +begin + Result := false; + + // Difference smoothing + if DataSmoothingChkLB.Checked[0] then + begin + if (LagEdit.Text = '') or not TryStrToInt(LagEdit.Text, n) or (n < 0) then + begin + AMsg := 'Valid non-negative integer required.'; + AControl := LagEdit; + Notebook.PageIndex := 0; + exit; + end; + + if (DiffOrderEdit.Text = '') or not TryStrToInt(DiffOrderEdit.Text, n) or (n < 0) then + begin + AMsg := 'Valid non-negative integer required.'; + AControl := DiffOrderEdit; + Notebook.PageIndex := 0; + exit; + end; + end; + + // Moving average smoothing + if DataSmoothingChkLB.Checked[1] then + begin + if (WeightOrder.Text = '') or not TryStrToInt(WeightOrder.Text, n) or (n < 0) then + begin + AMsg := 'Valid non-negative integer required.'; + AControl := WeightOrder; + Notebook.PageIndex := 1; + exit; + end; + for i := 1 to WeightsGrid.RowCount-1 do + begin + if (WeightsGrid.Cells[1, i] = '') or not TryStrToFloat(WeightsGrid.Cells[1, i], x) then + begin + AMsg := 'Valid number required.'; + AControl := WeightsGrid; + Notebook.PageIndex := 1; + WeightsGrid.Row := i; + WeightsGrid.Col := 1; + exit; + end; + end; + end; + + // Exponential smoothing + if DataSmoothingChkLB.Checked[2] then + begin + if (AlphaEdit.Text = '') then + begin + AMsg := 'Floating point input between 0 and 1 required.'; + AControl := AlphaEdit; + Notebook.PageIndex := 2; + exit; + end; + if not TryStrToFloat(AlphaEdit.Text, x) or (x < 0) or (x > 1) then + begin + AMsg := 'Floating point input between 0 and 1 required.'; + AControl := AlphaEdit; + Notebook.PageIndex := 2; + exit; + end; + end; + + // Fourier smoothing + if DataSmoothingChkLB.Checked[3] then + begin + if (StrengthEdit.Text = '') or not TryStrToInt(StrengthEdit.Text, n) or (n < 0) then + begin + AMsg := 'Valid non-negative integer required.'; + AControl := StrengthEdit; + Notebook.PageIndex := 3; + exit; + end; + end; + + // Polynomial smoothing + if DataSmoothingChkLB.checked[4] then + begin + if (PolyOrderEdit.Text = '') or not TryStrToInt(PolyOrderEdit.Text, n) or (n < 0) then + begin + AMsg := 'Valid non-negative integer required.'; + AControl := PolyOrderEdit; + Notebook.PageIndex := 4; + exit; + end; + end; + + Result := true; +end; + + +procedure TSmoothingForm.WeightOrderChange(Sender: TObject); +begin + WeightOrderEditingDone(nil); +end; + + +procedure TSmoothingForm.WeightOrderEditingDone(Sender: TObject); +var + i: Integer; + n: Integer; +begin + n := Length(FParams.MovingAvgParams.RawWeights); + SetLength(FParams.MovingAvgParams.RawWeights, WeightOrder.Value + 1); + if Length(FParams.MovingAvgParams.RawWeights) > n then + for i := n to High(FParams.MovingAvgParams.RawWeights) do + FParams.MovingAvgParams.RawWeights[i] := 0; + FParams.MovingAvgParams.NormalizeWeights; + + WeightsToGrid; +end; + + +procedure TSmoothingForm.WeightsGridEditingDone(Sender: TObject); +var + x: Double; +begin + if WeightsGrid.Col <> 1 then + exit; + if TryStrToFloat(WeightsGrid.Cells[1, WeightsGrid.Row], x) then + FParams.MovingAvgParams.RawWeights[WeightsGrid.Row-1] := x + else + ErrorMsg('Not a valid number.'); + + WeightsToGrid; +end; + + +procedure TSmoothingForm.WeightsGridPrepareCanvas(Sender: TObject; aCol, + aRow: Integer; aState: TGridDrawState); +begin + if not Notebook.Enabled then + WeightsGrid.Canvas.Font.Color := clGray; +end; + + +procedure TSmoothingForm.WeightsGridSelectEditor(Sender: TObject; aCol, + aRow: Integer; var Editor: TWinControl); +begin + if (ACol = 2) then Editor := nil; +end; + + +procedure TSmoothingForm.WeightsToGrid; +var + i, r: Integer; +begin + FParams.MovingAvgParams.NormalizeWeights; + Weightsgrid.RowCount := Length(FParams.MovingAvgParams.RawWeights) + WeightsGrid.FixedRows; + + i := 0; + for r := 1 to WeightsGrid.RowCount-1 do + begin + if i = 0 then + WeightsGrid.Cells[0, r] := '0 (center)' + else + WeightsGrid.Cells[0, r] := IntToStr(i); + WeightsGrid.Cells[1, r] := Format('%.2g', [FParams.MovingAvgParams.RawWeights[i]]); + WeightsGrid.Cells[2, r] := Format('%.3f', [FParams.MovingAvgParams.Weights[i]]); + inc(i); + end; + + Compute; +end; + + +end. + diff --git a/applications/lazstats/source/units/mathunit.pas b/applications/lazstats/source/units/mathunit.pas index 51cbd7691..47f84d97e 100644 --- a/applications/lazstats/source/units/mathunit.pas +++ b/applications/lazstats/source/units/mathunit.pas @@ -400,8 +400,10 @@ end; { Cumulative probability of the F distribution for given F value and with DF1 and DF2 degrees of freedom } function ProbF(F, DF1, DF2: Double): Double; +{ var b1, b2: Double; +} begin if f <= 0.0 then Result := 1.0 diff --git a/applications/lazstats/source/units/matrixunit.pas b/applications/lazstats/source/units/matrixunit.pas index fb3703adb..fff762dca 100644 --- a/applications/lazstats/source/units/matrixunit.pas +++ b/applications/lazstats/source/units/matrixunit.pas @@ -11,15 +11,15 @@ uses Globals; type - TDblVector = DblDyneVec; - TDblMatrix = DblDyneMat; - TDblCube = DblDyneCube; - TDblQuad = DblDyneQuad; + TDblVector = DblDyneVec; // array of double + TDblMatrix = DblDyneMat; // array of array of double + TDblCube = DblDyneCube; // array of array of array of double + TDblQuad = DblDyneQuad; // array of array of array of array of double - TIntVector = IntDyneVec; - TIntMatrix = IntDyneMat; - TIntCube = IntDyneCube; - TIntQuad = IntDyneQuad; + TIntVector = IntDyneVec; // array of integer + TIntMatrix = IntDyneMat; // array of array of integer + TIntCube = IntDyneCube; // array of array of array of integer + TIntQuad = IntDyneQuad; // array of array of array of array of integer EMatrix = class(Exception); @@ -136,6 +136,8 @@ implementation uses Math; +const + TWO_PI = 2.0 * pi; // redeclaration to avoid dependence on MathUnit operator + (A, B: TDblVector): TDblVector; var diff --git a/applications/lazstats/source/units/regressionunit.pas b/applications/lazstats/source/units/regressionunit.pas index ae715ab2f..0b37b46dd 100644 --- a/applications/lazstats/source/units/regressionunit.pas +++ b/applications/lazstats/source/units/regressionunit.pas @@ -178,14 +178,14 @@ begin // Product of transposed augmented with augmented matrix XT_X := XT * X; - AResults.XT_X := Matcopy(XT_X); + AResults.XT_X := MatCopy(XT_X); // Product of inverse of XT_X with XT_Y yields the coefficients. Intercept is last. inv_XT_X := MatInverse(XT_X); AResults.Coeffs := inv_XT_X * XT_Y; // THe XT_X matrix contains the sum of squares (diagonal elements), sum of - // cross-products (off-diagonal elements and the sum the variable values + // cross-products (off-diagonal elements and the sum of the variable values // (augmented column). SetLength(AResults.Means, mx+1); SetLength(AResults.Variances, mx+1); @@ -206,7 +206,10 @@ begin for j := 0 to mx-1 do begin VarCovar[i, j] := (XT_X[i, j] - XT_X[i, mx] * XT_X[j, mx] / nx) / (nx - 1); - Correlations[i, j] := VarCovar[i, j] / (StdDevs[i] * StdDevs[j]); + if StdDevs[i] <> 0 then + Correlations[i, j] := VarCovar[i, j] / (StdDevs[i] * StdDevs[j]) + else + Correlations[i, j] := 0.0; end; end;