scipy.optimize.leastsq无法适合简单模型

Question

我一直试图使用python的scipy.optimize库来估计模型的参数，但目前为止没有成功。我试图使用使用levenberg-marquardt算法的scipy.optimize.leastsq。不幸的是，即使我将初始参数的猜测设置得非常接近最佳拟合，它总是无法找到我的模型函数的最小值。实际上，它总是返回最初猜测的参数。所以，我认为我做错了什么。我的模型是一个简单的圆，但为了使事情更简单，只有半径是实际参数，数据中圆的中心是已知的并且是硬编码的。数据是一个10x10像素的浮点图像，中心为5,5和半径为4的圆。实际上，数据是使用我试图拟合的模型生成的。所以，存在一个完美契合。这里是我的代码：为了去除任何数据的依赖性，并允许代码安装scipy任何机器上运行开箱

import math 
import numpy 
import scipy.optimize 

# ======================================================================== 

g_data_width = 10 
g_data_height = 10 
g_xo = 5.0 
g_yo = 5.0 

def evaluate_model01(x,y,r): 
    x2 = x*x 
    y2 = y*y 
    r2 = r*r 
    v = 0.0 
    if(x2 + y2 <= r2): 
     v = 20.0 
    return v 

def model01(params,data_o): 
    data_r = numpy.zeros(g_data_height*g_data_width) 
    r = params[0] 
    for y in range(g_data_height): 
     for x in range(g_data_width): 
      xnew = x - g_xo 
      ynew = y - g_yo 
      data_r[y*g_data_width+x] = math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r)) 
    return data_r 

# ======================================================================== 

g_data_o = numpy.array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
         [ 0, 0, 0, 0, 0, 20, 0, 0, 0, 0], 
         [ 0, 0, 0, 20, 20, 20, 20, 20, 0, 0], 
         [ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0], 
         [ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0], 
         [ 0, 20, 20, 20, 20, 20, 20, 20, 20, 20], 
         [ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0], 
         [ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0], 
         [ 0, 0, 0, 20, 20, 20, 20, 20, 0, 0], 
         [ 0, 0, 0, 0, 0, 20, 0, 0, 0, 0]],dtype=numpy.float32) 
g_params = numpy.array([8.0]) 
print(scipy.optimize.leastsq(model01,g_params,args=(g_data_o),full_output=1)) 

# ======================================================================== 

我硬编码的数据。我不太明白的是model01函数应该返回。根据文件，它应该返回一个数组。一个什么样的数组？在我的代码中，我假设我必须为每个数据点返回一个残差数组。那是对的吗？我的数据是二维数组，因为它是一个图像，但我的残差是一个平坦的二维残差数组。这可以吗？有人能告诉我我究竟做错了什么吗？有人可以修改和修复代码吗？正如我上面提到的，代码应该在安装有scipy和numpy的任何机器上运行。如果scipy.optimize.leastsq不能实现我想要实现的功能，那么能否使用levenberg-marquardt算法推荐一些适合模型的其他库？

Answer 1

恐怕你不能用leastsq解决你的特殊问题。 leastsq是FORTRAN库minipack的环绕，它调用MINPACK的lmdif和lmder算法。重要的是，它基于最小二乘目标函数的雅可比矩阵和Hessian矩阵。你的目标函数不具有光滑的衍生物，由于：

if(x2 + y2 <= r2): 
    v = 20.0 

和

math.fabs(......) 

，所以leastsq总是返回初始启动参数。

您应该尝试使用一些不需要梯度/衍生的方法，例如Powell的方法fmin_powell或Nelder-Mead fmin。

关于您的model101()正在发生什么。它首先制作一个width*height的大小的一维数组，称为data_r。然后它遍历g_data，计算每个元素的math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r))并将该值放入一维数组data_r。最后，它返回date_r。

你的解释是正确的，model101()为你的二维数据返回平坦的1D残差阵列。