最近有个需求,要拟合一些数据采集通道的修正结果,由于数据量大,需要自己写软件处理,于是上网找了个现成的..bug较多,基本发现的bug都修复了,增加了一些新的功能.
只适用于 y = m0 + m1*x + m2*x^2 +....+mn*x^n 这一个函数的拟合.阶数可以自行定义,也可以根据判断最接近1的R来找到最合适的阶数.最大阶数为20.
一样会送上demo,以及全部源码.demo里的数据为r141b(制冷剂)的饱和压力温度.7次方的公式就可以比较完美的拟合了.
功能实现unit的pas,原作者的版权声明基本都在...修改者-我的声明也在.嘿嘿,貌似一次上传不完,分段了:
1 unit CYHSMatrixUnit;
2
3 interface
4
5 uses
6 Classes, Math, SysUtils;
7 {
8 TCMatrix：是为了用最小二乘法拟和曲线而开发的一个类，所以在这里的数组必须满足最小二乘法
9 的要求，x,y数组必须维数相同。
10 作者：yanghai0437 QQ:955743
11 时间：2003年5月19日
12 修改: solokey www.solokey.cn
13 时间: 2010年12月17日
14 本类中使用的数组下标均是0序的，也就是数组的第一个值的下标是0，第二个值的下标是1。
15 }
16 type
17 TDoubleArray = array of Double;
18
19 type
20 TCMatrix = class
21 private
22 Fdou_MaxtrixValue: Double; //用于返回行列式的值
23 FbIsNeedSort: Boolean; //是否需要排序
24 FnCoordinateSize: Integer; //坐标个数
25 FnMaxPower: Integer; //最高次幂
26 FxCoordinate: array of Double; //X坐标数组
27 FyCoordinate: array of Double; //Y坐标数组
28 FSumxCoordinate: array of Double; //最小二乘法中的X^n的和的值(n=0,...m+m) m是方程的次数
29 FSumyCoordinate: array of Double; //最小二乘法中的X^n*Y的和的值(n=0,...m) m是方程的次数
30 FLuSuccess: Boolean;
31 FEquationPara: TDoubleArray;
32 procedure Muti(MaxPower: Integer);
33 function SetXYCoordinate(pXCoordinate, pYCoordinate: array of Double): Boolean; //最小二乘法
34 public
35 constructor Create;
36 destructor Destory;
37 property LuSuccess: Boolean read FLuSuccess; //true表示进行Lu分解成功
38 property EquationPara:TDoubleArray read FEquationPara; //求出来的方程的系数 * FEquationPara[0]是常数项 * FEquationPara[1]是一次项系数...
39 function GetMaxPower: Integer; //得到方程的最高次幂
40 function GetCoordinateSize: Integer; //返回x,y坐标个数
41 function GetMatrixValue: Double; //返回行列式的值
42 function BeginComputer(const X, Y: array of Double; const MaxPower: Integer = 1;const GetBestResult: Boolean = False;const Sort: Boolean = False): Boolean;
43 class procedure SortArray(var pSort: array of Double); //为数组排序--冒泡排序法
44 class procedure SortTwoArray(var pSortx: array of Double; var pSortY: array of Double); //根据X坐标排序--冒泡排序法
45 class function MatrixLU(var Matrixs, SumyCoordinate, EquationParam: array of Double;
46 i_Power: integer; IsOnlyLu: Boolean; var MaxtrixValue: Extended): Boolean; //行列式 L U分解
47 class function CalY1(X: Double; const EquationParam: array of Double): Double;
48 class function CalY1Array(const X, EquationParam: array of Double; var Y1: array of Double): Boolean;
49 class function CalR(const X, Y, EquationParam: array of Double): Double; //相关系数
50 class function CalRMSE(const X, Y, EquationParam: array of Double): Double; //均方差
51 end;
52
53 const
54 MaxMi = 20; //用最小二乘法时是最高次幂 ；只求行列式的值时是最大阶数
55
56 implementation
57 //------------------------------------------------------------------------------
58 //功能：构造函数，初始化默认值
59 //参数：无
60 //返回值：无
61
62 constructor TCMatrix.Create;
63 begin
64 FLuSuccess := False;
65 FnMaxPower := 1;
66 FbIsNeedSort := False;
67 Fdou_MaxtrixValue := 0;
68 FnCoordinateSize := 0;
69 inherited;
70 end;
71
72 //------------------------------------------------------------------------------
73 //功能：得到在这次使用过程中的方程的最高次幂
74 //参数：无
75 //返回值：方程的最高次幂
76 //
77
78 function TCMatrix.GetMaxPower: Integer;
79 begin
80 Result := Length(FEquationPara) - 1;
81 end;
82
83 //------------------------------------------------------------------------------
84 //功能：制作坐标数组的副本。为这个类设置坐标数组，主要是为了在最小二乘法使用
85 //参数：pXCoordinate x 坐标，也是最小二乘法中的因变量，pYCoordinate ，y坐标
86 //返回值：如果制作副本成功返回true,否则返回false,
87 //如果两个数组的维数不同则不能进行最小二乘法，所以就会返回false
88
89 function TCMatrix.SetXYCoordinate(pXCoordinate: array of Double; pYCoordinate: array of Double): Boolean;//设置X,y坐标数组
90 var
91 nXCount, nYCount, i, nLowX, nLowY: Integer;
92 begin
93 Result := False;
94 //先释放坐标数组
95 SetLength(FxCoordinate, 0); //X坐标数组
96 SetLength(FyCoordinate, 0); //Y坐标数组
97 nXCount := Length(pXCoordinate);
98 nYCount := Length(pYCoordinate);
99 nLowX := Low(pXCoordinate);
100 nLowY := Low(pYCoordinate);
101 //判断是否有数据
102 if ((nXCount = 0) or (nXCount <> nYCount) or (nLowX <> nLowY)) then
103 Exit;
104 FnCoordinateSize := nXCount;
105 SetLength(FxCoordinate, nXCount); //X坐标数组
106 SetLength(FyCoordinate, nXCount); //Y坐标数组
107 for i := 0 to nXCount - 1 do
108 begin
109 FxCoordinate[i] := pXCoordinate[i + Low(pXCoordinate)];
110 FYCoordinate[i] := pYCoordinate[i + Low(pXCoordinate)];
111 end;
112 Result := True;
113 end;
114 //------------------------------------------------------------------------------
115 //功能：返回x,y坐标个数
116 //参数：
117 //返回值：返回x,y坐标个数
118
119 function TCMatrix.GetCoordinateSize: Integer;
120 begin
121 Result := FnCoordinateSize;
122 end;
123 //------------------------------------------------------------------------------
124 //功能：为double数组排序(升序)--冒泡排序法
125 //参数：
126 //返回值：
127
128 class procedure TCMatrix.SortArray(var pSort: array of Double);
129 var
130 nCount, nLow, i, j: Integer;
131 d_temp: Double; //临时变量
132 begin
133 nCount := High(pSort);
134 nLow := Low(pSort);
135 //数组中没有数据，返回
136 if (nCount = 0) then
137 Exit;
138 for i := nLow to nCount - 1 do
139 begin
140 for j := i + 1 to nCount do
141 begin
142 if (pSort[i] > pSort[j]) then
143 begin
144 d_temp := pSort[j];
145 pSort[j] := pSort[i];
146 pSort[i] := d_temp;
147 end;
148 end; //end for j
149 end; // end for i
150 end;
1
2 //------------------------------------------------------------------------------
3 //功能：为double数组排序。在坐标系中的排序(升序)，(x,y)代表一个点，所以y数组坐标将跟着
4 //x数组坐标变化,而--冒泡排序法
5 //参数： pSortx 坐标数组 pSortY 坐标数组
6 //返回值：
7 class procedure TCMatrix.SortTwoArray(var pSortx: array of Double; var pSortY: array of Double);
8 var
9 nCount, nLow, i, j: Integer;
10 d_tempx, d_tempy: Double; //临时变量
11 begin//根据X坐标排序--冒泡排序法
12 nCount := High(pSortx);
13 nLow := Low(pSortx);
14 i := High(pSortY);
15 j := Low(pSortY);
16 //数组中没有数据，或者两者的维数不相等，或者两者的初始维数序号不等，则不能排序，返回
17 if ((nCount = 0) or (i <> nCount) or (j <> nLow)) then
18 Exit;
19 for i := nLow to nCount - 1 do
20 begin
21 for j := i + 1 to nCount do
22 begin
23 if (pSortx[i] > pSortx[j]) then
24 begin //Y坐标根据X坐标变化
25 d_tempx := pSortx[j];
26 pSortx[j] := pSortx[i];
27 pSortx[i] := d_tempx;
28 d_tempy := pSorty[j];
29 pSorty[j] := pSorty[i];
30 pSorty[i] := d_tempy;
31 end;
32 end;
33 end;
34 end;
35 //------------------------------------------------------------------------------
36 //功能：最小二乘法
37 //参数： 具体请参考计算方法书中的最小二乘法章节
38 //返回值：将计算出来的值存入FSumxCoordinate ，FSumyCoordinate两个数组中
39 //
40
41 procedure TCMatrix.Muti(MaxPower: Integer);
42 var
43 dou_TempSum, dou_m: Double;
44 i, N, k: Integer;
45 TempX: array of Double;
46 begin
47 try
48 //重新分配数组空间
49 SetLength(FSumxCoordinate, 0); //最小二乘法中的X^n的和的值(n=0,...m+m) m是方程的次数
50 SetLength(FSumyCoordinate, 0); //最小二乘法中的X^n*Y的和的值(n=0,...m) m是方程的次数
51 N := High(FxCoordinate) + 1;
52 SetLength(FSumxCoordinate, (MaxPower + 1) * 2); //得到坐标的个数
53 for i := 0 to (MaxPower + 1) * 2 - 1 do
54 begin
55 dou_TempSum := 0;
56 SetLength(TempX, 0);
57 SetLength(TempX, N);
58 for k := 0 to N - 1 do
59 begin
60 dou_m := Power(FxCoordinate[k], i); //求X的i次方
61 TempX[k] := dou_m;
62 end;
63 dou_TempSum := dou_TempSum + Sum(TempX); //求X的i次方的和
64 //最小二乘法中的X^n的和的值(n=0,...m+m) m是方程的次数
65 FSumxCoordinate[i] := dou_TempSum;
66 end;
67 SetLength(FSumyCoordinate, MaxPower + 1);
68 for i := 0 to MaxPower do
69 begin
70 dou_TempSum := 0;
71 SetLength(TempX, 0);
72 SetLength(TempX, N);
73 for k := 0 to N - 1 do
74 begin
75 dou_m := Power(FxCoordinate[k], i); //求X(k)的i次方
76 dou_m := dou_m * FyCoordinate[k]; //求X(k)的i次方 乘以 y(k)
77 TempX[k] := dou_m;
78 end;
79 dou_TempSum := dou_TempSum + Sum(TempX); //求X(k)的i次方 乘以 y(k) 的和
80 //最小二乘法中的X^n*Y的和的值(n=0,...m) m是方程的次数
81 FSumyCoordinate[i] := dou_TempSum;
82 end;
83 finally
84 SetLength(TempX, 0);
85 TempX := nil;
86 end;
87 end;
1 //------------------------------------------------------------------------------
2 //功能：行列式 L U分解
3 //参数： Matrixs 矩阵行列式数组，i_Power 阶数－1，
4 //IsOnlyLu False 表示还要进行最小二乘法计算，得到拟和曲线系数,true表示只求行列式的值
5 //返回值： L U分解成功返回true,否则返回false
6 //
7 class function TCMatrix.MatrixLU(var Matrixs, SumyCoordinate, EquationParam: array of Double;
8 i_Power: integer; IsOnlyLu: Boolean; var MaxtrixValue: Extended): Boolean;
9 var
10 dou_Result, dou_temp: Extended;
11 dou_Matrixs: array[0..MaxMi, 0..MaxMi] of Extended; //定义二维数组 用于存储 转换后的行列式的值
12 i_Plusminus, j, i, i_x: integer;
13 L: array[0..MaxMi, 0..MaxMi] of Extended;
14 U: array[0..MaxMi, 0..MaxMi] of Extended;
15 dou_Y: array[0..MaxMi] of Extended;
16 Temp, TempB: Extended;
17 k, m, N, CC: Integer;
18 BB, Rows: Integer;
19 label
20 SS1, SS2, SS3;
21 begin
22 //*Matrixs 行列式数组
23 //i_Power 方程的次数 也就是行列式的阶-1，
24 //IsOnlyLu False 表示还要进行最小二乘法计算
25 Result := False;
26 MaxtrixValue := 0;
27 if i_Power > MaxMi then //方程的次数高于指定的最高次数
28 Exit;
29 dou_Result := 0; //要返回的值
30 i_Plusminus := 1; //计算时如果行列式进行了列与列的交换，则结果变负
31
32 //注意：行列式数组的下标是从0开始的
33 if not IsOnlyLu then
34 begin
35 for i := 0 to i_Power do
36 begin
37 for j := 0 to i_Power do
38 dou_Matrixs[i][j] := Matrixs[j + i];
39 end;
40 end
41 else//只用于与求行列式的值时，把行列式直接付值给数组
42 begin
43 i_x := (i_Power + 1) * (i_Power + 1);
44 //行数和列数不相等则不能计算行列式的值
45 i := High(Matrixs) + 1;
46 if (i_x <> i) then
47 begin
48 MaxtrixValue := dou_Result;
49 Exit;
50 end;
51 //把数组转换成矩阵行列式，矩阵行列式数组的下标是从0开始的
52 for i := 0 to i_Power do
53 begin
54 for j := 0 to i_Power do
55 dou_Matrixs[i][j] := Matrixs[j + i * (i_Power + 1)];
56 end;
57 end;
58 //判断是否需要换列--当第一行第一列=0时会产生除0错误，所以需要交换
59 if (abs(dou_Matrixs[0][0]) <= 0.000001) then
60 begin //需要交换
61 j := 0;
62 for i := 1 to i_Power do
63 begin
64 if abs(dou_Matrixs[0][i]) >= 0.000001 then //只在第一行中找不=0的值
65 begin
66 j := i; //找到第一行不为0的列
67 Break;
68 end;
69 end;
70 if j <> 0 then //找到了不=0的列的值
71 begin
72 i_Plusminus := i_Plusminus * (-1);
73 dou_temp := 0;
74 for i := 0 to i_Power do //列互换
75 begin
76 dou_temp := dou_Matrixs[i][0];
77 dou_Matrixs[i][0] := dou_Matrixs[i][j];
78 dou_Matrixs[i][j] := dou_temp;
79 end;
80 end
81 else //第一行全为0，则行列式的值=0
82 begin
83 MaxtrixValue := 0;
84 Exit;
85 end;
86 end;
87 //------------开始LU分解--------------------------------------------------//
88 BB := 1;
89 Rows := i_Power + 1; //行列式的阶数
90 //行列式分布如下
91 //a[0][0] a[0][1] a[0][2] ... a[0][Rows-1]
92 //a[1][0] a[1][1] a[1][2] ... a[1][Rows-1]
93 //.
94 //.
95 //.
96 //a[Rows-1][0] a[Rows-1][1] a[Rows-1][2] ... a[Rows-1][Rows-1]
97 for i := 0 to Rows - 1 do
98 begin
99 L[i][i] := 1; //获得L[i][i] = 1
100 U[0][i] := dou_Matrixs[0][i]; //获得U[0][i] = a[0][i]
101 end;
102 if Rows >= 2 then //阶数>=2
103 begin
104 for i := 1 to Rows - 1 do
105 L[i][0] := dou_Matrixs[i][0] / U[0][0]; //获得L[i][0] = a[i][0] / u[0][0]
106 end;
107 if (Rows >= 2) then
108 begin
109 CC := 1; //从第二行开始求 U[][]
110
111 SS2:
112 for i := CC to Rows - 1 do
113 begin
114 CC := i + 1; //下次循环从i+1行开始
115 for j := i to Rows - 1 do
116 begin
117 TempB := 0;
118 for k := 0 to i - 1 do
119 TempB := TempB + L[i][k] * U[k][j]; //∑(k=(0,..i-1))(L[i][k] * U[k][j])
120 U[i][j] := dou_Matrixs[i][j] - TempB; //a[i][j] -∑(k=(0,..i-1))(L[i][k] * U[k][j])
121 if ((i = j) and (abs(U[i][j]) <= 0.000001)) then //分解时对角线上的U不能等于0
122 begin
123 MaxtrixValue := 0;
124 Exit;
125 end;
126 end; //该行的U[][]已经求完了
127 if (i < Rows - 1) then //求下行的U时要用到下一行的L
128 begin
129 BB := i + 1;
130 goto SS1; //所以先去计算下行的L[][]
131 end
132 else
133 goto SS3; //所有的L U都计算完了
134 end; //end for i
135 end; //if(Rows >= 2 ) then
136 if (Rows >= 3) then //只有>=3阶的行列式才会再计算L值 <=2时L值可以直接得到
137 begin
138
139 SS1:
140 for m := BB to Rows - 1 do //从要计算的行开始计算
141 begin
142 for N := 0 to m - 1 do
143 begin
144 Temp := 0;
145 for k := 0 to N - 1 do
146 Temp := Temp + L[m][k] * U[k][N]; //∑(k=(0,..j-1) L[i][k] * U[k][j])
147 if U[N][N] <> 0 then
148 L[m][N] := (dou_Matrixs[m][N] - Temp) / U[N][N]; //(a[i][j] -∑(k=(0,..j-1) L[i][k] * U[k][j]) / u[j][j]
149 end; // end for N
150 if m <= Rows - 1 then //该行的L值已经求完
151 goto SS2; //去计算下行的U
152 end; // end for m
153 end; // end if
154
155 SS3:
156 dou_Result := U[0][0];
157 if Rows >= 2 then
158 begin
159 //行列式的值|A |= |L|*|U|= U[1][1]*U[2][2]*...U[m][m] (m=行列式的阶数) */
160 for i := 1 to Rows - 1 do
161 dou_Result := dou_Result * U[i][i];
162 end;
163 dou_Result := dou_Result * i_Plusminus; //得到行列式的值
164 MaxtrixValue := dou_Result;
165 //------------进行LU分解完成-------------------------------------------------------//
166 if not IsOnlyLu then
167 begin
168 dou_Y[0] := SumyCoordinate[0]; //得到y[0]
169 for i := 0 to i_Power - 1 do
170 begin
171 Temp := 0;
172 for k := 0 to i do
173 Temp := Temp + dou_Y[k] * L[i + 1][k]; // ∑(k=(0,..i) L[i][k] * y[k]
174 TempB := SumyCoordinate[i + 1]; //得到b[i]
175 dou_Y[i + 1] := TempB - Temp; //得到y[i]
176 end;
177 //求出方程式的(解)系数
178 if U[i_Power][i_Power] <> 0 then
179 Temp := dou_Y[i_Power] / U[i_Power][i_Power]; //得到x[n]
180
181 EquationParam[i_Power] := Temp; //将x[n]存入pEquationPara数组中
182
183 for i := i_Power - 1 downto 0 do
184 begin
185 Temp := 0;
186 for k := i + 1 to i_Power do
187 Temp := Temp + U[i][k] * EquationParam[k]; //∑(k=(n,..i+1) (U[i][k] * x[k])
188 if U[i][i] <> 0 then
189 TempB := (dou_Y[i] - Temp) / U[i][i]; //(y[i] - ∑(k=(n,..i+1) (U[i][k] * x[k]) ) / U[i][i]
190 EquationParam[i] := TempB; //将x[i]存入pEquationPara数组中
191 end;
192 end;
193 Result := True;
194 end;
1 //------------------------------------------------------------------------------
2 //功能：一个执行计算的过程 ，在执行前需要先设置最高次幂，是否排序，和坐标数组3项，
3 // 然后才能进行计算
4 //参数：
5 //
6 //返回值：
7 //
8
9 function TCMatrix.BeginComputer(const X, Y: array of Double;
10 const MaxPower: Integer = 1; const GetBestResult: Boolean = False;
11 const Sort: Boolean = False): Boolean;
12 var
13 dt: Extended;
14 i, j, MP, BestIndex: Integer;
15 EquationParaAry: array[0..MaxMi] of array of Double;
16 EquationParaResult: array[0..MaxMi] of Boolean;
17 R: array[0..MaxMi] of Double;
18 BestR: Double;
19 begin
20 Result := False;
21
22 if not SetXYCoordinate(X, Y) then
23 Exit;
24 FnMaxPower := MaxPower;
25 if FnMaxPower > MaxMi then
26 FnMaxPower := MaxMi;
27
28 FbIsNeedSort := Sort;
29 if (FbIsNeedSort) then //判断数组是否需要排序
30 SortTwoArray(FxCoordinate, FyCoordinate);
31
32 SetLength(FEquationPara, 0);
33 SetLength(FEquationPara, FnMaxPower + 1);
34 for i := 0 to High(FEquationPara) do
35 FEquationPara[i] := 0; //初始化系数数组
36 if FnCoordinateSize = 1 then //如果只有一组数据
37 begin
38 FEquationPara[0] := Y[0] - X[0];
39 FEquationPara[1] := 1;
40 Exit;
41 end;
42 if GetBestResult then
43 MP := MaxMi
44 else
45 MP := FnMaxPower;
46 FillChar(EquationParaResult, SizeOf(EquationParaResult), 0);
47 FillChar(R, SizeOf(R), 0);
48 try
49 for i := MP downto 0 do
50 begin
51 SetLength(FEquationPara, i + 1);
52 SetLength(EquationParaAry[i], i + 1);
53 R[i] := 0;
54 for j := 0 to High(FEquationPara) do
55 FEquationPara[j] := 0; //初始化系数数组
56 //通过输入的坐标求出正规方程组
57 Muti(i);
58 EquationParaResult[i] := MatrixLU(FSumxCoordinate, FSumyCoordinate, FEquationPara, i, False, dt);
59 if (not GetBestResult) and EquationParaResult[i] then
60 begin
61 Result := True;
62 Exit;
63 end;
64 if EquationParaResult[i] then
65 begin
66 for j := 0 to High(FEquationPara) do
67 EquationParaAry[i][j] := FEquationPara[j];
68 R[i] := CalR(FxCoordinate, FyCoordinate, EquationParaAry[i]);
69 end;
70 end;
71 if GetBestResult then
72 begin
73 BestIndex := -1;
74 for i := 0 to MP do
75 begin
76 if not EquationParaResult[i] then
77 Continue;
78 if BestIndex = -1 then
79 begin
80 BestIndex := i;
81 BestR := R[i];
82 end;
83 if BestIndex <> -1then // new
84 begin
85 if R[i] > BestR then
86 begin
87 BestIndex := i;
88 BestR := R[i];
89 end;
90 end;
91 end;
92 Result := BestIndex <> -1;
93 if Result then
94 for j := 0 to High(EquationParaAry[BestIndex]) do
95 begin
96 SetLength(FEquationPara, Length(EquationParaAry[BestIndex]));
97 FEquationPara[j] := EquationParaAry[BestIndex][j];
98 end;
99 end;
100 finally
101 for i := Low(EquationParaAry) to High(EquationParaAry) do
102 begin
103 SetLength(EquationParaAry[i], 0); EquationParaAry[i] := nil;
104 end;
105 end;
106 //pCoefficientx 是正规方程组左边的值 m_i_MaxPower 是拟合方程的最高次数
107 //False 表示求拟合方程 ，此处必须是false
108 end;
109
110 //------------------------------------------------------------------------------
111 //功能：计算行列式的值后，可以用这个函数得到行列式的值
112 //
113 //参数：
114 //
115 //返回值：返回行列式的值
116 //
117
118 function TCMatrix.GetMatrixValue: Double;
119 begin
120 Result := Fdou_MaxtrixValue;
121 end;
122
123 destructor TCMatrix.Destory;
124 begin
125 SetLength(FxCoordinate, 0);
126 FxCoordinate := nil;
127 SetLength(FyCoordinate, 0);
128 FyCoordinate := nil;
129 SetLength(FSumxCoordinate, 0);
130 FSumxCoordinate := nil;
131 SetLength(FSumyCoordinate, 0);
132 FSumyCoordinate := nil;
133 SetLength(FEquationPara, 0);
134 FEquationPara := nil;
135 end;
136
137 class function TCMatrix.CalR(const X, Y, EquationParam: array of Double): Double;
138 var
139 R1, R2, R3, Y1: array of Double;
140 XLen, YLen, i: Integer;
141 YAvg, Y1Avg: Double;
142 a, b, c, d: Extended;
143 begin
144 Result := 0;
145 XLen := Length(X);
146 YLen := Length(Y);
147 if (XLen <> YLen) or (XLen = 0) or (Low(X) <> Low(Y)) or (Length(EquationParam) = 0) then
148 Exit;
149 try
150 SetLength(R1, XLen);
151 SetLength(R2, XLen);
152 SetLength(R3, XLen);
153 SetLength(Y1, XLen);
154
155 CalY1Array(X, EquationParam, Y1);
156
157 YAvg := Mean(Y);
158 Y1Avg := Mean(Y1);
159 for i := Low(X) to High(X) do
160 begin
161 R1[i] := Y[i] * Y1[Low(Y1) + i - Low(X)];
162 R2[i] := Sqr(Y[i]);
163 R3[i] := Sqr(Y1[Low(Y1) + i - Low(X)]);
164 end;
165 a := Sum(R1) - XLen * YAvg * Y1Avg;
166 b := Sum(R2);
167 c := Sum(R3);
168 d := Sqrt(Abs(b - XLen * Sqr(YAvg))) * Sqrt(Abs(c - XLen * Sqr(Y1Avg))); //new
169 if d <> 0 then
170 Result := a / d;
171 if Result > 1 then
172 Result := 1;
173 if Result < -1 then
174 Result := -1;
175 finally
176 SetLength(R1, 0); R1 := nil;
177 SetLength(R2, 0); R2 := nil;
178 SetLength(R3, 0); R3 := nil;
179 SetLength(Y1, 0); Y1 := nil;
180 end;
181 end;
182
183 class function TCMatrix.CalY1(X: Double;
184 const EquationParam: array of Double): Double;
185 var
186 i: Integer;
187 begin
188 Result := 0;
189 for i := High(EquationParam) downto 0 do
190 begin
191 Result := Result + EquationParam[i] * Power(X, i);
192 end;
193 end;
194
195 class function TCMatrix.CalY1Array(const X,
196 EquationParam: array of Double;var Y1: array of Double): Boolean;
197 var
198 i: Integer;
199 begin
200 Result := False;
201 if (Length(X) <> Length(Y1))then
202 Exit;
203 for i := Low(X) to High(X) do
204 begin
205 Y1[Low(Y1) + i - Low(X)] := CalY1(X[i], EquationParam);
206 end;
207 Result := True;
208 end;
209
210 class function TCMatrix.CalRMSE(const X, Y,
211 EquationParam: array of Double): Double;
212 var
213 XLen, YLen, i: Integer;
214 Y1: array of Double;
215 Mean, StdMean: Extended;
216 begin
217 XLen := Length(X);
218 YLen := Length(Y);
219 if (XLen <> YLen) or (XLen = 0) or (Low(X) <> Low(Y)) or (Length(EquationParam) = 0) then
220 Exit;
221 SetLength(Y1, XLen);
222 try
223 CalY1Array(X, EquationParam, Y1);
224 for i := Low(X) to High(X) do
225 begin
226 Y1[Low(Y1) + i - Low(X)] := Y[i] - Y1[Low(Y1) + i - Low(X)];
227 end;
228 MeanAndStdDev(Y1, Mean, StdMean);
229 Result := StdMean;
230 finally
231 SetLength(Y1, 0); Y1 := nil;
232 end;
233 end;
234
235 end.
上面是功能代码...分了好多段..
demo unit pas
1 unit Unit1;
2
3 interface
4
5 uses
6 Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
7 StdCtrls, CYHSMatrixUnit, Math, Spin, ExtCtrls, TeeProcs, TeEngine, Chart,
8 Series;
9
10 type
11 TForm1 = class(TForm)
12 Button1: TButton;
13 Memo1: TMemo;
14 Memo2: TMemo;
15 Memo3: TMemo;
16 Label1: TLabel;
17 Label2: TLabel;
18 Label3: TLabel;
19 Button2: TButton;
20 Button3: TButton;
21 Memo4: TMemo;
22 Label4: TLabel;
23 SpinEdit1: TSpinEdit;
24 Chart1: TChart;
25 Series1: TFastLineSeries;
26 Series2: TFastLineSeries;
27 CheckBox1: TCheckBox;
28 procedure Button1Click(Sender: TObject);
29 procedure Button2Click(Sender: TObject);
30 procedure Button3Click(Sender: TObject);
31 private
32 { Private declarations }
33 public
34 { Public declarations }
35 end;
36
37 var
38 Form1: TForm1;
39
40 implementation
41
42 {$R *.DFM}
43
44 procedure TForm1.Button1Click(Sender: TObject);
45 var
46 Obj: TCMatrix;
47 X, Y: array of Double;
48 i, j, Code: Integer;
49 Value, p: Extended;
50 str: string;
51 begin
52 Obj := TCMatrix.Create;
53 try
54 SetLength(X, Memo1.Lines.Count);
55 SetLength(Y, Memo2.Lines.Count);
56 if High(X) <> High(Y) then
57 Exit;
58 for i := 0 to Memo1.Lines.Count - 1 do
59 begin
60 Val(Trim(Memo1.Lines[i]), Value, Code);
61 if Code <> 0 then
62 Exit;
63 X[i] := Value;
64 Val(Trim(Memo2.Lines[i]), Value, Code);
65 if Code <> 0 then
66 Exit;
67 Y[i] := Value;
68 end;
69 Chart1.SeriesList.ClearValues;
70 Memo3.Clear;
71 if not Obj.BeginComputer(X, Y, SpinEdit1.Value, CheckBox1.Checked, False) then
72 Exit;
73 Label4.Caption := 'R:'+ FloatToStr(Obj.CalR(X, Y, Obj.EquationPara)) + #13#10 +
74 'RMSE:' + FloatToStr(Obj.CalRMSE(X, Y, Obj.EquationPara));
75 for i := High(Obj.EquationPara) downto 0 do
76 begin
77 Memo3.Lines.Add(Format('x^%d: ',[i]) + FloatToStr(Obj.EquationPara[i]));
78 end;
79 Memo4.Clear;
80 Memo4.Lines.BeginUpdate;
81 for i := 0 to High(X) do
82 begin
83 //Value := Obj.CalY1(X[i], Obj.EquationPara);
84 Value := 0;
85 for j := High(Obj.EquationPara) downto 0 do
86 begin
87 Value := Value + Obj.EquationPara[j] * Power(X[i], j);
88 end;
89 p := 0;
90 if Value <> 0 then
91 p := Abs((Y[i] - Value) / Value) * 100;
92 str := Format('%-30s %.4f%%', [FloatToStr(Value), p]);
93 Memo4.Lines.Add(str);
94 Series1.AddXY(X[i], Y[i]);
95 Series2.AddXY(X[i], Value);
96 end;
97 Memo4.Lines.EndUpdate;
98 finally
99 Obj.Free;
100 SetLength(X, 0);
101 X := nil;
102 SetLength(Y, 0);
103 Y := nil;
104 end;
105
106 end;
107
108 procedure TForm1.Button2Click(Sender: TObject);
109 begin
110 Memo1.Clear;
111 end;
112
113 procedure TForm1.Button3Click(Sender: TObject);
114 begin
115 Memo2.Clear;
116 end;
117
118 end.
demo unit dfm 看来也要分段了.
1 object Form1: TForm1
2 Left = 173
3 Top = 35
4 Width = 968
5 Height = 724
6 Caption = 'Form1'
7 Color = clBtnFace
8 Font.Charset = DEFAULT_CHARSET
9 Font.Color = clWindowText
10 Font.Height = -11
11 Font.Name = 'MS Sans Serif'
12 Font.Style = []
13 OldCreateOrder = False
14 PixelsPerInch = 96
15 TextHeight = 13
16 object Label1: TLabel
17 Left = 112
18 Top = 8
19 Width = 5
20 Height = 13
21 Caption = 'x'
22 end
23 object Label2: TLabel
24 Left = 264
25 Top = 8
26 Width = 5
27 Height = 13
28 Caption = 'y'
29 end
30 object Label3: TLabel
31 Left = 416
32 Top = 16
33 Width = 30
34 Height = 13
35 Caption = 'Param'
36 end
37 object Label4: TLabel
38 Left = 592
39 Top = 14
40 Width = 30
41 Height = 13
42 Caption = 'Result'
43 end
44 object Button1: TButton
45 Left = 856
46 Top = 8
47 Width = 75
48 Height = 25
49 Caption = 'Calculate'
50 TabOrder = 3
51 OnClick = Button1Click
52 end
53 object Memo1: TMemo
54 Left = 104
55 Top = 40
56 Width = 129
57 Height = 297
58 Lines.Strings = (
59 '-63.15'
60 '-62.15'
61 '-61.15'
62 '-60.15'
63 '-59.15'
64 '-58.15'
65 '-57.15'
66 '-56.15'
67 '-55.15'
68 '-54.15'
69 '-53.15'
70 '-52.15'
71 '-51.15'
72 '-50.15'
73 '-49.15'
74 '-48.15'
75 '-47.15'
76 '-46.15'
77 '-45.15'
78 '-44.15'
79 '-43.15'
80 '-42.15'
81 '-41.15'
82 '-40.15'
83 '-39.15'
84 '-38.15'
85 '-37.15'
86 '-36.15'
87 '-35.15'
88 '-34.15'
89 '-33.15'
90 '-32.15'
91 '-31.15'
92 '-30.15'
93 '-29.15'
94 '-28.15'
95 '-27.15'
96 '-26.15'
97 '-25.15'
98 '-24.15'
99 '-23.15'
100 '-22.15'
101 '-21.15'
102 '-20.15'
103 '-19.15'
104 '-18.15'
105 '-17.15'
106 '-16.15'
107 '-15.15'
108 '-14.15'
109 '-13.15'
110 '-12.15'
111 '-11.15'
112 '-10.15'
113 '-9.15'
114 '-8.15'
115 '-7.15'
116 '-6.15'
117 '-5.15'
118 '-4.15'
119 '-3.15'
120 '-2.15'
121 '-1.15'
122 '-0.15'
123 '0.85'
124 '1.85'
125 '2.85'
126 '3.85'
127 '4.85'
128 '5.85'
129 '6.85'
130 '7.85'
131 '8.85'
132 '9.85'
133 '10.85'
134 '11.85'
135 '12.85'
136 '13.85'
137 '14.85'
138 '15.85'
139 '16.85'
140 '17.85'
141 '18.85'
142 '19.85'
143 '20.85'
144 '21.85'
145 '22.85'
146 '23.85'
147 '24.85'
148 '25.85'
149 '26.85'
150 '27.85'
151 '28.85'
152 '29.85'
153 '30.85'
154 '31.85'
155 '32.85'
156 '33.85'
157 '34.85'
158 '35.85'
159 '36.85'
160 '37.85'
161 '38.85'
162 '39.85'
163 '40.85'
164 '41.85'
165 '42.85'
166 '43.85'
167 '44.85'
168 '45.85'
169 '46.85'
170 '47.85'
171 '48.85'
172 '49.85'
173 '50.85'
174 '51.85'
175 '52.85'
176 '53.85'
177 '54.85'
178 '55.85'
179 '56.85'
180 '57.85'
181 '58.85'
182 '59.85'
183 '60.85'
184 '61.85'
185 '62.85'
186 '63.85'
187 '64.85'
188 '65.85'
189 '66.85'
190 '67.85'
191 '68.85'
192 '69.85'
193 '70.85'
194 '71.85'
195 '72.85'
196 '73.85'
197 '74.85'
198 '75.85'
199 '76.85'
200 '77.85'
201 '78.85'
202 '79.85'
203 '80.85'
204 '81.85'
205 '82.85'
206 '83.85'
207 '84.85'
208 '85.85'
209 '86.85'
210 '87.85'
211 '88.85'
212 '89.85'
213 '90.85'
214 '91.85'
215 '92.85'
216 '93.85'
217 '94.85'
218 '95.85'
219 '96.85'
220 '97.85'
221 '98.85'
222 '99.85'
223 '100.85')
1 ScrollBars = ssVertical
2 TabOrder = 5
3 end
4 object Memo2: TMemo
5 Left = 256
6 Top = 40
7 Width = 137
8 Height = 297
9 Lines.Strings = (
10 '0.00054 '
11 '0.00059 '
12 '0.00064 '
13 '0.00070 '
14 '0.00076 '
15 '0.00082 '
16 '0.00089 '
17 '0.00097 '
18 '0.00104 '
19 '0.00113 '
20 '0.00122 '
21 '0.00132 '
22 '0.00142 '
23 '0.00153 '
24 '0.00164 '
25 '0.00177 '
26 '0.00190 '
27 '0.00204 '
28 '0.00219 '
29 '0.00235 '
30 '0.00251 '
31 '0.00269 '
32 '0.00288 '
33 '0.00308 '
34 '0.00329 '
35 '0.00351 '
36 '0.00375 '
37 '0.00399 '
38 '0.00426 '
39 '0.00453 '
40 '0.00482 '
41 '0.00513 '
42 '0.00545 '
43 '0.00579 '
44 '0.00615 '
45 '0.00653 '
46 '0.00692 '
47 '0.00734 '
48 '0.00777 '
49 '0.00823 '
50 '0.00871 '
51 '0.00921 '
52 '0.00973 '
53 '0.01028 '
54 '0.01086 '
55 '0.01146 '
56 '0.01208 '
57 '0.01274 '
58 '0.01343 '
59 '0.01414 '
60 '0.01489 '
61 '0.01567 '
62 '0.01648 '
63 '0.01733 '
64 '0.01821 '
65 '0.01912 '
66 '0.02008 '
67 '0.02107 '
68 '0.02211 '
69 '0.02318 '
70 '0.02430 '
71 '0.02546 '
72 '0.02666 '
73 '0.02791 '
74 '0.02921 '
75 '0.03055 '
76 '0.03195 '
77 '0.03340 '
78 '0.03489 '
79 '0.03645 '
80 '0.03805 '
81 '0.03972 '
82 '0.04144 '
83 '0.04322 '
84 '0.04506 '
85 '0.04697 '
86 '0.04893 '
87 '0.05097 '
88 '0.05307 '
89 '0.05524 '
90 '0.05748 '
91 '0.05979 '
92 '0.06217 '
93 '0.06463 '
94 '0.06716 '
95 '0.06978 '
96 '0.07247 '
97 '0.07524 '
98 '0.07810 '
99 '0.08104 '
100 '0.08407 '
101 '0.08719 '
102 '0.09040 '
103 '0.09370 '
104 '0.09709 '
105 '0.10058 '
106 '0.10417 '
107 '0.10785 '
108 '0.11164 '
109 '0.11553 '
110 '0.11952 '
111 '0.12363 '
112 '0.12784 '
113 '0.13216 '
114 '0.13660 '
115 '0.14115 '
116 '0.14581 '
117 '0.15060 '
118 '0.15550 '
119 '0.16053 '
120 '0.16568 '
121 '0.17096 '
122 '0.17637 '
123 '0.18191 '
124 '0.18759 '
125 '0.19339 '
126 '0.19934 '
127 '0.20542 '
128 '0.21165 '
129 '0.21802 '
130 '0.22454 '
131 '0.23120 '
132 '0.23802 '
133 '0.24498 '
134 '0.25211 '
135 '0.25938 '
136 '0.26682 '
137 '0.27442 '
138 '0.28218 '
139 '0.29011 '
140 '0.29821 '
141 '0.30647 '
142 '0.31491 '
143 '0.32353 '
144 '0.33232 '
145 '0.34129 '
146 '0.35045 '
147 '0.35978 '
148 '0.36931 '
149 '0.37902 '
150 '0.38893 '
151 '0.39902 '
152 '0.40932 '
153 '0.41981 '
154 '0.43051 '
155 '0.44141 '
156 '0.45251 '
157 '0.46382 '
158 '0.47534 '
159 '0.48708 '
160 '0.49903 '
161 '0.51120 '
162 '0.52359 '
163 '0.53620 '
164 '0.54904 '
165 '0.56211 '
166 '0.57541 '
167 '0.58894 '
168 '0.60271 '
169 '0.61671 '
170 '0.63096 '
171 '0.64545 '
172 '0.66018 '
173 '0.67516 '
174 '0.69040 ')
175 ScrollBars = ssVertical
176 TabOrder = 6
177 end
178 object Memo3: TMemo
179 Left = 408
180 Top = 40
181 Width = 161
182 Height = 297
183 Lines.Strings = (
184 'Memo3')
185 TabOrder = 7
186 end
187 object Button2: TButton
188 Left = 136
189 Top = 8
190 Width = 75
191 Height = 25
192 Caption = 'clear memox'
193 TabOrder = 0
194 OnClick = Button2Click
195 end
196 object Button3: TButton
197 Left = 288
198 Top = 8
199 Width = 75
200 Height = 25
201 Caption = 'clear memoy'
202 TabOrder = 1
203 OnClick = Button3Click
204 end
205 object Memo4: TMemo
206 Left = 600
207 Top = 40
208 Width = 345
209 Height = 297
210 Lines.Strings = (
211 'Memo4')
212 ScrollBars = ssVertical
213 TabOrder = 8
214 end
215 object SpinEdit1: TSpinEdit
216 Left = 464
217 Top = 8
218 Width = 65
219 Height = 22
220 MaxValue = 20
221 MinValue = 0
222 TabOrder = 2
223 Value = 1
224 end
225 object Chart1: TChart
226 Left = 96
227 Top = 352
228 Width = 833
229 Height = 329
230 Legend.Visible = False
231 Title.Text.Strings = (
232 'TChart')
233 Title.Visible = False
234 View3D = False
235 TabOrder = 9
236 object Series1: TFastLineSeries
237 Marks.Arrow.Visible = True
238 Marks.Callout.Brush.Color = clBlack
239 Marks.Callout.Arrow.Visible = True
240 Marks.Visible = False
241 LinePen.Color = clRed
242 XValues.Name = 'X'
243 XValues.Order = loAscending
244 YValues.Name = 'Y'
245 YValues.Order = loNone
246 end
247 object Series2: TFastLineSeries
248 Marks.Arrow.Visible = True
249 Marks.Callout.Brush.Color = clBlack
250 Marks.Callout.Arrow.Visible = True
251 Marks.Visible = False
252 SeriesColor = clBlue
253 LinePen.Color = clBlue
254 XValues.Name = 'X'
255 XValues.Order = loAscending
256 YValues.Name = 'Y'
257 YValues.Order = loNone
258 end
259 end
260 object CheckBox1: TCheckBox
261 Left = 800
262 Top = 12
263 Width = 49
264 Height = 17
265 Caption = 'Best'
266 TabOrder = 4
267 end
268 end