问题描述
我正在使用 python 批量处理一些数据并绘制它。我可以使用 scipy.curve_fit、一个双指数函数和一些合理的初始猜测很好地拟合它。这是一个代码片段:
def biexpfunc(x,a,b,c,d,e):
y_new = []
for i in range(len(x)):
y = (a * np.exp(b*x[i])) + (c * np.exp(d*x[i])) + e
y_new.append(y)
return y_new
x = np.linspace(0,160,100)
y = biexpfunc(x,50,-0.2,-0.1,10)
jitter_y = y + 0.5 *np.random.rand(len(y)) - 0.1
plt.scatter(x,jitter_y)
sigma = np.ones(len(x))
sigma[[0,-1]] = 0.01
popt,pcov = curve_fit(biexpfunc,x,jitter_y,p0 = (50,10),sigma = sigma)
x_fit = np.linspace(0,x[-1])
y_fit = biexpfunc(x_fit,*popt)
plt.plot(x_fit,y_fit,'r--')
plt.show()
我知道如何对给定的 x 值进行插值以找到 y(通过将其放回函数中),但是如何为给定的 y 值找到 x?我觉得必须有一种不需要重新安排和定义新函数的明智方法(部分原因是数学不是我的强项,我不知道如何!)。如果曲线很好地拟合了数据,有没有办法简单地读出一个值?任何帮助将不胜感激!
解决方法
事实证明,您的问题与曲线拟合无关,但实际上与求根有关。 Scipy.optimize
具有此任务的 whole arsenal of functions。选择和配置正确的有时很困难。我可能不是这里最好的向导,但由于没有其他人加强...
根查找尝试确定 x
为零的 f(x)
值。要找到一个 x0
,其中 f(x0)
是某个 y0
值,我们只需将函数转换为 g(x) = f(x)-y0
。
由于您的函数是单调的,对于给定的 y 值,预计不会超过一个根。我们也知道要搜索的 x 区间,所以 bisect
似乎是一个合理的策略:
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit,bisect
def biexpfunc(x,a,b,c,d,e):
return (a * np.exp(b*x)) + (c * np.exp(d*x)) + e
np.random.seed(123)
x = np.linspace(0,160,100)
y = biexpfunc(x,50,-0.2,-0.1,10)
jitter_y = y + 0.5 *np.random.rand(len(y)) - 0.1
fig,ax = plt.subplots(figsize=(10,8))
ax.scatter(x,jitter_y,marker="x",color="blue",label="raw data")
#your curve fit routine
sigma = np.ones(len(x))
sigma[[0,-1]] = 0.01
popt,pcov = curve_fit(biexpfunc,x,p0 = (50,10),sigma = sigma)
x_fit = np.linspace(x.min(),x.max(),100)
y_fit = biexpfunc(x_fit,*popt)
ax.plot(x_fit,y_fit,'r--',label="fit")
#y-value for which we want to determine the x-value(s)
y_test=55
test_popt = popt.copy()
test_popt[-1] -= y_test
#here,the bisect method tries to establish the x for which f(x)=0
x_test=bisect(biexpfunc,x.min(),args=tuple(test_popt))
#we calculate the deviation from the expected y-value
tol_test,= np.abs(y_test - biexpfunc(np.asarray([x_test]),*popt))
#and mark the determined point in the graph
ax.axhline(y_test,ls="--",color="grey")
ax.axvline(x_test,color="grey")
ax.plot(x_test,y_test,c="tab:orange",marker="o",markersize=15,alpha=0.5)
ax.annotate(f"X: {x_test:.2f},Y: {y_test:.2f}\ntol: {tol_test:.4f}",xy=(x_test,y_test),xytext=(50,50),textcoords="offset points",arrowprops=dict(facecolor="tab:orange",shrink=0.05),)
ax.legend(title="root finding: bisect")
plt.show()
另一种确定更复杂函数根的方法是,令人惊讶的是,root
。脚本基本相同,只是root-routine略有不同,比如我们可以选择寻根方式:
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit,root
def biexpfunc(x,label="fit")
#y-value for which we want to determine the x-value(s)
y_test=55
test_popt = popt.copy()
test_popt[-1] -= y_test
#calculate corresponding x-value with root finding
r=root(biexpfunc,x.mean(),args=tuple(test_popt),method="lm")
x_test,= r.x
tol_test,= np.abs(y_test - biexpfunc(r.x,*popt))
#mark point in graph
ax.axhline(y_test,shrink=0.05))
ax.legend(title="root finding: lm")
plt.show()
在这种情况下,图表看起来相同。这不一定适用于每个功能;就像曲线拟合一样,正确的方法可以显着改善结果。