使用 curve_fit 与多个方程和输入进行联合拟合 - 这甚至可能吗？
Joint-fit using curve_fit with multiple equations and inputs - Is it even possible?
我在 and April 中问了类似的问题，@Miłosz Wieczór 和 @Joe 很友好地表现出了兴趣。现在，我面临着一个类似但不同的问题，因为我需要与多个方程和输入以获得两个参数 fc 和 alpha 的通用解。我的代码（基于前面问题的答案）如下：
import numpy as np
from numpy import linalg
import math
from scipy.optimize import curve_fit, least_squares, minimize
ya_exp = np.array([1, 1.3, 1.7, 2.1, 2.7, 3.5, 4.5, 5.8, 7.5, 9.7, 12, 16, 21, 27, 34, 44, 57, 73, 94, 120, 156, 250000])
yb_exp = np.array([1, 1.3, 1.7, 2.1, 2.7, 3.5, 4.5, 5.8, 7.5, 9.7, 12, 16, 21, 27, 34, 44, 57, 73, 94, 120, 156])
xa1_exp = np.array([4.68, 4.70, 4.71, 4.72, 4.74, 4.75, 4.76, 4.77, 4.79, 4.80,
4.82, 4.83, 4.85, 4.87, 4.89, 4.90, 4.96, 4.99, 5.02, 5.06,
5.11, 6.23])
xb1_exp = np.array([0.018, 0.023, 0.019, 0.023, 0.022, 0.023, 0.023, 0.023, 0.023, 0.024,
0.025, 0.028, 0.032, 0.033, 0.034, 0.037, 0.040, 0.043, 0.045, 0.047,
0.049])
xa2_exp = np.array([7.01, 7.03, 7.04, 7.04, 7.04, 7.10, 7.13, 7.13, 7.16, 7.14,
7.19, 7.18, 7.19, 7.22, 7.24, 7.28, 7.32, 7.35, 7.40, 7.45,
7.49, 10.1])
xb2_exp = np.array([0.008, 0.009, 0.008, 0.009, 0.008, 0.010, 0.010, 0.010, 0.011, 0.012,
0.016, 0.017, 0.020, 0.023, 0.027, 0.029, 0.036, 0.040, 0.043, 0.046,
0.052])
xa3_exp = np.array([5.67, 5.67, 5.68, 5.69, 5.72, 5.74, 5.74, 5.76, 5.76, 5.79,
5.81, 5.81, 5.83, 4.86, 5.89, 5.91, 5.96, 6.00, 6.04, 6.10,
6.14, 7.56])
xb3_exp = np.array([0.011, 0.011, 0.012, 0.011, 0.012, 0.012, 0.012, 0.012, 0.015, 0.017,
0.021, 0.026, 0.028, 0.030, 0.036, 0.039, 0.046, 0.050, 0.056, 0.059,
0.063])
xa1_zero = np.min(xa1_exp)
xa1_inf = np.max(xa1_exp)
xa2_zero = np.min(xa2_exp)
xa2_inf = np.min(xa2_exp)
xa3_zero = np.min(xa3_exp)
xa3_inf = np.min(xa3_exp)
ig_fc = 500
ig_alpha = 0.35
def CCXA(f_exp, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
RI = np.sqrt(R ** 2 + I ** 2)
return RI
def CCXB(f_exp, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
iQ = I / R
return iQ
poptXA, pcovXA = curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha))
poptXB, pcovXB = curve_fit(CCXB, yb_exp, xb_exp, p0=(ig_fc, ig_alpha))
def objective(e_exp, f_exp):
poptXA, pcovXA = curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha))
poptXB, pcovXB = curve_fit(CCXB, yb_exp, xb_exp, p0=(ig_fc, ig_alpha))
err_total = np.sum(np.sqrt(np.diag(pcovXA))) + np.sum(np.sqrt(np.diag(pcovXB)))
delta = linalg.norm(poptXB - poptXA)
return err_total, delta
test = objective(xa_exp, ya_exp)
我的第一个问题是我不确定如何让 CCXA 和 CCXB 在全局范围内搜索和定位 xa_exp 和 xb_exp，因为不同的变量由其他名称定义：xa1_exp、xa2_exp 和 xa3_exp 以及 xb1_exp、xb2_exp 和 xb2_exp。同样，我对 fc 和 alpha 作为可优化参数感兴趣。 curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha)) 能够针对 fc 和 alpha 进行优化，但依赖于全局范围内的 xa_exp 和 xb_exp。当它们不同时如何将它们传递给函数？另外，请注意 ya_exp、xa1_exp、xa2_exp 和 xa3_exp 的长度是 22，而 yb_exp、[=23 的长度=]、xb2_exp 和 xb3_exp 是 21。
我的第二个问题是我不知道如何编写一个 function，它使用 curve_fit 作为包含所有值的全局联合拟合。换句话说，我想为所有值找到 fc 和 alpha 的最佳通用拟合，以便我得到一个全局（或通用？）拟合而不是六个独立拟合。 objective提供test[0]=poptXB[0]-poptXA[0]，而test[1]是popXA和poptXB的范数，但returns都没有提供[=14=的拟合值]和alpha，这就是我要找的。
这可能吗？
编辑 2020 年 8 月 26 日（和 2020 年 8 月 30 日）
我遇到了另一个关于关节配合的 并相应地调整了我的代码：
from lmfit import minimize, Parameters, fit_report
params = Parameters()
params.add('fc', value=500)
params.add('alpha', value=0.2)
def CCXA(f_exp, xa_zero, xa_inf, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
RI = np.sqrt(R ** 2 + I ** 2)
return RI
def CCXB(f_exp, xa_zero, xa_inf, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
iQ = I / R
return iQ
def fit_function(params, ya_data=None, xa1_data=None, xa2_data=None, xa3_data=None, yb_data=None, xb1_data=None, xb2_data=None, xb3_data=None):
xa_data = np.array([xa1_data, xa2_data, xa3_data])
modxa = np.array([CCXA(ya_data, np.min(i), np.max(i), params['fc'], params['alpha']) for i in xa_data])
#sigma = np.array([1 / np.sqrt(i+1) for i in xa_data])
#sigma = np.full((1,22), 0.5)
#resxa = np.array([(i - j)/k for i, j, k in zip(xa_data, modxa, sigma)])
resxa = np.array([i - j for i, j in zip(xa_data, modxa)])
xb_data = np.array([xb1_data, xb2_data, xb3_data])
modxb = np.array([CCXB(yb_data, np.min(i), np.max(i), params['fc'], params['alpha']) for i in xa_data])
#sigmb = np.array([1 / np.sqrt(i+1) for i in xb_data])
#sigmb = np.full((1,21), 1)
#resxb = np.array([(i - j)/k for i, j, k in zip(xb_data, modxb, sigmb)])
resxb = np.array([i - j for i, j in zip(xb_data, modxb)])
return np.concatenate((resxa.ravel(), resxb.ravel()))
对 CCXA 和 CCXB 的细微修改加上 fit_function 的引入，使用 lmfit 为 fc 和 [=15= 提供的值解决]:
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 23
# data points = 129
# variables = 2
chi-square = 8.89790022
reduced chi-square = 0.07006221
Akaike info crit = -340.945624
Bayesian info crit = -335.225999
[[Variables]]
fc: 1149.59953 +/- 572.031178 (49.76%) (init = 500)
alpha: 0.64118507 +/- 0.04236375 (6.61%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
C(fc, alpha) = 0.874
由于我是 lmfit 的新手，我想知道它是否应该这样使用。我做的是正确的还是完全错误的？
直接使用least_squares可能更容易。它需要一个残差向量，您只需在其中指定 l.h.s。你的两个方程
是的，您对 lmfit 和 np.concatenate() 的使用在我看来是正确的，它可以找到能够最小化 (model1, dataset1) 和 (model2, dataset2) 残差的一组参数值。我不确定我是否理解您的模型和 3 x 值。我认为这可能会更有效地完成，但这是一个不同的问题，即装配机制是否按照您的意图进行 - 我认为它们是。
下一个（更棘手的）问题是您是否要强调一个数据集中的失配比另一个更严重。如所写，您的函数断言所有数据点都具有相同的权重 - 这是一种很好的默认方法，但情况可能并非总是如此。
我在 fc 和 alpha 的通用解。我的代码（基于前面问题的答案）如下：
import numpy as np
from numpy import linalg
import math
from scipy.optimize import curve_fit, least_squares, minimize
ya_exp = np.array([1, 1.3, 1.7, 2.1, 2.7, 3.5, 4.5, 5.8, 7.5, 9.7, 12, 16, 21, 27, 34, 44, 57, 73, 94, 120, 156, 250000])
yb_exp = np.array([1, 1.3, 1.7, 2.1, 2.7, 3.5, 4.5, 5.8, 7.5, 9.7, 12, 16, 21, 27, 34, 44, 57, 73, 94, 120, 156])
xa1_exp = np.array([4.68, 4.70, 4.71, 4.72, 4.74, 4.75, 4.76, 4.77, 4.79, 4.80,
4.82, 4.83, 4.85, 4.87, 4.89, 4.90, 4.96, 4.99, 5.02, 5.06,
5.11, 6.23])
xb1_exp = np.array([0.018, 0.023, 0.019, 0.023, 0.022, 0.023, 0.023, 0.023, 0.023, 0.024,
0.025, 0.028, 0.032, 0.033, 0.034, 0.037, 0.040, 0.043, 0.045, 0.047,
0.049])
xa2_exp = np.array([7.01, 7.03, 7.04, 7.04, 7.04, 7.10, 7.13, 7.13, 7.16, 7.14,
7.19, 7.18, 7.19, 7.22, 7.24, 7.28, 7.32, 7.35, 7.40, 7.45,
7.49, 10.1])
xb2_exp = np.array([0.008, 0.009, 0.008, 0.009, 0.008, 0.010, 0.010, 0.010, 0.011, 0.012,
0.016, 0.017, 0.020, 0.023, 0.027, 0.029, 0.036, 0.040, 0.043, 0.046,
0.052])
xa3_exp = np.array([5.67, 5.67, 5.68, 5.69, 5.72, 5.74, 5.74, 5.76, 5.76, 5.79,
5.81, 5.81, 5.83, 4.86, 5.89, 5.91, 5.96, 6.00, 6.04, 6.10,
6.14, 7.56])
xb3_exp = np.array([0.011, 0.011, 0.012, 0.011, 0.012, 0.012, 0.012, 0.012, 0.015, 0.017,
0.021, 0.026, 0.028, 0.030, 0.036, 0.039, 0.046, 0.050, 0.056, 0.059,
0.063])
xa1_zero = np.min(xa1_exp)
xa1_inf = np.max(xa1_exp)
xa2_zero = np.min(xa2_exp)
xa2_inf = np.min(xa2_exp)
xa3_zero = np.min(xa3_exp)
xa3_inf = np.min(xa3_exp)
ig_fc = 500
ig_alpha = 0.35
def CCXA(f_exp, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
RI = np.sqrt(R ** 2 + I ** 2)
return RI
def CCXB(f_exp, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
iQ = I / R
return iQ
poptXA, pcovXA = curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha))
poptXB, pcovXB = curve_fit(CCXB, yb_exp, xb_exp, p0=(ig_fc, ig_alpha))
def objective(e_exp, f_exp):
poptXA, pcovXA = curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha))
poptXB, pcovXB = curve_fit(CCXB, yb_exp, xb_exp, p0=(ig_fc, ig_alpha))
err_total = np.sum(np.sqrt(np.diag(pcovXA))) + np.sum(np.sqrt(np.diag(pcovXB)))
delta = linalg.norm(poptXB - poptXA)
return err_total, delta
test = objective(xa_exp, ya_exp)
我的第一个问题是我不确定如何让 CCXA 和 CCXB 在全局范围内搜索和定位 xa_exp 和 xb_exp，因为不同的变量由其他名称定义：xa1_exp、xa2_exp 和 xa3_exp 以及 xb1_exp、xb2_exp 和 xb2_exp。同样，我对 fc 和 alpha 作为可优化参数感兴趣。 curve_fit(CCXA, ya_exp, xa_exp, p0=(ig_fc, ig_alpha)) 能够针对 fc 和 alpha 进行优化，但依赖于全局范围内的 xa_exp 和 xb_exp。当它们不同时如何将它们传递给函数？另外，请注意 ya_exp、xa1_exp、xa2_exp 和 xa3_exp 的长度是 22，而 yb_exp、[=23 的长度=]、xb2_exp 和 xb3_exp 是 21。
我的第二个问题是我不知道如何编写一个 function，它使用 curve_fit 作为包含所有值的全局联合拟合。换句话说，我想为所有值找到 fc 和 alpha 的最佳通用拟合，以便我得到一个全局（或通用？）拟合而不是六个独立拟合。 objective提供test[0]=poptXB[0]-poptXA[0]，而test[1]是popXA和poptXB的范数，但returns都没有提供[=14=的拟合值]和alpha，这就是我要找的。
这可能吗？
编辑 2020 年 8 月 26 日（和 2020 年 8 月 30 日）
我遇到了另一个关于关节配合的
from lmfit import minimize, Parameters, fit_report
params = Parameters()
params.add('fc', value=500)
params.add('alpha', value=0.2)
def CCXA(f_exp, xa_zero, xa_inf, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
RI = np.sqrt(R ** 2 + I ** 2)
return RI
def CCXB(f_exp, xa_zero, xa_inf, fc, alpha):
x = np.log(f_exp/fc)
R = xa_zero + 1/2 * (xa_inf - xa_zero) * (1 + np.sinh((1 - alpha) * x) / (np.cosh((1 - alpha) * x) + np.sin(1/2 * alpha * math.pi)))
I = 1/2 * (xa_inf - xa_zero) * np.cos(alpha * math.pi / 2) / (np.cosh((1 - alpha) * x) + np.sin(alpha * math.pi / 2))
iQ = I / R
return iQ
def fit_function(params, ya_data=None, xa1_data=None, xa2_data=None, xa3_data=None, yb_data=None, xb1_data=None, xb2_data=None, xb3_data=None):
xa_data = np.array([xa1_data, xa2_data, xa3_data])
modxa = np.array([CCXA(ya_data, np.min(i), np.max(i), params['fc'], params['alpha']) for i in xa_data])
#sigma = np.array([1 / np.sqrt(i+1) for i in xa_data])
#sigma = np.full((1,22), 0.5)
#resxa = np.array([(i - j)/k for i, j, k in zip(xa_data, modxa, sigma)])
resxa = np.array([i - j for i, j in zip(xa_data, modxa)])
xb_data = np.array([xb1_data, xb2_data, xb3_data])
modxb = np.array([CCXB(yb_data, np.min(i), np.max(i), params['fc'], params['alpha']) for i in xa_data])
#sigmb = np.array([1 / np.sqrt(i+1) for i in xb_data])
#sigmb = np.full((1,21), 1)
#resxb = np.array([(i - j)/k for i, j, k in zip(xb_data, modxb, sigmb)])
resxb = np.array([i - j for i, j in zip(xb_data, modxb)])
return np.concatenate((resxa.ravel(), resxb.ravel()))
对 CCXA 和 CCXB 的细微修改加上 fit_function 的引入，使用 lmfit 为 fc 和 [=15= 提供的值解决]:
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 23
# data points = 129
# variables = 2
chi-square = 8.89790022
reduced chi-square = 0.07006221
Akaike info crit = -340.945624
Bayesian info crit = -335.225999
[[Variables]]
fc: 1149.59953 +/- 572.031178 (49.76%) (init = 500)
alpha: 0.64118507 +/- 0.04236375 (6.61%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
C(fc, alpha) = 0.874
由于我是 lmfit 的新手，我想知道它是否应该这样使用。我做的是正确的还是完全错误的？
直接使用least_squares可能更容易。它需要一个残差向量，您只需在其中指定 l.h.s。你的两个方程
是的，您对 lmfit 和 np.concatenate() 的使用在我看来是正确的，它可以找到能够最小化 (model1, dataset1) 和 (model2, dataset2) 残差的一组参数值。我不确定我是否理解您的模型和 3 x 值。我认为这可能会更有效地完成，但这是一个不同的问题，即装配机制是否按照您的意图进行 - 我认为它们是。
下一个（更棘手的）问题是您是否要强调一个数据集中的失配比另一个更严重。如所写，您的函数断言所有数据点都具有相同的权重 - 这是一种很好的默认方法，但情况可能并非总是如此。