最小二乘法
数学和统计上面一个基本方法是，根据最小二衬发拟合平面上的点集。其拟合的图形通常是基本类型函数，如：线性函数、多项式、三角多项式等。由于数据有测量误差或者试验误差，我们不要求数据通过所有数据点。实际上，我们需要的是在所有数据点的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)是不一样的，后者会报错。
