最小二乘法Python实现

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