Delphi版的线性回归分析

文章《C语言版的线性回归分析函数》发布后，不少朋友留言或给我来信，询问能否提供Delphi版的线性回归分析代码，因C语言版是我以前DOS下的老代码稍作整理后发布的，所以没有现成的Delphi代码，今天比较闲，于是将C代码改写为Delphi代码贴在下面，有关的回归公式说明及例子图示可参见《C语言版的线性回归分析函数》，这里不再累叙，由于改写时间仓促，可能有错误，请不吝指出，亦可来信建议：maozefa@hotmail.com

线性回归分析代码：

unitRegression;

interface

usesSysUtils;

type
PEquationsData = ^ TEquationsData;
TEquationsData = array[ 0 .. 0 ]ofDouble;
// 线性回归
TLinearRegression = class (TObject)
private
FData:PEquationsData;
FAnswer:PEquationsData;
FSquareSum:Double;
FSurplusSum:Double;
FRowCount:Integer;
FColCount:Integer;
FModify:Boolean;
functionGetAnswer(Index:Integer):Double;
functionGetItem(ARow,ACol:Integer):Double;
procedureSetItem(ARow,ACol:Integer; const Value:Double);
procedureSetColCount( const Value:Integer);
procedureSetRowCount( const Value:Integer);
procedureSetSize( const ARowCount,AColCount:Integer);
procedureSetModify( const Value:Boolean);
functionGetCorrelation:Double;
functionGetDeviatSum:Double;
functionGetFTest:Double;
functionGetSurplus:Double;
functionGetVariance:Double;
functionGetStandardDiffer:Double;
functionGetEstimate(ARow:Integer):Double;
public
constructorCreate( const AData; const ARowCount,AColCount:Integer);overload;
destructorDestroy; override ;
// 计算回归方程
procedureCalculation;
// 设置回归数据
// AData[ARowCount*AColCount]二维数组；X1i,X2i,...Xni,Yi(i=0toARowCount-1)
// ARowCount：数据行数；AColCount数据列数
procedureSetData( const AData; const ARowCount,AColCount:Integer);
// 数据列数(自变量个数+Y)
propertyColCount:IntegerreadFColCountwriteSetColCount;
// 数据行数
propertyRowCount:IntegerreadFRowCountwriteSetRowCount;
// 原始数据
propertyData[ARow,ACol:Integer]:DoublereadGetItemwriteSetItem; default ;
propertyModify:BooleanreadFModify;
// 回归系数数组(B0,B1...Bn)
propertyAnswer[Index:Integer]:DoublereadGetAnswer;
// Y估计值
propertyEstimate[ARow:Integer]:DoublereadGetEstimate;
// 回归平方和
propertyRegresSquareSum:DoublereadFSquareSum;
// 剩余平方和
propertySurplusSquareSum:DoublereadFSurplusSum;
// 离差平方和
propertyDeviatSquareSum:DoublereadGetDeviatSum;
// 回归方差
propertyRegresVariance:DoublereadGetVariance;
// 剩余方差
propertySurplusVariance:DoublereadGetSurplus;
// 标准误差
propertyStandardDiffer:DoublereadGetStandardDiffer;
// 相关系数
propertyCorrelation:DoublereadGetCorrelation;
// F检验
propertyF_Test:DoublereadGetFTest;
end;

// 解线性方程。AData[count*(count+1)]矩阵数组；count：方程元数；
// Answer[count]：求解数组。返回：True求解成功，否则无解或者无穷解
functionLinearEquations( const AData;Count:Integer;varAnswer:arrayofDouble):Boolean;

implementation

const
SMatrixSizeError = ' Regressiondatamatrixcannotbelessthan2*2 ' ;
SIndexOutOfRange = ' indexoutofrange ' ;
SEquationNoSolution = ' EquationnosolutionorInfiniteSolutions ' ;

functionLinearEquations( const AData;Count:Integer;varAnswer:arrayofDouble):Boolean;
var
j,m,n,ColCount:Integer;
tmp:Double;
Data,d:PEquationsData;
begin
Result: = False;
if Count < 2 thenExit;

