用边界可变的多个域函数拟合数据

Question

正如标题所提到的，我很难将数据点拟合到具有3个域的函数，这些域的边界是我的函数的参数。下面是我正在处理的函数：

global sigma_m
sigma_m=2*10**(-12)

global sigma_f
sigma_f=10**3

def Conductivity (phi,phi_c,t,s):
    sigma=[0]*(len(phi))
    for i in range (0,len(phi)):
        if phi[i]<phi_c:
            sigma[i]=sigma_m*(phi_c-phi[i])**(-s)
        elif phi[i]==phi_c:
            sigma[i]=sigma_f*(sigma_m/sigma_f)**(t/(t+s))
        else:
            sigma[i]=sigma_f*(phi[i]-phi_c)**t      
    return sigma 

我的数据点是：

phi_data=[0,0.005,0.007,0.008,0.017,0.05,0.085,0.10]
sigma_data=[2.00E-12,2.50E-12,3.00E-12,9.00E-04,1.00E-01,1.00E+00,2.00E+00,3.00E+00]

我的约束是phi_c、s和t必须严格大于0(在实践中，phi_c很少大于0.001，但高于0.5，s通常在0.5到1.5之间，t通常在1.5到6之间)。

我的目标是拟合我的数据点，并让我的拟合给出phi_c、s和t的值。可以估计s和t来帮助代码(在我显示的特定数据点集中，t应该在2左右，s应该在0.5左右)。除了我在上面提到的值范围之外，phi_c是完全未知的。

我使用了来自scipy的curve_fit和来自lmfit的模型，但它们都提供了非常小的phi_c值(例如10**(-16)或类似的小值，这让我相信程序希望phi_c为负值)。

下面是我使用curve_fit时的代码：

popt, pcov = curve_fit(Conductivity, phi_data, sigma_data, p0=[0.01,2,0.5], bounds=(0,[0.5,10,3]))

以下是我在使用lmfit中的Model时的代码：

t_estimate=0.5
s_estimate=2
phi_c_estimate=0.005
condmodel = Model(Conductivity)
params = condmodel.make_params(phi_c=phi_c_estimate,t=t_estimate,s=s_estimate)
result = condmodel.fit(sigma_data, params, phi=phi_data)
params['phi_c'].min = 0
params['phi_c'].max = 0.1

这两种选择在绘制时都可以很好地拟合，但phi_c的估计值远非可信。如果你有任何想法，我可以做更好的适合，请让我知道！

PS:我读了一篇很有希望的文章，关于使用symfit包分别拟合不同区域上的数据，不幸的是symfit包不适合我。它一直在卸载我的scipy版本，然后重新安装一个旧版本，然后它告诉我它需要一个新版本的scipy才能工作。

编辑:我设法让symfit包正常工作。下面是我的完整代码：

from symfit import parameters, variables, Fit, Piecewise, exp, Eq
import numpy as np
import matplotlib.pyplot as plt

global sigma_m
sigma_m=2*10**(-12)

global sigma_f
sigma_f=10**3

phi, sigma = variables ('phi, sigma')
t, s, phi_c = parameters('t, s, phi_c')

phi_c.min = 0.001
phi_c.max = 0.1

sigma1 = sigma_m*(phi_c-phi)**(-s)
sigma2 = sigma_f*(phi-phi_c)**t
model = {sigma: Piecewise ((sigma1, phi <= phi_c), (sigma2, phi > phi_c))}

constraints = [Eq(sigma1.subs({phi: phi_c}), sigma2.subs({phi: phi_c}))]

phi_data=np.array([0,0.005,0.007,0.008,0.017,0.05,0.085,0.10])
sigma_data=np.array([2.00E-12,2.50E-12,3.00E-12,9.00E-04,1.00E-01,1.00E+00,2.00E+00,3.00E+00])

fit = Fit(model, phi=phi_data, sigma=sigma_data, constraints=constraints)
fit_result = fit.execute()
print(fit_result)

不幸的是，我得到了以下错误：

  File "D:\Programs\Anaconda\lib\site-packages\sympy\printing\pycode.py", line 236, in _print_ComplexInfinity
    return self._print_NaN(expr)

  File "D:\Programs\Anaconda\lib\site-packages\sympy\printing\pycode.py", line 74, in _print_known_const
    known = self.known_constants[expr.__class__.__name__]

KeyError: 'ComplexInfinity'

我的编码知识非常有限，我不知道这意味着什么，我应该做什么才能不再有这个错误。如果你有什么想法，请告诉我。

Answer 1

我不确定我是否有一个单一的答案，但这将是太长了，不适合在评论中。

首先，切换函数形式的模型特别具有挑战性。但是，更重要的是，您的表单具有

 elif phi[i]==phi_c:

对于作为变量的浮点数，这基本上永远不会是真的。您的意思可能不是“完全相等”，而是“非常接近”，这可能是

 elif abs(phi[i] - phi_c) < 1.0e-5:  

或者别的什么..。

但是，将其从for循环转换为使用numpy.where()可能也是值得一试的。

其次，为了确保函数的连续性，你的不同形式实际上在边界上的值是相同的，这一点一点也不清楚。你可能会想要检查一下。

第三，具有幂和指数的模型特别具有挑战性，因为功率的微小变化可能会对结果价值产生巨大影响。也很容易得到“负值提升到非整数值”，这当然很复杂。

第四，那些sigma_m和sigma_f常量看起来很容易造成麻烦。你绝对应该用你的起始参数值来评估你的模型，看看你是否能用你的模型和合理的起始值来重现你的数据。我怀疑您需要更改您的起始值。