参考官网:Scipy.
 
最小二乘最小化问题
问题形式
SciPy能够解决有约束性非线性最小二乘法问题
 其中
    
     
      
       
        f
       
       
        (
       
       
        x
       
       
        )
       
      
      
       f(x)
      
     
    f(x)是残值,
    
     
      
       
        ρ
       
       
        (
       
       
        .
       
       
        )
       
      
      
       ρ(.)
      
     
    ρ(.)是损失函数。一个线性损失函数给出了一个标准的最小二乘法问题。此外,允许对一些
    
     
      
       
        
         x
        
        
         j
        
       
      
      
       x_j
      
     
    xj进行上下限的约束。所有针对最小二乘法的方法都利用了一个叫做雅可比的偏导矩阵 。强烈建议以解析式的方式计算这个矩阵,并将其传递给最小二乘法,否则,它将通过有限差分来估计,这需要大量的额外时间,而且在困难情况下可能非常不准确。
一个例子
考虑如下式子(酶促反应):
 
 
    
     
      
       
        
         y
        
        
         i
        
       
      
      
       y_i
      
     
    yi是测量值,
    
     
      
       
        
         u
        
        
         i
        
       
      
      
       u_i
      
     
    ui是独立变量。待拟合的未知系数为
    
     
      
       
        x
       
       
        =
       
       
        (
       
       
        
         x
        
        
         0
        
       
       
        ,
       
       
        
         x
        
        
         1
        
       
       
        ,
       
       
        
         x
        
        
         2
        
       
       
        ,
       
       
        
         x
        
        
         3
        
       
       
        
         )
        
        
         T
        
       
      
      
       x=(x_0,x_1,x_2,x_3)^T
      
     
    x=(x0,x1,x2,x3)T
 建议先计算雅可比矩阵:
 
 为了找到一个有物理意义的解决方案,避免潜在的除以零问题的出现,并确保收敛到全局最小值,我们施加了一些限制条件:
    
     
      
       
        0
       
       
        ≤
       
       
        
         x
        
        
         j
        
       
       
        ≤
       
       
        100
       
       
        ,
       
       
        j
       
       
        =
       
       
        0
       
       
        ,
       
       
        1
       
       
        ,
       
       
        2
       
       
        ,
       
       
        3
       
      
      
       0≤x_j≤100,j=0,1,2,3
      
     
    0≤xj≤100,j=0,1,2,3.
 下面的代码实现了最小二乘法的估计,最后绘制了原始数据和拟合模型函数。
import numpy as np
from scipy.optimize import least_squares
def model(x, u):
    return x[0] * (u ** 2 + x[1] * u) / (u ** 2 + x[2] * u + x[3])
def fun(x, u, y):
    return model(x, u) - y
def jac(x, u, y):
    J = np.empty((u.size, x.size))
    den = u ** 2 + x[2] * u + x[3]
    num = u ** 2 + x[1] * u
    J[:, 0] = num / den
    J[:, 1] = x[0] * u / den
    J[:, 2] = -x[0] * num * u / den ** 2
    J[:, 3] = -x[0] * num / den ** 2
    return J
#数据如下
u = np.array([4.0, 2.0, 1.0, 5.0e-1, 2.5e-1, 1.67e-1, 1.25e-1, 1.0e-1,
              8.33e-2, 7.14e-2, 6.25e-2])
y = np.array([1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2,
              4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2])
#设置初始值
x0 = np.array([2.5, 3.9, 4.15, 3.9])
res = least_squares(fun, x0, jac=jac, bounds=(0, 100), args=(u, y), verbose=1)
print(res.x)
 
其结果为:
`ftol` termination condition is satisfied.
Function evaluations 131, initial cost 4.4383e+00, final cost 1.5375e-04, first-order optimality 4.52e-08.
[0.192806   0.19130332 0.12306046 0.13607205]
 
绘制拟合图像
import matplotlib.pyplot as plt
u_test = np.linspace(0, 5)
y_test = model(res.x, u_test)
plt.plot(u, y, 'o', markersize=4, label='data')
plt.plot(u_test, y_test, label='fitted model')
plt.xlabel("u")
plt.ylabel("y")
plt.legend(loc='lower right')
plt.show()
 










