线性最小二乘法的Python程序

线性最小二乘法系数与常数的求解公式

例子

#线性最小二乘法
X = [0.698132, 0.959931, 1.134464, 1.570796, 1.919862]
Y = [0.188224, 0.209138, 0.230052, 0.250965, 0.313707]

n = len(X)
xy, x, y, x2,  = 0, 0, 0, 0

for i in range(n):
    xy += X[i]*Y[i]  #xy的和
    x += X[i]  #x的和
    y += Y[i]  #y的和
    x2 += X[i]**2  #x的平方的和

a1 = (n*xy - x*y) / (n*x2 - x**2)  #系数K
a0 = (x2*y - x*xy) / (n*x2 - x**2)  #常数b
print('a1 =',a1)
print('a0 =',a0)

结果：

a1 = 0.09609143373278023 
a0 = 0.11766514898834005

因此线性最小二乘拟合为
T = 0.09609143373278023 * θ + 0.11766514898834005

特殊点的线性最小二乘法，即截距为0

如果曲线过（0，0）点，线性最小二乘法的拟合曲线的截距应为0，那么其斜率的计算又是另一种公式了。




代码：

#线性最小二乘法
X = [0, 0.183, 0.36, 0.5324, 0.702, 0.867, 1.0244, 1.1774, 1.329, 1.479, 1.5, 1.56]
Y = [0, 306, 612, 917, 1223, 1529, 1835, 2140, 2446, 2752, 2767, 2896]

n = len(X)
xy, x, y, x2,  = 0, 0, 0, 0

if X.index(0) == Y.index(0):
    a0 = 0  #常数b为0,即截距为0
    for i in range(n):
        xy += X[i]*Y[i]  #xy的和
        x2 += X[i]**2  #x的平方的和
    a1 = xy / x2  #系数K
    print('a1 =',a1)
    print('a0 =',a0)
else:
    for i in range(n):
        xy += X[i]*Y[i]  #xy的和
        x += X[i]  #x的和
        y += Y[i]  #y的和
        x2 += X[i]**2  #x的平方的和

    a1 = (n*xy - x*y) / (n*x2 - x**2)  #系数K
    a0 = (x2*y - x*xy) / (n*x2 - x**2)  #常数b
    print('a1 =',a1)
    print('a0 =',a0)

结果：

a1 = 1828.3740666629756 
a0 = 0

上面的程序只能在有（0,0）的情况下正常运行，如果还存在（0，y）和（x，0）就不能正常运行。
所以进行了改进：

#线性最小二乘法
X = [0, 0.183, 0.36, 0.5324, 0.702, 0.867, 1.0244, 1.1774, 1.329, 1.479, 1.5, 1.56]
Y = [0, 306, 612, 917, 1223, 1529, 1835, 2140, 2446, 2752, 2767, 2896]

n = len(X)
xy, x, y, x2,  = 0, 0, 0, 0
flag = 1
#判断趋势线过不过(0,0)点
for i,j in zip(X,Y):
    if i ==0 and j == 0:
        flag = 0
#求解系数
if flag == 0:
    a0 = 0  #常数b为0,即截距为0
    for i in range(n):
        xy += X[i]*Y[i]  #xy的和
        x2 += X[i]**2  #x的平方的和
    a1 = xy / x2  #系数K
    print('a1 =',a1)
    print('a0 =',a0)
else:
    for i in range(n):
        xy += X[i]*Y[i]  #xy的和
        x += X[i]  #x的和
        y += Y[i]  #y的和
        x2 += X[i]**2  #x的平方的和

    a1 = (n*xy - x*y) / (n*x2 - x**2)  #系数K
    a0 = (x2*y - x*xy) / (n*x2 - x**2)  #常数b
    print('a1 =',a1)
    print('a0 =',a0)

结果：

a1 = 1828.3740666629756
a0 = 0