Scipy 最小化、fmin、leastsq 类型问题（用序列设置数组元素）、不合适答案

【问题标题】：Scipy minimize, fmin, leastsq type problems (setting array element with sequence), bad fitScipy 最小化、fmin、leastsq 类型问题（用序列设置数组元素）、不合适
【发布时间】：2013-08-09 13:42:58
【问题描述】：

我在使用 scipy.optimize.fmin 和 scipy.optimize.minimize 函数时遇到问题。我检查并确认传递给函数的所有参数都是 numpy.array 类型，以及错误函数的返回值。此外，carreau 函数返回一个标量值。

一些额外参数（例如大小）的原因是：我需要使用给定模型（Carreau）拟合数据。数据是在不同的温度下获取的，并用偏移因子（也由模型拟合）进行校正，我最终得到了几组数据，这些数据都应该用于计算 same 4常数（参数 p）。

我读到我无法将数组列表传递给 fmin 函数，因此我必须将所有数据连接到 x_data_lin 中，并使用 size 参数跟踪不同的集合。 t 保存不同的测试温度，而 t_0 是一个包含参考温度的单元素数组。

我很肯定（三重检查）传递给函数的所有参数以及结果都是一维数组。除此之外的代码如下：

import numpy as np
import scipy.optimize
from scipy.optimize import fmin as simplex

def err_func2(p, x, y, t, t_0, size):
    result = array([])
    temp = 0
    for i in range(0, int(len(size)-1)):
        for j in range(int(temp), int(temp+size[i])):
            result = np.append(result, (carreau(p, x[j], t[i], t_0[0])-y[i]))
        temp += size[i]
    return result

p1 = simplex(err_func2, initial_guess,
            args=(x_data_lin, y_data_lin, t_list, t_0, size), full_output=0)

这是错误：

Traceback (most recent call last):
  File "C:\Python27\Scripts\projects\Carreau - WLF\carreau_model_fit.py", line 146, in <module>
    main()
  File "C:\Python27\Scripts\projects\Carreau - WLF\carreau_model_fit.py", line 105, in main
    args=(x_data_lin, y_data_lin, t_list, t_0, size), full_output=0)
  File "C:\Python27\lib\site-packages\scipy\optimize\optimize.py", line 351, in fmin
    res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
  File "C:\Python27\lib\site-packages\scipy\optimize\optimize.py", line 415, in _minimize_neldermead
    fsim[0] = func(x0)
ValueError: setting an array element with a sequence.

值得注意的是，在传递数组列表时，我得到了 leastsq 函数。不幸的是，它在拟合数据方面做得很差。但是，由于我花了很多时间和研究来达到这一点，我将发布代码如下。如果有人有兴趣查看所有代码，我很乐意将其发布，如果您可以推荐我上传一些文件（因为它包含另一个导入的脚本，当然还有示例数据）：

##def error_function(p, x, y, t, t_0):
##    result = array([])
##    for index in range(len(x)):
##        result = np.append(result,(carreau(p, x[index],
##                                           t[index], t_0) - y[index]))
##    return result

##    p1, success = scipy.optimize.leastsq(error_function, initial_guess,
##                                         args=(x_list, y_list, t_list, t_0),
##                                         maxfev=10000)

