关于最小二乘法，帮忙看看这段代码啥意思？ 啊-----加分，100分吧 (10分)

L

lcl_003

Unregistered / Unconfirmed
GUEST, unregistred user!
这段代码是关于最小二乘法的，有几 个地方不太明白主要是调用的时候

procedure MinSqrMul(m, n: Integer; A: TMatrix; b: array of Double; var x: array of Double);
var
i, j, k: Integer;
ATA: array of array of Double;
ATb: array of Double;
Elem: Double;
begin
SetLength(ATA, n, n);
SetLength(ATb, n);
for i := 0 to n - 1 do
for j := 0 to n - 1 do begin
ATA[i, j] := 0;
for k := 0 to m - 1 do
ATA[i, j] := ATA[i, j] + A[k, i] * A[k, j];
end;
for i := 0 to n - 1 do begin
ATb := 0;
for j := 0 to m - 1 do
ATb := ATb + A[j, i] * b[j];
end;
for i := 0 to n - 1 do begin
Elem := ATA[i, i];
for j := i to n - 1 do
ATA[i, j] := ATA[i, j] / Elem;
ATb := ATb / Elem;
for k := i + 1 to n - 1 do begin
Elem := -ATA[k, i];
for j := [red]i + 1[/red] to n - 1 do
ATA[k, j] := ATA[k, j] + ATA[i, j] * Elem;
ATb[k] := ATb[k] + ATb * Elem;
end;
end;
x[n - 1] := ATb[n - 1];
for i := n - 2 downto 0 do begin
for j := i + 1 to n - 1 do
ATb := ATb - ATA[i, j] * x[j];
x := ATb;
end;
end;


以下是调用
procedure TForm1.FormClick(Sender: TObject);
var
A: TMatrix;
b: array [0..2] of Double;
x: array [0..1] of Double;
begin
SetLength(A, 3, 2);
A[0, 0] := 1; A[0, 1] := 2; b[0] := 3;
A[1, 0] := 1; A[1, 1] := 2; b[1] := 4;
A[2, 0] := 2; A[2, 1] := 1; b[2] := 3;
MinSqrMul(3, 2, A, b, x);
ShowMessage(Format('%f, %f', [x[0], x[1]]));
end;

我现在不明白调用的时候各个数据是什么。MinSqrMul(3, 2, A, b, x); 中
A应该是各点吧，那3，2，b,x是什么呢？
 
3，2是维数？（没有细看代码）
x是返还结果
 
这是一个最小二乘法，没错。
按照上面的调用（3，2），应该是一个线性回归.
一般的线性公式如下：Y=aX+b
主要是根据 各个点(Xi,Yi)求出 a 和 b
不难看出这里应该是 3 个点，那么如何调用呢？

假设 我们的三个点为 [X1,Y1],[X2,Y2],[X3,Y3]
MinSqrMul(3, 2, A, b, x)中的：
A 是一个三行二列的数组（3X2矩阵） [[ X1, 1.0 ],[ x2 , 1.0 ],[ X3 , 1.0 ]]
b 为 [ Y1 , Y2 , Y3 ]
x 是返回值，为 [a,b] 为里的 a,b 是公式中 “Y=aX+b” 的 a,b
a 是斜率， b 是载距

但这里的算法比较粗糙，在计算逆矩阵时，没有考虑 A'*A 对角线元素是否为0，
主要表现在进行以下除法时没有判断 Elem 的值：
ATA[i, j] := ATA[i, j] / Elem;
ATb := ATb / Elem;
应该先判断，然后进行必要的行交换。
以下的调用赋值可能有问题，没法解释：
A[0, 0] := 1; A[0, 1] := 2; b[0] := 3;
A[1, 0] := 1; A[1, 1] := 2; b[1] := 4;
A[2, 0] := 2; A[2, 1] := 1; b[2] := 3;
至少要将 A 的某一列全赋为 1 ，
如果第一列全赋为 1 ，那么返回的 x 为 [b,a]
如果第二列全赋为 1 ，那么返回的 x 为 [a,b]

以下赋值应该得到结果 x=[1,0]
A[0, 0] := 1; A[0, 1] := 1; b[0] := 1;
A[1, 0] := 2; A[1, 1] := 1; b[1] := 2;
A[2, 0] := 3; A[2, 1] := 1; b[2] := 3;
 
