如何在python中拟合阶跃函数 [英] how to fit a step function in python
问题描述
我有一个关于使用曲线拟合等 scipy 例程拟合阶跃函数的问题.我无法将其矢量化，例如:
I have a question about fitting a step function using scipy routines like curve_fit. I have trouble making it vectorized, for example:
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
xobs=np.linspace(0,10,100)
yl=np.random.rand(50); yr=np.random.rand(50)+100
yobs=np.concatenate((yl,yr),axis=0)
def model(x,rf,T1,T2):
#1: x=np.vectorize(x)
if x<rf:
ret= T1
else:
ret= T2
return ret
#2: model=np.vectorize(model)
popt, pcov = curve_fit(model, xobs, yobs, [40.,0.,100.])
它说
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
如果我添加 #1 或 #2，它会运行但并不真正适合数据:
If I add #1 or #2 it runs but doesn't really fit the data:
OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning)
[ 40. 50.51182064 50.51182064] [[ inf inf inf]
[ inf inf inf]
[ inf inf inf]]
有人知道怎么解决吗?谢谢
Anybody know how to fix that? THX
推荐答案
这是我所做的.我保留了 xobs 和 yobs:
Here's what I did. I retained xobs and yobs:
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
xobs=np.linspace(0,10,100)
yl=np.random.rand(50); yr=np.random.rand(50)+100
yobs=np.concatenate((yl,yr),axis=0)
现在，必须生成 Heaviside 函数.为了给你这个函数的概述，考虑 Heaviside 函数的半最大值约定:
Now, Heaviside function must be generated. To give you an overview of this function, consider the half-maximum convention of Heaviside function:
在 Python 中，这相当于:def f(x): return 0.5 * (np.sign(x) + 1)
In Python, this is equivalent to: def f(x): return 0.5 * (np.sign(x) + 1)
示例图为:
xval = sorted(np.concatenate([np.linspace(-5,5,100),[0]])) # includes x = 0
yval = f(xval)
plt.plot(xval,yval,'ko-')
plt.ylim(-0.1,1.1)
plt.xlabel('x',size=18)
plt.ylabel('H(x)',size=20)
现在，绘制 xobs 和 yobs 给出:
Now, plotting xobs and yobs gives:
plt.plot(xobs,yobs,'ko-')
plt.ylim(-10,110)
plt.xlabel('xobs',size=18)
plt.ylabel('yobs',size=20)
请注意，比较两个图，第二个图移动了 5 个单位，最大值从 1.0 增加到 100.我推断第二个图的函数可以表示如下:
Notice that comparing the two figures, the second plot is shifted by 5 units and the maximum increases from 1.0 to 100. I infer that the function for the second plot can be represented as follows:
或在 Python 中:(0.5 * (np.sign(x-5) + 1) * 100 = 50 * (np.sign(x-5) + 1)
or in Python: (0.5 * (np.sign(x-5) + 1) * 100 = 50 * (np.sign(x-5) + 1)
组合图产生(其中Fit代表上述拟合函数)
Combining the plots yields (where Fit represents the above fitting function)
情节证实了我的猜测是正确的.现在，假设您不知道这个正确的拟合函数是如何产生的，则创建一个广义拟合函数: def f(x,a,b,c): return a * (np.sign(xb) +c)，理论上a = 50，b = 5，c = 1.
The plot confirms that my guess is correct. Now, assuming that YOU DO NOT KNOW how did this correct fitting function come about, a generalized fitting function is created: def f(x,a,b,c): return a * (np.sign(x-b) + c), where theoretically, a = 50, b = 5, and c = 1.
继续估算:
popt,pcov=curve_fit(f,xobs,yobs,bounds=([49,4.75,0],[50,5,2])).
现在，bounds = ([每个参数的下限(a,b,c)],[每个参数的上限]).从技术上讲，这意味着 49 <a <50, 4.75 ＜b <5，且0＜c 2.
Now, bounds = ([lower bound of each parameter (a,b,c)],[upper bound of each parameter]). Technically, this means that 49 < a < 50, 4.75 < b < 5, and 0 < c < 2.
以下是我对 popt 和 pcov 的结果:
Here are MY results for popt and pcov:
pcov 表示估计的 popt 协方差.对角线提供了参数估计的方差 [来源].
pcov represents the estimated covariance of popt. The diagonals provide the variance of the parameter estimate [Source].
结果表明参数估计pcov接近理论值.
Results show that the parameter estimates pcov are near the theoretical values.
基本上，一个广义的 Heaviside 函数可以表示为:a * (np.sign(x-b) + c)
Basically, a generalized Heaviside function can be represented by: a * (np.sign(x-b) + c)
这是将生成参数估计和相应协方差的代码:
Here is the code that will generate parameter estimates and the corresponding covariances:
import numpy as np
from scipy.optimize import curve_fit
xobs = np.linspace(0,10,100)
yl = np.random.rand(50); yr=np.random.rand(50)+100
yobs = np.concatenate((yl,yr),axis=0)
def f(x,a,b,c): return a * (np.sign(x-b) + c) # Heaviside fitting function
popt, pcov = curve_fit(f,xobs,yobs,bounds=([49,4.75,0],[50,5,2]))
print 'popt = %s' % popt
print 'pcov = \n %s' % pcov
最后，请注意 popt 和 pcov 的估计值不同.
Finally, note that the estimates of popt and pcov vary.
这篇关于如何在python中拟合阶跃函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
