【发布时间】：2020-08-25 11:27:29
【问题描述】：
我在January 和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、xb1_exp、xb2_exp和xb3_exp的长度是@9876543354 @。
我的第二个问题是我不知道如何编写一个使用curve_fit 作为包含所有值的全局联合拟合的function。换句话说，我想为所有值找到fc 和alpha 的最佳通用拟合，以便我得到一个全局（或通用？）拟合而不是六个独立拟合。 objective提供test[0]=poptXB[0]-poptXA[0]，而test[1]是popXA和poptXB的范数，但returns都没有提供fc和alpha的拟合值，这正是我所寻求的。
这可能吗？
编辑 26.08.2020（和 30.08.2020）
我遇到了另一个question，关于关节配合并相应地调整了我的代码：
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 和alpha 提供的值：
[[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 的新手，我想知道它是否应该这样使用。我的做法是正确的还是完全错误的？
【问题讨论】：
标签： python scipy scipy-optimize lmfit