ColCount: = Count + 1 ;
GetMem(Data,Count * ColCount * Sizeof(Double));
GetMem(d,ColCount * Sizeof(Double));
try
Move(AData,Data ^ ,Count * ColCount * Sizeof(Double));
for m: = 0 toCount - 2 do
begin
n: = m + 1 ;
// 如果主对角线元素为0，行交换
while (n < Count)and(Data ^ [m * ColCount + m] = 0.0 ) do
begin
if Data ^ [n * ColCount + m] <> 0.0 then
begin
Move(Data ^ [m * ColCount + m],d ^ ,ColCount * Sizeof(Double));
Move(Data ^ [n * ColCount + m],Data ^ [m * ColCount + m],ColCount * Sizeof(Double));
Move(d ^ ,Data ^ [n * ColCount + m],ColCount * Sizeof(Double));
end;
Inc(n);
end;
// 行交换后，主对角线元素仍然为0，无解
if Data ^ [m * ColCount + m] = 0.0 thenExit;
// 消元
for n: = m + 1 toCount - 1 do
begin
tmp: = Data ^ [n * ColCount + m] / Data ^ [m * ColCount + m];
for j: = mtoCount do
Data ^ [n * ColCount + j]: = Data ^ [n * ColCount + j] - tmp * Data ^ [m * ColCount + j];
end;
end;
FillChar(d ^ ,Count * Sizeof(Double), 0 );
// 求得count-1的元
Answer[Count - 1 ]: = Data ^ [(Count - 1 ) * ColCount + Count] /
Data ^ [(Count - 1 ) * ColCount + Count - 1 ];
// 逐行代入求各元
for m: = Count - 2 downto 0 do
begin
for j: = Count - 1 downtom + 1 do
d ^ [m]: = d ^ [m] + Answer[j] * Data ^ [m * ColCount + j];
Answer[m]: = (Data ^ [m * ColCount + Count] - d ^ [m]) / Data ^ [m * ColCount + m];
end;
Result: = True;
finally
FreeMem(d);
FreeMem(Data);
end;
end;

{TLinearRegression}

procedureTLinearRegression.Calculation;
var
m,n,i,count:Integer;
dat:PEquationsData;
a,b,d:Double;
begin
if (FRowCount < 2 )or(FColCount < 2 )then
raiseException.Create(SMatrixSizeError);
if notFModifythenExit;
GetMem(dat,FColCount * (FColCount + 1 ) * Sizeof(Double));
try
count: = FColCount - 1 ;
dat ^ [ 0 ]: = FRowCount;
for n: = 0 tocount - 1 do
begin
a: = 0.0 ;
b: = 0.0 ;
for m: = 0 toFRowCount - 1 do
begin
d: = FData ^ [m * FColCount + n];
a: = a + d;
b: = b + d * d;
end;
dat ^ [n + 1 ]: = a;
dat ^ [(n + 1 ) * (FColCount + 1 )]: = a;
dat ^ [(n + 1 ) * (FColCount + 1 ) + n + 1 ]: = b;
for i: = n + 1 tocount - 1 do
begin
a: = 0.0 ;
for m: = 0 toFRowCount - 1 do
a: = a + FData ^ [m * FColCount + n] * FData ^ [m * FColCount + i];
dat ^ [(n + 1 ) * (FColCount + 1 ) + i + 1 ]: = a;
dat ^ [(i + 1 ) * (FColCount + 1 ) + n + 1 ]: = a;
end;
end;
b: = 0.0 ;
for m: = 0 toFRowCount - 1 do
b: = b + FData ^ [m * FColCount + count];
dat ^ [FColCount]: = b;
for n: = 0 tocount - 1 do
begin
a: = 0.0 ;
for m: = 0 toFRowCount - 1 do
a: = a + FData ^ [m * FColCount + n] * FData ^ [m * FColCount + count];
dat ^ [(n + 1 ) * (FColCount + 1 ) + FColCount]: = a;
end;
if notLinearEquations(dat ^ ,FColCount,FAnswer ^ )then
raiseException.Create(SEquationNoSolution);
FSquareSum: = 0.0 ;
FSurplusSum: = 0.0 ;
b: = b / FRowCount;
for m: = 0 toFRowCount - 1 do
begin
a: = FAnswer ^ [ 0 ];
for i: = 1 tocount do
a: = a + FData ^ [m * FColCount + i - 1 ] * FAnswer[i];
FSquareSum: = FSquareSum + (a - b) * (a - b);
d: = FData ^ [m * FColCount + count];
FSurplusSum: = FSurplusSum + (d - a) * (d - a);
end;
SetModify(False);
finally
FreeMem(dat);
end;
end;