jsxjd, 多谢你的解答，不过还是有点不明白：
》A 是一个三行二列的数组（3X2矩阵） [[ X1, 1.0 ],[ x2 , 1.0 ],[ X3 , 1.0 ]]
》b 为 [ Y1 , Y2 , Y3 ]

也就是说A是个数组，代表三个点的横坐标？b代表纵坐标？那A为什么用数组表示，我表示
横坐标用一个数就可以了啊，干吗用数组？

》x 是返回值，为 [a,b] 为里的 a,b 是公式中 “Y=aX+b” 的 a,b
》 a 是斜率， b 是载距

x是返回值，返回什么？好象它没用啊。
 
上面的代码也适用于多元回归,不一定是线性的：
比如 Z = aX +bY +c
更一般地，
Y = An*Xn + An-1*Xn-1 +......+ A1*X1 + A0
 
也就是说A是个数组，代表三个点的横坐标？b代表纵坐标？那A为什么用数组表示，我表示
横坐标用一个数就可以了啊，干吗用数组？

x是返回值，返回什么？好象它没用啊。

我的问题还没回答啊

 
哦，我现在要割据20个点（不规则）拟合出一条直线，如下：
。
。
。
。
。

在这二十个点中间根据最小二乘法画一条直线，然后求直线aX+b=y的a,b值。

我刚试了试上面的算法，结果是对的。但是只是做了两个点的。我是这么写的：
procedure TForm1.Button1Click(Sender: TObject);
var
A: TMatrix;
b: array [0..1] of Double;
x: array [0..1] of Double;
begin
SetLength(A, 2, 2);
A[0, 0] := 1; A[0, 1] := 2; b[0] := 9;
A[1, 0] := 2; A[1, 1] := 3 b[1] := 17;

MinSqrMul(2, 2, A, b, x);
ShowMessage(Format('%f, %f', [x[0], x[1]]));
end;

这样得出的方程是x+2y=9,2x+3y=17吗？可我要求的是a,b啊，如果是20个点该怎么做？
 
没人会啊，jsxjd, 在不在
 
你的写法好象不对，我马上举个例子给你。
 
好的，多谢。
 
首先要知道在以下声明中：
MinSqrMul(m, n: Integer; A: TMatrix; b: array of Double; var x: array of Double);

m是点的数量。
n是所求参数的个数，
A 是自变量的值（矩阵 n X m ）
b 是应变量的值（n 维的列向量）

假设现在有 4 个平面上的点(x,y)：
[0.0,1.0],[1.0,3.0],[2.0,5.0],[3.0,7.0]
很容易自己画一下，就可以得到 Y=2X+1,因为四个点都在同一
直线上，回归的结果也应该是这样（不考虑计算误差）。

我不知道你的TMatrix是怎么声明的，它很可能是二维动态浮点数数组。
在这里应是4行，2列，且如下(是各自变量的值，应将常数项 b 理解为
1*b，这儿也有一个自变量，只是不变，总是为1而已)：
[ [0.0,1.0],
[1.0,1.0],
[2.0,1.0],
[3.0,1.0] ]
“b” 就是各点的 Y值。 [1.0,3.0,5.0,7.0]

m 应为 4(四个点), 因为我们在 aX+b=Y 中要求两个参数，所有 n=2。

上面所举例子的计算结果应为 x=[2.0,1.0]，你可以试一下。

该过程的不足在前面已经讲过，如要真正使用该过程，应进行交换，
还要考虑逆矩阵不存在的情况，这时有无数个解。
 
jsxjd, 非常感谢。还要请教：）

形如ax+b=y，a,b是要求的数，那你下面的话：
“
应将常数项 b 理解为
1*b，这儿也有一个自变量，只是不变，总是为1而已)：
[ [0.0,1.0],
[1.0,1.0],
[2.0,1.0],
[3.0,1.0] ]
“b” 就是各点的 Y值。 [1.0,3.0,5.0,7.0]
”

你上面说“应将常数项 b 理解为
1*b，这儿也有一个自变量，只是不变，总是为1而已)”但是最后有说
““b” 就是各点的 Y值。 [1.0,3.0,5.0,7.0]”有点矛盾？糊涂了。

