线性最小二乘法系数与常数的求解公式
例子
#线性最小二乘法
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