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)
Now, Heaviside function must be generated. To give you an overview of this function, consider the half-maximum convention of Heaviside function:

In Python, this is equivalent to: def f(x): return 0.5 * (np.sign(x) + 1)
A sample plot would be:
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)

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)

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:

or in Python: (0.5 * (np.sign(x-5) + 1) * 100 = 50 * (np.sign(x-5) + 1)
Combining the plots yields (where Fit
represents the above fitting function)

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
.
Proceed to estimation:
popt,pcov=curve_fit(f,xobs,yobs,bounds=([49,4.75,0],[50,5,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.
Here are MY results for popt
and pcov
:

pcov
represents the estimated covariance of popt. The diagonals provide the variance of the parameter estimate [Source].
Results show that the parameter estimates pcov
are near the theoretical values.
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
Finally, note that the estimates of popt
and pcov
vary.
pythonscipy
jobs
-->yobs
?[40.,0.,100.]
represent ?