使用修改的贝塞尔函数进行 Gnuplot 拟合

Question

我想将我的数据拟合到 gnuplot 中，并使用其中包含第二类修正贝塞尔函数的拟合关系。假设它看起来像这样:

f(x)= A*x + b*besselk(1,b)

(我用 matlab 或 Octave 音阶编写，我想通过拟合找到的系数是 A 和 b，所以其中一个在贝塞尔中)。但问题是 gnuplot 没有第二类修改的贝塞尔函数。

有人知道我该怎么做吗？

Answer 1

最佳拟合策略通常是将非线性或非多项式问题简化为线性或多项式问题。特别是，线性问题总是只有一个解。因此，我们理想地适合 f(x) = A*x + B 其中 B = b * besy1(b) - 这是第二类贝塞尔函数，请参阅在下面编辑第二类贝塞尔函数的修改，该函数在 Gnuplot 中不可用。你这样做:

fit A*x + B "datafile" via A, B

获得 B 后，您可以找到与 y = x * besy1(x) 与 的交集相对应的 b >B 位于 x = b。由于 besy1(x) 是振荡的，您可能会得到多个结果，但根据给定数据的范围，您可以选择正确的结果。假设您从拟合中得到 B = 1.2，则 [0:10] 区间中的交集如下:

plot [0:10] x*besy1(x), 1.2

如果您感兴趣的区域位于 x = 4.65 附近，其中有一个交点的大致位置，则查找确切的交点。在该区域中，x * besy1(x) 和 B 之间的距离将接近于零，因此距离的平方可以通过具有明确定义的最小值的抛物线来近似:

plot [4.6:4.7] (x*besy1(x)-1.2)**2

您的最佳x = b 就是该最小值的位置。您可以将其导出为数据并拟合抛物线 f(x) = a2*x**2 + b2*x + c2，最小值为 f'(x) = 0 ，即x = -b2/(2.*a2):

set table "data_minimum"
plot [4.6:4.7] (x*besy1(x)-1.2)**2
unset table
fit [4.6:4.7] a2*x**2 + b2*x + c2 "data_minimum" via a2,b2,c2
print -b2/2./a2

这给出了最小值位置的 x = 4.65447163370989，它对应于 B = b*besy1(b) 中的最佳 b。

其准确性将取决于二次拟合的好坏，而二次拟合的好坏又取决于 x 值的最小值范围有多窄。在本例中，范围 [4.6:4.7] 导致二次拟合相当好，但并不完美(您可以进一步缩小范围):

plot [4.6:4.7] "data_minimum" t "data", a*x**2+b*x+c t "quadratic fit"

编辑

对于第二类修改后的 Bessel 函数或 Gnuplot 中不可用的其他复杂函数，您可以使用外部解析器。例如，请参阅我关于如何使用外部 python 代码解析函数的答案:Passing Python functions to Gnuplot .

您可以使用scipy来访问您的函数，修改我的其他答案中的Python脚本(文件名test.py):

import sys
from scipy.special import kn as kn

n=float(sys.argv[1])
x=float(sys.argv[2])

print kn(n,x)

并在 Gnuplot 中使用它作为

kn(n,x) = real(system(sprintf("python test.py %g %g", n, x)))

那么上述所有过程只需将 besy1(x) 替换为 kn(1,x) 即可。