constructorTLinearRegression.Create( const AData; const ARowCount,
AColCount:Integer);
begin
SetData(AData,ARowCount,AColCount);
end;

destructorTLinearRegression.Destroy;
begin
SetSize( 0 , 0 );
end;

functionTLinearRegression.GetAnswer(Index:Integer):Double;
begin
if (Index < 0 )or(Index >= FColCount)then
raiseException.Create(SIndexOutOfRange);
if notAssigned(FAnswer)then
Result: = 0.0
else
Result: = FAnswer ^ [Index];
end;

functionTLinearRegression.GetCorrelation:Double;
begin
Result: = DeviatSquareSum;
if Result <> 0.0 then
Result: = Sqrt(FSquareSum / Result);
end;

functionTLinearRegression.GetDeviatSum:Double;
begin
Result: = FSquareSum + FSurplusSum;
end;

functionTLinearRegression.GetEstimate(ARow:Integer):Double;
var
I:Integer;
begin
if (ARow < 0 )or(ARow >= FRowCount)then
raiseException.Create(SIndexOutOfRange);
Result: = Answer[ 0 ];
for I: = 1 toColCount - 1 do
Result: = Result + FData ^ [ARow * FColCount + I - 1 ] * Answer[I];
end;

functionTLinearRegression.GetFTest:Double;
begin
Result: = SurplusVariance;
if Result <> 0.0 then
Result: = RegresVariance / Result;
end;

functionTLinearRegression.GetItem(ARow,ACol:Integer):Double;
begin
if (ARow < 0 )or(ARow >= FRowCount)or(ACol < 0 )or(ACol >= FColCount)then
raiseException.Create(SIndexOutOfRange);
Result: = FData ^ [ARow * FColCount + ACol];
end;

functionTLinearRegression.GetStandardDiffer:Double;
begin
Result: = Sqrt(SurplusVariance);
end;

functionTLinearRegression.GetSurplus:Double;
begin
if FRowCount - FColCount < 1 then
Result: = 0.0
else
Result: = FSurplusSum / (FRowCount - FColCount);
end;

functionTLinearRegression.GetVariance:Double;
begin
if FColCount < 2 then
Result: = 0.0
else
Result: = FSquareSum / (FColCount - 1 );
end;

procedureTLinearRegression.SetColCount( const Value:Integer);
begin
if Value < 2 then
raiseException.Create(SMatrixSizeError);
SetSize(FRowCount,Value);
end;

procedureTLinearRegression.SetData( const AData; const ARowCount,AColCount:Integer);
begin
if (ARowCount < 2 )or(AColCount < 2 )then
raiseException.Create(SMatrixSizeError);
SetSize(ARowCount,AColCount);
Move(AData,FData ^ ,FRowCount * FColCount * Sizeof(Double));
end;

procedureTLinearRegression.SetItem(ARow,ACol:Integer; const Value:Double);
begin
if (ARow < 0 )or(ARow >= FRowCount)or(ACol < 0 )or(ACol >= FColCount)then
raiseException.Create(SIndexOutOfRange);
if FData ^ [ARow * (FColCount) + ACol] <> Valuethen
begin
FData ^ [ARow * (FColCount) + ACol]: = Value;
SetModify(True);
end;
end;