:(我打算发布一张最小平方拟合的图形数据的图片，但我没有必要的 10 分。

后期编辑：我现在已经让 optimize.curvefit 和 optimize.leastsq 工作（可能不太巧合地给出相同的答案），但曲线很糟糕。我一直在试图弄清楚optimize.minimize，但这有点让人头疼。单纯形（fmin，Nelder_Mead，无论你想怎么称呼它）都会运行，但会产生一个疯狂的答案。我以前从未处理过非线性优化问题，我真的不知道该往哪个方向发展。

这是有效的curve_fit代码：

def temp_shift(t_s, t, t_0):
    """ This function calculates the a_t temperature shift factor for polymer
    viscosity curves. Variable is the standard temperature, t_s
    """
    C_1 = 8.86
    C_2 = 101.6
    return(np.exp(
        (C_1*(t_0-t_s) / (C_2+(t_0-t_s))) - (C_1*(t-t_s) / (C_2 + (t-t_s)))
        ))

def pass_data(t, t_0):
    def carreau_2(x, p0, p1, p2, p3):
        visc_0 = p0
        m = p1
        n = p2
        t_s = p3
        a_T = temp_shift(p3, t, t_0)
        return (visc_0 * a_T / (1 + m * x * a_T)**n)
    return carreau_2

initial_guess = array([20000, 3, 0.94, -20])

p1, conv = scipy.optimize.curve_fit(pass_data(t_all, t_0), x_data_lin,
                                   y_data_lin, initial_guess)

以下是一些示例数据：

x_data_lin = array([0.01998, 0.04304, 0.2004, 0.43160, 0.92870, 2.0000, 4.30900,
                    9.28500, 15.51954, 21.94936, 37.52960, 90.41786, 204.35230,
                    331.58495, 811.92250, 1694.55309, 3464.27648, 8826.65738,
                    14008.00242])   

y_data_lin = array([13520.00000, 13740.00000, 12540.00000, 9384.00000, 5201,
                    3232.00000, 2094.00000, 1484.00000, 999.00000, 1162.05088
                    942.56946, 705.62489, 429.47341, 254.15136, 185.22916, 
                    122.07113, 76.46324, 47.85064, 25.74315, 18.84875])

t_all = array([190, 190, 190, 190, 190, 190, 190, 190, 190, 190, 190, 190, 
               190, 190, 190, 190, 190, 190, 190])

t_0 = 80

这是curve_fit结果的图片（现在我有10个点并且可以发布！）。请注意，绘制了 3 条曲线，因为我在 3 个不同的温度下使用了 3 组数据来优化曲线。聚合物具有剪切速率 - 粘度关系保持不变的特性，只是移动了一个温度因子 a_T：

我非常感谢有关如何改善拟合，或如何定义函数以使 optimize.minimize 起作用以及哪种方法（Nelder-Mead、Powel、BFGS）可能起作用的任何建议。

要添加的另一个编辑：我让 Nelder-Mead（optimize.fmin，和默认的 optimize.minimize）函数工作 - 我将在下面包含修改后的错误函数。之前，我只是简单地对结果数组求和并返回它。这导致了极负值（显然，因为函数的目标是最小化）。在求和之前对结果求平方解决了这个问题。请注意，按照 JaminSore 的建议，我还完全更改了函数以利用 numpy 的数组广播（感谢 Jamin！）

def err_func2(p, x, y, t, t_0):
    return ((carreau(p, x, t, t_0)-y)**2).sum()

不幸的是，Nelder-Mead 函数给出的结果与 leastsq 和 curve_fit 相同。您可以在上图中看到它不是最佳拟合；事实上，在这一点上，Microsoft Excel 的求解器函数在数据上做得更好。

至少，我希望这个线程可以对初学者在未来进行 scipy.optimize 有用，因为我花了很长时间才发现所有这些。

【问题讨论】：

标签： python numpy scipy

【解决方案1】：

与leastsq 不同，fmin 只能处理返回标量的错误函数，因此如果可能，您必须重写错误函数以使其返回标量。这是一个简单的工作示例。

导入必要的库

import numpy as np
from scipy.optimize import fmin

定义一个辅助函数（稍后会看到）

def prob(a, b):
    return (1 + np.exp(b - a))**-1

模拟一些数据

true_ = np.random.normal(size = 100) #parameters we're trying to recover
b = np.random.normal(size = 20)

exp_ = prob(true_[:, None], b) #expected
a_s, b_s = true_.shape[0], b.shape[0]
noise = np.random.uniform(size = (a_s, b_s))
response = (noise > (1 - exp_)).astype(int)

定义我们的错误函数（我使用的是lambdas，但实际上不建议这样做）

# sum of the squared residuals
err_func = lambda a : ((prob(a[:, None], b) - response) ** 2).sum()
result = fmin(err_func, np.zeros_like(true_)) #solve

如果我在错误函数定义的末尾删除.sum()，我会得到同样的错误。

【讨论】：

我以为只有minimize 函数需要标量误差函数。至少，文档对我来说似乎并不那么清楚。 optimize.fmin 的文档没有说标量，而optimize.minimize 的文档却说。无论如何，谢谢你的提示！现在我只需要弄清楚如何将这个算法变成一个标量函数。
没问题。要将您的函数转换为标量，简单的return (result**2).sum()。此外，for i in range(len(x)) 几乎总是不必要的，就像循环遍历 numpy 数组列表一样。使用 numpy 的broadcasting rules 可能会大大简化和加快您的问题。如果您发现您必须循环遍历数组列表，请查看this video 了解更多 Pythonic 方法。
如果你能够向量化你的函数，因为你试图估计 4 个常数，curve_fit 可能是要走的路。
我试过curve_fit，但它给了我与optimize.leastsq相同的结果（见上面的编辑）。

【解决方案2】：

好的，现在我终于知道答案了！首先，最后一块，然后是回顾。拟合问题不是curve_fit、leastsq、Nelder_Mead 或Powell（我尝试过的方法）的错。它与误差的相对权重有关。由于该数据是对数尺度的，因此高 y 值附近的拟合误差非常昂贵，而低 y 值附近的误差微不足道。为了纠正这个问题，我通过除以数据的 y 值来使错误相对，如下所示：

def err_func2(p, x, y, t, t_0):
    return (((carreau(p, x, t, t_0)-y)/y)**2).sum()

现在，对每个相对误差进行平方、求和，然后最小化，得到以下拟合（将 optimize.minimize 与 Powell 方法一起使用，尽管其他方法也应该相同。）

所以现在回顾一下在这个线程中找到的答案：

处理曲线拟合的最简单方法（或者至少对我来说，最简单）是将所有数据收集到 1D numpy.arrays 中。然后，您可以依靠 numpy 的数组广播来执行所有操作。这意味着算术运算的处理方式与向量点积相同。例如，array_1 = [a,b]，array_2 = [c,d]，然后 array_1 + array_2 = [a+c, b+d]。这适用于加法、减法、乘法、除法和幂：array+1array_2 = [ac, b**d]。
对于optimize.leastsq函数，需要让目标函数返回一个数组；即return result，其中结果是一个数组。对于 optimize.curve_fit，您还返回一个数组。在这种情况下，传递额外的参数（想想其他常量）有点复杂，但您可以使用嵌套函数来完成，正如我在上面的 pass_data 函数中演示的那样。
对于 optimize.minimize，您需要返回一个标量 - 即单个数字。我认为，您还可以返回一组答案，但正如我之前提到的，我通过将所有数据放入一维数组来避免这种情况。要获得这个标量，您可以简单地对结果进行平方和求和（就像我在err_func2 下的这篇文章中所写的那样）对数据进行平方非常重要，否则负错误会接管并导致结果标量非常负。
最后，如前所述，当您的数据跨越多个尺度（105、104、10**3 等）时，可能需要对错误进行归一化。我通过将每个错误除以 y 值来做到这一点。

所以...我猜是这样？最后？

【讨论】：

更常见的是通过除以方差来归一化误差。
好吧，在这种情况下，平均值不是很有意义（双关语），也无助于解决问题。考虑预测值 10000 和 100，实际值 9800 和 98，平均值为 5000。误差分别为 200 和 2。方差为 2.35E7。权重分别为 200/2.35E7 和 2/2.35E7。您可以看到，数字越大，权重仍然越大（显着性的 100 倍）。
另一方面，200/9800 和 2/98 的权重相同。因此，除以实际数据值应该使误差在整个范围内具有相同的权重。就像我提到的，数据跨越多个尺度这一事实确实会影响误差计算。