我想将我的数据拟合到 gnuplot 中,拟合关系中包含第二类修改后的 Bessel 函数。所以让我们说它看起来像这样:
f(x)= A*x + b*besselk(1,b)
(我是用 matlab 或 octave 编写的,而我想通过拟合找到的系数是 A 和 b,所以其中一个在贝塞尔中)。但问题是 gnuplot 没有修改过的第二类贝塞尔函数。
有人知道我该怎么做吗?
我想将我的数据拟合到 gnuplot 中,拟合关系中包含第二类修改后的 Bessel 函数。所以让我们说它看起来像这样:
f(x)= A*x + b*besselk(1,b)
(我是用 matlab 或 octave 编写的,而我想通过拟合找到的系数是 A 和 b,所以其中一个在贝塞尔中)。但问题是 gnuplot 没有修改过的第二类贝塞尔函数。
有人知道我该怎么做吗?
最佳拟合策略通常是将您的非线性或非多项式问题简化为线性或多项式问题。特别是,线性问题总是只有一个解。所以我们最适合f(x) = A*x + B
在哪里B = b * besy1(b)
- 这是第二类贝塞尔函数,请参阅下面的编辑修改第二类贝塞尔函数,这在 Gnuplot 中不可用。你这样做:
fit A*x + B "datafile" via A, B
拥有B
后,您可以找到与atb
的交集对应的。因为是振荡的,所以您可能有多个结果,但根据您的数据提供的范围,您可以选择正确的一个。假设你从拟合中得到,那么区间内的交点如下:y = x * besy1(x)
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 = b*besy1(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 代码解析函数的答案:将Python 函数传递给 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)
.