lazarus-ccr/applications/lazstats/source/units/regressionunit.pas

unit RegressionUnit;

{$mode objfpc}{$H+}
{$modeswitch advancedrecords}

interface

uses
  Classes, SysUtils,
  MatrixUnit;

type
  TBivariateRegressionResults = record
    Slope, Intercept: Double;
    XMean, YMean: Double;
    XVariance, YVariance: Double;
    XStdDev, YStdDev: Double;
    StdErrorPredicted: Double;
    R: Double;
    t: Double;
    SXX, SXY, SYY: Double;
    Count, DF: Integer;
    function ConfidenceLimits(x: Double; Upper: Boolean): Double;
  end;

procedure BivariateRegression(const xData, yData: TDblVector;
  AConfLevel: Double; out AResults: TBivariateRegressionResults);

type
  TMultipleRegressionResults = record
  end;

procedure MultipleRegression(const xData: TDblMatrix; const yData: TDblVector;
  AConfLevel: Double; out AResults: TMultipleRegressionResults);


implementation

uses
  MathUnit;

{===============================================================================
                        Bivariate regression

  Calculates a bivariate regression y = a x + b

  x is a vector of independent observations, y of dependent observations.
  Both vectors must have the same length.

  It is assumed that xData and yData contain at least 3 elements and
  have the same count of elements.
===============================================================================}
procedure BivariateRegression(const xData, yData: TDblVector;
  AConfLevel: Double; out AResults: TBivariateRegressionResults);
var
  i: Integer;
begin
  with AResults do
  begin
    Count := Length(xData);

    // Calculate means, variances, stddevs
    VecMeanVarStdDevSS(xData, XMean, XVariance, XStdDev, SXX);
    VecMeanVarStdDevSS(yData, YMean, YVariance, YStdDev, SYY);

    SXY := 0;
    for i := 0 to Count-1 do
      SXY := SXY + xData[i] * yData[i];

    R := (SXY - XMean * YMean * Count) / ((Count - 1) * XStdDev * YStdDev);
    StdErrorPredicted := sqrt(1.0 - sqr(R)) * YStdDev * sqrt((Count - 1) / (Count - 2));

    Slope := R * YStdDev / XStdDev;
    Intercept := YMean - Slope * XMean;

    DF := Count - 2;
    t := InverseT(AConfLevel, DF);
  end;
end;


function TBivariateRegressionResults.ConfidenceLimits(x: Double; Upper: Boolean): Double;
var
  yPred, seData: Double;
begin
  yPred := Intercept + Slope * x;
  seData := StdErrorPredicted * sqrt(1.0 + 1/Count + sqr(x - XMean)/SXX);
  if Upper then
    Result := yPred + t*seData
  else
    Result := yPred - t*seData;
end;


{===============================================================================
                              Multiple regression
===============================================================================}
procedure MultipleRegression(const xData: TDblMatrix; const yData: TDblVector;
  AConfLevel: Double; out AResults: TMultipleRegressionResults);
begin
end;



end.