最小二乘法
数学和统计上面一个基本方法是,根据最小二衬发拟合平面上的点集。其拟合的图形通常是基本类型函数,如:线性函数、多项式、三角多项式等。由于数据有测量误差或者试验误差,我们不要求数据通过所有数据点。实际上,我们需要的是在所有数据点的y值,和逼近曲线相应点处的y值俩者之间误差的平方和最小意义下的最佳曲线。
具体原理推导不详细说,在线性代数的书里基本都会有介绍。下面介绍定理:
若A是秩为n的m×n的矩阵,则正规方程组:
A T A X = A T Y A^{T}AX = A^{T}Y ATAX=ATY
有唯一解:
x = ( A T A ) − 1 A T Y ( 1 ) x = (A^T A)^{-1}A^T Y (1) x=(ATA)−1ATY(1)
且x为方程组 AX = Y最小二乘解。
按书写习惯,X和Y都是列向量。参考最小二乘法求解 。由于输入数据格式不一样,刚开始看的时候没怎么看的懂。毕竟输入的是行向量,公式有点不同: ω = ( X X T ) − 1 X Y T ( 2 ) \omega = (XX^T)^{-1}XY^T (2) ω=(XXT)−1XYT(2)。这里生成的数据x和y都是行向量,所以使用的是公式(2),代码如下:
"""最小二乘法"""
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()
选择性注释掉 f i t f u n ( ) fit_fun() fitfun()函数中的返回值,选择不同的曲线作为被拟合对象。
然后按照我自己的理解,修改了其中的算法,才搞清楚。代码如下(注意运算按公式(1),列向量计算):
"""最小二乘法"""
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()
当然输出没什么不一样:
其中,有个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
官方文档:optimize.leastsq()用法 。重要的有3个参数,也就是前三个。这里参数传递有些注意点:首先是第一个参数是残差函数;第二个是初始参数值,数据类型是数组;第三个是额外参数,所有额外参数都可以放入这个元组中,我传入了x与y的值,这里要注意的是,x和y是数组,并且二者之间数据的size必须一样。我将yreshape了一下,但注意reshape(100,)和reshape(100,1)是不一样的,后者会报错。