procedureTLinearRegression.SetModify( const Value:Boolean);
begin
if FModify <> Valuethen
begin
FModify: = Value;
if FModifythen
begin
FillChar(FAnswer ^ ,FColCount * Sizeof(Double), 0 );
FSquareSum: = 0.0 ;
FSurplusSum: = 0.0 ;
end;
end;
end;

procedureTLinearRegression.SetRowCount( const Value:Integer);
begin
if Value < 2 then
raiseException.Create(SMatrixSizeError);
SetSize(Value,FColCount);
end;

procedureTLinearRegression.SetSize( const ARowCount,AColCount:Integer);
begin
if (FRowCount = ARowCount)and(FColCount = AColCount)then
Exit;
if Assigned(FData)then
begin
FreeMem(FData);
FData: = nil;
FreeMem(FAnswer);
FAnswer: = nil;
FModify: = False;
end;

FRowCount: = ARowCount;
FColCount: = AColCount;

if (FRowCount = 0 )or(FColCount = 0 )thenExit;

GetMem(FData,FRowCount * FColCount * Sizeof(Double));
FillChar(FData ^ ,FRowCount * FColCount * Sizeof(Double), 0 );
GetMem(FAnswer,FColCount * Sizeof(Double));
SetModify(True);
end;

end.

因为一元线性回归分析本是多元线性回归分析的一个特例，因此原C代码中的一元线性回归函数取消，一元线性回归和多元线性回归都使用TLinearRegression类。下面是Pascal控制台应用程序例子：

programLinearRegression;

{$APPTYPECONSOLE}

uses
SysUtils,
Regression in ' ....pasRegression.pas ' ;

const
data1:array[ 1 .. 12 , 1 .. 2 ]ofDouble = (
// XY
( 187.1 , 25.4 ),
( 179.5 , 22.8 ),
( 157.0 , 20.6 ),
( 197.0 , 21.8 ),
( 239.4 , 32.4 ),
( 217.8 , 24.4 ),
( 227.1 , 29.3 ),
( 233.4 , 27.9 ),
( 242.0 , 27.8 ),
( 251.9 , 34.2 ),
( 230.0 , 29.2 ),
( 271.8 , 30.0 )
);

data:array[ 1 .. 15 , 1 .. 5 ]ofDouble = (
// X1X2X3X4Y
( 316 , 1536 , 874 , 981 , 3894 ),
( 385 , 1771 , 777 , 1386 , 4628 ),
( 299 , 1565 , 678 , 1672 , 4569 ),
( 326 , 1970 , 785 , 1864 , 5340 ),
( 441 , 1890 , 785 , 2143 , 5449 ),
( 460 , 2050 , 709 , 2176 , 5599 ),
( 470 , 1873 , 673 , 1769 , 5010 ),
( 504 , 1955 , 793 , 2207 , 5694 ),
( 348 , 2016 , 968 , 2251 , 5792 ),
( 400 , 2199 , 944 , 2390 , 6126 ),
( 496 , 1328 , 749 , 2287 , 5025 ),
( 497 , 1920 , 952 , 2388 , 5924 ),
( 533 , 1400 , 1452 , 2093 , 5657 ),
( 506 , 1612 , 1587 , 2083 , 6019 ),
( 458 , 1613 , 1485 , 2390 , 6141 )
);

procedureDisplay(s: string ;R:TLinearRegression);
var
i:Integer;
v,o:Double;
begin
Writeln(s);
Writeln( ' 回归方程式: ' );
Write( ' Y= ' ,R.Answer[ 0 ]: 1 : 5 );
for i: = 1 toR.ColCount - 1 do
Write( ' + ' ,R.Answer[i]: 1 : 5 , ' *X ' ,i);
Writeln;
Writeln( ' 回归显著性检验: ' );
Writeln( ' 回归平方和： ' ,R.RegresSquareSum: 12 : 4 , ' 回归方差： ' ,R.RegresVariance: 12 : 4 );
Writeln( ' 剩余平方和： '