1

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

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))

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

scipy我对数据进行了硬编码，以消除任何数据依赖关系，并允许代码在任何已安装的机器上开箱即用。我不太明白该model01函数应该返回什么。根据文档，它应该返回一个数组。什么数组？在我的代码中，我假设我必须为每个数据点返回一个残差数组。那是对的吗？我的数据是一个二维数组，因为它是一个图像，但我的残差是一个展平的二维残差数组。那样行吗？有人能告诉我我到底做错了什么吗？有人可以修改和修复代码吗？scipy正如我上面提到的，代码应该在任何安装并numpy安装的机器上开箱即用。如果我想实现的目标是不可能的scipy.optimize.leastsq，您能否推荐一些其他适合使用 levenberg-marquardt 算法的模型的库？

4

1 回答 1

1

恐怕您无法使用leastsq. leastsq是一个围绕 FORTRAN 库的包装minipack，它调用MINPACK'slmdif和lmder算法。重要的是，它基于最小二乘目标函数的 Jacobian 和 Hessian。您的目标函数没有平滑导数，原因是：

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

和

math.fabs(......)

，因此leastsq将始终返回初始启动参数。

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

关于你的model101(). 它首先创建一个大小为width*height的一维数组g_data，称为data_r。然后它迭代g_data，计算math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r))每个元素并将该值放入一维数组data_r中。最后，它返回date_r。

您的解释是正确的，model101()为您的 2D 数据返回展平的 1D 残差数组。

于 2013-10-01T16:02:23.727 回答