如何在 Python 2.7 中同时优化和查找两个方程的系数？

Question

我有想要拟合两个方程的数据集：

y1 = a1 + a2 * T / 2 + a3 * T^2 / 3 + a4 * T^3 / 4 + a5 * T^4 / 5 + a6 / T
y2 = a1 * lnT + a2 * T + a3 * T^2 / 2 + a4 * T^3 / 3 + a5 * T^4 / 4 + a7

这两个多项式共享一些参数（a1 到 a5），所以我想同时拟合这两个方程。

我尝试用 scipy.optimize.curve_fit 来做到这一点：

import numpy as np
from scipy.optimize import curve_fit

def func(T, a1, a2, a3, a4, a5, a6, a7):
    y1 = a1 + a2 * T / 2 + a3 * T**2 / 3 + a4 * T**3 / 4 + a5 * T**4/5 + a6/T
    y2 = a1*np.log(T) + a2*T + a3 * T**2/2 + a4 * T**3/4 + a5 * T**4/4 + a7
    return np.stack((y1, y2), axis = 1)

T = np.linspace(300, 1000, 20)
ydata_1 = np.array([
    0.02139265,  0.40022353,  0.70653103,  0.95896469,  1.17025634,
    1.34944655,  1.50316659,  1.63641239,  1.75303086,  1.85603601,
    1.94782051,  2.03030092,  2.10501971,  2.17321829,  2.23589026,
    2.29382086,  2.34761661,  2.39772787,  2.44446625,  2.48801814])

ydata_2 = np.array([
    15.73868267,  16.14232408,  16.50633034,  16.83724622,
    17.14016153,  17.41914701,  17.67752993,  17.91807535,
    18.14310926,  18.35460465,  18.55424316,  18.74346017,
    18.92347836,  19.09533317,  19.25989235,  19.41787118,
    19.56984452,  19.71625632,  19.85742738,  19.99356154])

ydata = np.stack((ydata_1, ydata_2), axis = 1)
popt, pconv = curve_fit(f = func, xdata = T, ydata = ydata)

但是我得到了错误：

minpack.error: Result from function call is not a proper array of floats.

我什至不确定这是否是解决问题的正确方法。

Answer 1

您可以尝试在二维空间中为您的 y 值最小化 L_2 范数（即最小二乘拟合）：

from scipy.optimize import minimize

def func(params):
    a1, a2, a3, a4, a5, a6, a7 = params
    y1 = a1 + a2 * T / 2 + a3 * T**2 / 3 + a4 * T**3 / 4 + a5 * T**4/5 + a6/T
    y2 = a1*np.log(T) + a2*T + a3 * T**2/2 + a4 * T**3/4 + a5 * T**4/4 + a7
    return np.sum((y1 - ydata_1) ** 2 + (y2 - ydata_2) ** 2)

T = np.linspace(300, 1000, 20)
ydata_1 = np.array([
    0.02139265,  0.40022353,  0.70653103,  0.95896469,  1.17025634,
    1.34944655,  1.50316659,  1.63641239,  1.75303086,  1.85603601,
    1.94782051,  2.03030092,  2.10501971,  2.17321829,  2.23589026,
    2.29382086,  2.34761661,  2.39772787,  2.44446625,  2.48801814])

ydata_2 = np.array([
    15.73868267,  16.14232408,  16.50633034,  16.83724622,
    17.14016153,  17.41914701,  17.67752993,  17.91807535,
    18.14310926,  18.35460465,  18.55424316,  18.74346017,
    18.92347836,  19.09533317,  19.25989235,  19.41787118,
    19.56984452,  19.71625632,  19.85742738,  19.99356154])

# choose reasonable values for your 7 parameters here,
# i.e. close to the "right" answer, this may take a few tries
first_guess = [a1_0, a2_0, a3_0, a4_0, a5_0, a6_0, a7_0]  

# here we run the minimisation
res = minimize(func, first_guess)

# this is an array of your best fit values for a1-a7
best_fit = res.x

但是，@Stelios 似乎是正确的，因为您很难很好地适应您的特定模型。

Answer 2

（扩展评论）

a) curve_fit 用于拟合基于 single 数据集的 single 函数。在您的情况下，您有两个函数可以基于两个数据集进行拟合。这原则上需要从头开始设置优化问题，即定义一个单个目标函数（有或没有约束）。例如，目标函数可以是两次拟合的残差平方之和。然后，您将使用优化求解器（例如 scipy.optimimize.minimize）来找到最佳变量。

b) 您的模型（拟合函数）可能会在优化中引入数值困难。例如，y1 的变量a5 和a6 分别是T**4 和1/T 的因子，对于T=10**3，它们对应于10**12 和10**-3 的值。这是一个巨大的比例差异，接近硬件精度，这对我来说建议你应该重新考虑你的模型。