最小二乘法和高斯牛顿方程

2021-07-11

最小二乘法和高斯牛顿方程的Python实现

应用背景：消除误差或忽略无关细节，从干扰数据中提取信号或找出趋势，将大量数据降低到可管理的数量或用简单的近似来代替复杂函数。

这里我选用的是Python，需要用到numpy来做矩阵运算，matplotlib来绘制图像，math来引入自然底数e。

主要逻辑定义在GaussNetonMehthod中，因为需要在迭代体中计算雅可比矩阵函数和残差函数向量，所以将他们的参数作为GaussNetonMehthod的参数传入。

为了方便修改模型，我将雅可比矩阵函数和模型函数单独定义。又因为这两个函数结构简单，所以我使用lambda函数的形式，避免过多出现一句话的函数。

应用实例

根据美国1815年至1885年数据，估计人口模型中的参数A和B。如下表所示，已知年份和人口总量，及人口模型方程，求方程中的参数。

Year	Population
1815	8.3
1825	11.0
1835	14.7
1845	19.7
1855	26.7
1865	35.2
1875	44.4
1885	55.9

import matplotlib.pyplot as plt

import numpy as np

from numpy.linalg import inv

import math

from math import e

 

\# f = x1*exp(x2*t)

func = lambda x,t:x[0]*pow(e,x[1]*t)

\# J = dr/dx = -df/dx = [exp(x2*t) , t*x1*exp(x2*t)]

jaco = lambda x,t:[pow(e,x[1]*t),x[0]*t*pow(e,x[1]*t)]

 

def GaussNewtonMethod(t,y,x0,maxK,epsl):

  CalculatedPoints = x0

  for k in range(maxK):

    \#计算Jacobian矩阵

    Jacobian = []

    for i in t:

      Jacobian.append(jaco(CalculatedPoints,i))

    

    \#计算residual向量

    res_vec = []

    for i in range(len(y)):

      res_vec.append(func(CalculatedPoints,t[i]) - y[i])

 

    \#计算deltaX

    res_vec = np.array(res_vec)

    Jacobian = np.array(Jacobian)

    JacobianT = Jacobian.transpose()

    JTJ = JacobianT.dot(Jacobian)

    JTJ_inv = inv(JTJ)

    temp = JTJ_inv.dot(JacobianT)

    deltax = temp.dot(res_vec)

    

    \#计算||deataX||

    dX = np.linalg.norm(deltax,axis=0)

    if dX < epsl:

      break

    \#x(k+1)=x(k)-dX

    CalculatedPoints = CalculatedPoints-deltax

 

  print('经过{}次迭代取得最优解'.format(k+1))

  print('此时,\nX1 = {}\nX2 = {}'.format(CalculatedPoints[0],CalculatedPoints[1]))

  print('拟合曲线为 y = %.3f*exp(%.3f * t)'%(CalculatedPoints[0],CalculatedPoints[1]))

  \#print(Jacobian)

  fit = []

  for i in t:

    fit.append(func(CalculatedPoints,i))

  fig, axs = plt.subplots()

  axs.scatter(1805+10*t, y)

  axs.plot(1805+10*t, fit)

  fig.suptitle('Gauss - Newton')

  plt.xlabel('Year')

  plt.ylabel('United States Population(in millions)')

  plt.show()

 

def main():

  t = np.array([1,2,3,4,5,6,7,8]) #将1815年视作1，以此类推

  y = np.array([8.3,11.0,14.7,19.7,26.7,35.2,44.4,55.9])#人口数

  x0 = [6,0.3]

  maxK = 1000

  epsl = 1e-5

  GaussNewtonMethod(t,y,x0,maxK,epsl)

 

if __name__ == '__main__':

  main()