在 Python 中将数据拟合到具有多个变量的互补误差函数

Posted 2023-03-24

【中文标题】在 Python 中将数据拟合到具有多个变量的互补误差函数

【英文标题】：Fitting data to a complementary error function with multiple variables in Python

【发布时间】：2021-05-14 06:34:28

【问题描述】：

我无法将实验数据拟合到 Python 3.7.4 中的互补误差函数。

import matplotlib.pyplot as plt
import math 
import numpy as np
from scipy import optimize
from scipy import integrate

with open('C:Path\\Data\\test.txt', 'r') as f:
    lines = f.readlines()
    x = [float(line.split(',')[0]) for line in lines]
    y = [float(line.split(',')[1]) for line in lines]


int_start = 35
int_end = 75

start = float(x[int_start])
end = float(x[int_end]) 
  
x_data = np.linspace(start, end, (int_end-int_start)+1)
y_data = y[int_start: int_end+1]


def integrand(x, a, b, c):
    return a*np.exp(((-1)*(x-b)**2)/(2*(c**2)))

def cerf(x, a, b, c):
    return integrate.quad(integrand, x, np.inf)

params, params_covariance = optimize.curve_fit(cerf, x_data, y_data)

plt.plot(x_data, y_data, 'x', label='Data')
plt.plot(x_data, integrand(x_data, params[0], params[1], params[2]), '-', label="fit")

plt.legend(loc='best')
plt.show()

更准确地说，我想将我的数据拟合到由带有参数a、b、c 的 integrand 函数和执行实际积分的 cerf 函数组成的互补误差函数。积分应该从 x（函数的参数）到 +infinity。之后，我想使用来自scipy 的标准curve_fit。但是现在我得到一个值错误如下：

> ValueError                                Traceback (most recent call last)
<ipython-input-33-8130a3eb44bb> in <module>
     29     return integrate.quad(integrand, x, np.inf)
     30 
---> 31 params, params_covariance = optimize.curve_fit(cerf, x_data, y_data)

~\AppData\Roaming\Python\Python37\site-packages\scipy\integrate\quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    346 
    347     # check the limits of integration: \int_a^b, expect a < b
--> 348     flip, a, b = b < a, min(a, b), max(a, b)
    349 
    350     if weight is None:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

如果有人知道如何使用 x 参数作为积分的下界来拟合函数，我将非常感激。

数据如下：

0.20,0.40
0.21,0.40
0.22,0.40
0.23,0.40
0.24,0.40
0.25,0.40
0.26,0.40
0.27,0.40
0.28,0.40
0.29,0.40
0.30,0.40
0.31,0.40
0.32,0.40
0.33,0.40
0.34,0.40
0.35,0.40
0.36,0.40
0.37,0.40
0.38,0.40
0.39,0.40
0.40,0.40
0.41,0.40
0.42,0.39
0.43,0.39
0.44,0.38
0.45,0.38
0.46,0.37
0.47,0.37
0.48,0.35
0.49,0.34
0.50,0.33
0.51,0.31
0.52,0.30
0.53,0.28
0.54,0.26
0.55,0.24
0.56,0.21
0.57,0.19
0.58,0.16
0.59,0.14
0.60,0.12
0.61,0.10
0.62,0.09
0.63,0.07
0.64,0.06 
0.65,0.05
0.66,0.04
0.67,0.03
0.68,0.02
0.69,0.02
0.70,0.01
0.71,0.01
0.72,0.00
0.73,0.00
0.74,0.00
0.75,0.00
0.76,0.00
0.77,0.00
0.78,-0.00
0.79,0.00
0.80,0.00
0.81,-0.00
0.82,-0.00
0.83,-0.00
0.84,0.00
0.85,-0.00
0.86,0.00
0.87,0.00
0.88,0.00
0.89,-0.00
0.90,0.00

【问题讨论】：

您能否在问题中添加一个具有代表性的数据示例。

【参考方案1】：

根据文档，scipy.integrate.quad 不接受数组，也不能调用带参数的函数。因此，我们必须在scipy.curve_fit 所寻址的函数中构造一个辅助函数f：

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize, integrate

#define a function that integrates or evaluates f depending on the Boolean flag func_integr
def cerf(x, a, b, c, func_integr=True):
    f = lambda x: a*np.exp(((-1)*(x-b)**2)/(2*(c**2)))
    #flag is preset to True, so will return the integrated values
    if func_integr:
        return np.asarray([integrate.quad(f, i, np.inf)[0] for i in x])
    #unless the flag func_integr is set to False, then it will return the function values
    else:
        return f(x)

#read file
arr=np.genfromtxt("test.txt", delimiter=",")
x_data = arr[:, 0]
y_data = arr[:, 1]

#provide reasonable start values...
start_p = [1, 0, -1]
#...for scipy.curve_fit
params, params_covariance = optimize.curve_fit(cerf, x_data, y_data, p0=start_p)
print(params)
#[ 2.26757666  0.56501062 -0.0704476 ]


#plot our results
plt.plot(x_data, y_data, 'x', label='Data')
plt.plot(x_data, cerf(x_data, *params), '-', label="fit")    
plt.legend(loc='best')
plt.show()

样本输出：

这种方法不是最快的 - 每个 x 值都是单独集成的。也许还有其他 scipy.integrate 函数可以与 numpy 数组一起使用；我不知道。 评估f 而不是它的积分值的部分在这里并不是真正需要的。但我最初使用它来验证cerf 是否按预期工作，所以我将它留在了脚本中。

【参考方案2】：

曲线的形状看起来更像逻辑而不是高斯类型。

【讨论】：

不，它应该是互补误差函数，因为当我用刀刃切割高斯光束时，我正在测量光电二极管上的电压。我想找出光束的浪费。这里也有：sfu.ca/~gchapman/e894/e894la1-extra4.pdf