使用修改后的 Bessel 函数进行 Gnuplot 拟合

Question

我想在 gnuplot 中拟合我的数据，拟合关系中包含第二类修改后的 Bessel 函数。所以假设它看起来像这样：

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 后，您可以在x = b 处找到y = x * besy1(x) 与B 的交集对应的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 访问您的函数，从我的其他答案（文件名test.py）修改Python 脚本：

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)。