最⼩⼆乘法Python 实现
最⼩⼆乘法
数学和统计上⾯⼀个基本⽅法是，根据最⼩⼆衬发拟合平⾯上的点集。其拟合的图形通常是基本类型函数，如：线性函数、多项式、三⾓多项式等。由于数据有测量误差或者试验误差，我们不要求数据通过所有数据点。实际上，我们需要的是在所有数据点的y值，和逼近曲线相应点处的y值俩者之间误差的平⽅和最⼩意义下的最佳曲线。
具体原理推导不详细说，在线性代数的书⾥基本都会有介绍。下⾯介绍定理：
若A是秩为n的m×n的矩阵，则正规⽅程组：
有唯⼀解：
且x为⽅程组 AX = Y最⼩⼆乘解。
按书写习惯，X和Y都是列向量。参考 。由于输⼊数据格式不⼀样，刚开始看的时候没怎么看的懂。毕竟输⼊的是⾏向量，公式有点不同：。这⾥⽣成的数据x和y都是⾏向量，所以使⽤的是公式（2），代码如下：
A AX =T A Y
T x =(A A )A Y （1）
T −1T ω=(XX )XY （2）T −1T
"""最⼩⼆乘法"""
import numpy as np
import matplotlib.pyplot as plt
def fun2ploy(x,n):
'''
数据转化为[x^0,x^1,x^2,...x^n]
⾸列变1
'''
lens = len(x)
X = np.ones([1,lens])
for i in range(1,n):
X = np.vstack((X,np.power(x,i)))#按⾏堆叠
return X
def leastseq_byploy(x,y,ploy_dim):
'''
最⼩⼆乘求解
'''
#散点图
plt.scatter(x,y,color="r",marker='o',s = 50)
X = fun2ploy(x,ploy_dim);
#直接求解
Xt = X.transpose();#转置变成列向量
XXt=X.dot(Xt);#矩阵乘
XXtInv = np.linalg.inv(XXt)#求逆
XXtInvX = XXtInv.dot(X)
coef = XXtInvX.dot(y.T)
y_est = Xt.dot(coef)
return y_est,coef
def fit_fun(x):
'''
要拟合的函数
'''
#return np.power(x,5)
return np.sin(x)
#    return 5*x+3
#    return np.sqrt(25-pow(x-5,2))
if __name__ == '__main__':
data_num = 100;
ploy_dim =10;#拟合参数个数，即权重数量
noise_scale = 0.2;
## 数据准备
x = np.array(np.linspace(-2*np.pi,2*np.pi,data_num))  #数据
y = fit_fun(x)+noise_scale*np.random.rand(1,data_num)  #添加噪声
# 最⼩⼆乘拟合
[y_est,coef] = leastseq_byploy(x,y,ploy_dim)
#显⽰拟合结果
org_data = plt.scatter(x,y,color="r",marker='o',s = 50)
est_data = plt.plot(x,y_est,color="g",linewidth= 3)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Fit funtion with leastseq method")
plt.legend(["Noise data","Fited function"]);
plt.show()
选择性注释掉函数中的返回值，选择不同的曲线作为被拟合对象。
然后按照我⾃⼰的理解，修改了其中的算法，才搞清楚。代码如下（注意运算按公式（1），列向量计算）：
fit un ()f
"""最⼩⼆乘法"""
import numpy as np
import matplotlib.pyplot as plt
def fun2ploy(x,n):
'''
数据转化为[x^0,x^1,x^2,...x^n]
⾸列变1
'''
lens = len(x)
X = np.ones([1,lens])
for i in range(1,n):
X = np.vstack((X,np.power(x,i)))#按⾏堆叠
Xt = X.transpose()
return  Xt
def leastseq_byploy(x,y,ploy_dim):
'''
最⼩⼆乘求解
'''
#散点图
plt.scatter(x,y,color="r",marker='o',s = 50)
X = fun2ploy(x,ploy_dim);
#直接求解
Xt = X.transpose();#转置变成列向量
XtX=Xt.dot(X);#矩阵乘
XtXInv = np.linalg.inv(XtX)#求逆
XtXInvX = XtXInv.dot(Xt)
coef = XtXInvX.dot(y.T)#权重值
y_est = X.dot(coef)
return y_est,coef
def fit_fun(x):
'''
要拟合的函数
'''
return np.power(x,5)
#    return np.sin(x)
#    return 3+ 5*x
if __name__ == '__main__':
data_num = 100;#数据数量，也就是x个数
ploy_dim =10;#拟合参数个数，即权重数量
noise_scale = 0.2;#噪声参数
## 数据准备
x = np.array(np.linspace(-2*np.pi,2*np.pi,data_num))  #数据
y = fit_fun(x)+noise_scale*np.random.rand(1,data_num)  #添加噪声
# 最⼩⼆乘拟合
[y_est,coef] = leastseq_byploy(x,y,ploy_dim)
#显⽰拟合结果
org_data = plt.scatter(x,y,color="r",marker='o',s = 50)
est_data = plt.plot(x,y_est,color="g",linewidth= 3)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Fit funtion with leastseq method")
plt.legend(["Noise data","Fited function"]);
plt.show()
当然输出没什么不⼀样：
linspace函数python
其中，有个lpoy_ dim参数，当时并不懂为什么有这么个万⼀，后来想明⽩了，最⼩⼆乘法就是⽤⼀个多项式函数去逼近它，这⾥将它设为10，表⽰⽤9次的多项式来逼近对象。调整这个参数，发现结果还是蛮有意思的。
scipy库封装的最⼩⼆乘法使⽤
scipy在numpy基础上增加了许多数学计算。其中拟合和优化模块optimize提供了许多数值计算模块。⽽我们可以使⽤optimize.leastsq()函数来使⽤最⼩⼆乘法。该函数需要传⼊：计算误差的函数、待确定参数（权重）初始值。给个例⼦：
from scipy import optimize
import numpy as np
def fit_fun(x):
return 10+5*x
x = np.linspace(-2*np.pi,2*np.pi,100)
y = fit_fun(x)+0.2*np.random.rand(1,100)
y1 = y.reshape(100,)
def residuals(p,y,x): #
"计算以p为参数的直线和原始数据之间的误差"
k, b = p
return (k*x + b)-y
p_init = np.random.randn(2)  # 随机初始化多项式参数
r = optimize.leastsq(residuals,p_init,args=(y1, x))
k, b = r[0]
print("k =",k, "b =",b)
输出为：
k = 5.00274134694685 b = 10.09289512565644
官⽅⽂档： 。重要的有3个参数，也就是前三个。这⾥参数传递有些注意点:⾸先是第⼀个参数是残差函数；第⼆个是初始参数值，数据类型是数组；第三个是额外参数，所有额外参数都可以放⼊这个元组中，我传⼊了x与y的值，这⾥要注意的是，x和y是数组，并且⼆者之间数据的size必须⼀样。我将yreshape了⼀下，但注意reshape(100,)和reshape(100,1)是不⼀样的，后者会报错。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。