看了你的解释你看我的理解对不。
按我的情况：20点拟合一条直线。那m是20，n是2，A是20个点x,y坐标的矩阵，对吗？
 
以下是调用(20个点，应说样本数为20)
procedure TForm1.FormClick(Sender: TObject);
var
A: TMatrix;
b: array of Double;
x: array of Double;
begin
SetLength(A, 20, 2);
A[0, 0] := 1; A[0, 1] := 2; b[0] := 3;
A[1, 0] := 1; A[1, 1] := 2; b[1] := 4;
A[2, 0] := 2; A[2, 1] := 1; b[2] := 3;
........
A[19, 0] := 2; A[19, 1] := 1; b[19] := 3;

SetLength(b, 20);
SetLength(x, 2);
MinSqrMul(20, 2, A, b, x);
ShowMessage(Format('%f, %f', [x[0], x[1]]));
end;
 
应将常数项 b 理解为
1*b，这儿也有一个自变量，只是不变，总是为1而已)：
>>>>>>>>>>>>>>>> 这里的b,是 y=aX+b的b ,是要 求的，最后在 x ([a,b])中返回。
==================================
““b” 就是各点的 Y值。 [1.0,3.0,5.0,7.0]”
>>>>>>>>>>>>>>>>>>>> 这里的 b 是你以下过程声明中的 b
MinSqrMul(m, n: Integer; A: TMatrix; b: array of Double; var x: array of Double);
^<<<<<<<<<<<<<<<<<<<<<<<<<
===========================================
按我的情况：20点拟合一条直线。那m是20，n是2，A是20个点x,y坐标的矩阵，对吗？
!!!!!!! 不对 !!!!!!!!!
假设有 20 个点 [x1,y1],[x2,y2],...[x20,y20]
你的函数声明中的 A 如下：
[ [x1,1.0],
[x2,1.0],
[x3,1.0],
................
[x20,1.0] ]
============== A 中的第二列全为 1.0 ============================

你的函数声明中的 b 如下：
[ y1,y2,y3,...,y20]
m=20,n=2 是正确的。

 
jsxjd, 我按你的例子写了，结果正确，但是为什么

》
假设有 20 个点 [x1,y1],[x2,y2],...[x20,y20]
你的函数声明中的 A 如下：
[ [x1,1.0],
[x2,1.0],
[x3,1.0],
................
[x20,1.0] ]
============== A 中的第二列全为 1.0 ============================

》
为什么这里第二列全为1。0？

procedure TForm1.Button1Click(Sender: TObject);
var
A: TMatrix;
b: array [0..3] of double;
x: array [0..1] of double;
begin
SetLength(A, 4, 2);
A[0, 0] := 0; A[0, 1] := 1; b[0] := 1;
A[1, 0] := 1; A[1, 1] := 1; b[1] := 3;
A[2, 0] := 2; A[2, 1] := 1; b[2] := 5;
A[3, 0] := 3; A[3, 1] := 1; b[3] := 7;
MinSqrMul(4, 2, A, b, x);
ShowMessage(Format('%f, %f', [x[0], x[1]]));
end;

结果2，1
 
哦，知道了，看了你的第一个帖子

>>至少要将 A 的某一列全赋为 1 ，
如果第一列全赋为 1 ，那么返回的 x 为 [b,a]
如果第二列全赋为 1 ，那么返回的 x 为 [a,b]
>>
可是这是为什么?
 
jsxjd,非常感谢， 留个qq吧，我的：33847540
 
第二列放 1 ，公式是这样的 aX+b=y
第一列放 1, 公式变成 b+aX =y
 
哦，明白。明天给开贴
 
接受答案了.
 
您必须登录或注册才可回复。

Similar threads

L
JohnsonGuo在不在？关于最小二乘法(10分)
回复
13
查看
161
lcl_003
L
T
关于最小二乘法解线性方程组的源程序(30分)
回复
14
查看
95
JohnsonGuo
J
L
平方根法解方程组的程序，知错不知改。请教。 (50分)
回复
3
查看
133
lj.ah
L
H
帮忙看看，用字符检索汉字?(100分)
回复
2
查看
125
cxz9
C
I
数组问题，不知道哪里错了，不好意思只有5分了。(5分)
回复
4
查看
79
ISY
I
顶部